Skip to content
Snippets Groups Projects
Commit a9a10f23 authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

Refactored Jacobi solver

parent cb9918bb
No related branches found
No related tags found
No related merge requests found
......@@ -10,8 +10,9 @@ SET( headers BICGStab.h
CWYGMRES_impl.h
GMRES.h
GMRES_impl.h
Jacobi.h
LinearResidueGetter.h
LinearResidueGetter_impl.h
LinearResidueGetter_impl.h
SOR.h
SOR_impl.h
TFQMR.h
......
......@@ -6,34 +6,24 @@
email : tomas.oberhuber@fjfi.cvut.cz
***************************************************************************/
/***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/
#pragma once
#pragma once
#include <math.h>
#include <TNL/Object.h>
#include <TNL/Solvers/Preconditioners/Dummy.h>
#include <TNL/Solvers/IterativeSolver.h>
#include <TNL/Solvers/Linear/Preconditioners/Dummy.h>
#include <TNL/Solvers/Linear/LinearResidueGetter.h>
namespace TNL {
namespace Solvers {
namespace Linear {
namespace Solvers {
namespace Linear {
template< typename Matrix,
typename Preconditioner = DummyPreconditioner< typename Matrix :: RealType,
typename Matrix :: DeviceType,
typename Matrix :: IndexType> >
typename Preconditioner = Preconditioners::Dummy< typename Matrix::RealType,
typename Matrix::DeviceType,
typename Matrix::IndexType> >
class Jacobi : public Object,
public IterativeSolver< typename Matrix :: RealType,
typename Matrix :: IndexType >
public IterativeSolver< typename Matrix::RealType,
typename Matrix::IndexType >
{
public:
......@@ -44,14 +34,14 @@ class Jacobi : public Object,
typedef Preconditioner PreconditionerType;
Jacobi():omega(0){}
Jacobi() : omega(0) {}
String getType() const
{
return String( "Jacobi< " ) + this -> matrix -> getType() + ", " + this -> preconditioner -> getType() + " >";
return String( "Jacobi< " ) + this->matrix->getType() + ", " + this->preconditioner->getType() + " >";
}
static void configSetup( tnlConfigDescription& config,
static void configSetup( Config::ConfigDescription& config,
const String& prefix = "" )
{
config.addEntry< double >( prefix + "jacobi-omega", "Relaxation parameter of the weighted/damped Jacobi method.", 1.0 );
......@@ -59,21 +49,21 @@ class Jacobi : public Object,
bool setup( const Config::ParameterContainer& parameters,
const String& prefix = "" )
{
IterativeSolver< RealType, IndexType >::setup( parameters, prefix );
this->setOmega( parameters.getParameter< double >( prefix + "jacobi-omega" ) );
if( this->omega <= 0.0 || this->omega > 2.0 )
{
cerr << "Warning: The Jacobi method parameter omega is out of interval (0,2). The value is " << this->omega << " the method will not converge." << endl;
}
return true;
{
IterativeSolver< RealType, IndexType >::setup( parameters, prefix );
this->setOmega( parameters.getParameter< double >( prefix + "jacobi-omega" ) );
if( this->omega <= 0.0 || this->omega > 2.0 )
{
std::cerr << "Warning: The Jacobi method parameter omega is out of interval (0,2). The value is " << this->omega << " the method will not converge." << std::endl;
}
return true;
}
void setOmega( const RealType& omega )
{
this->omega = omega;
}
const RealType& getOmega() const
{
return omega;
......@@ -81,24 +71,24 @@ class Jacobi : public Object,
void setMatrix( const MatrixType& matrix )
{
this -> matrix = &matrix;
this->matrix = &matrix;
}
void setPreconditioner( const Preconditioner& preconditioner )
{
this -> preconditioner = &preconditioner;
this->preconditioner = &preconditioner;
}
template< typename Vector,
typename ResidueGetter = LinearResidueGetter< Matrix, Vector > >
bool solve( const Vector& b, Vector& x, Vector& aux)
{
const IndexType size = matrix -> getRows();
const IndexType size = matrix->getRows();
this -> resetIterations();
this -> setResidue( this -> getConvergenceResidue() + 1.0 );
this->resetIterations();
this->setResidue( this->getConvergenceResidue() + 1.0 );
RealType bNorm = b. lpNorm( ( RealType ) 2.0 );
RealType bNorm = b.lpNorm( ( RealType ) 2.0 );
aux = x;
while( this->nextIteration() )
{
......@@ -106,10 +96,10 @@ class Jacobi : public Object,
matrix->performJacobiIteration( b, row, x, aux, this->getOmega() );
for( IndexType row = 0; row < size; row ++ )
matrix->performJacobiIteration( b, row, aux, x, this->getOmega() );
this -> setResidue( ResidueGetter :: getResidue( *matrix, x, b, bNorm ) );
this->setResidue( ResidueGetter::getResidue( *matrix, x, b, bNorm ) );
}
this -> setResidue( ResidueGetter :: getResidue( *matrix, x, b, bNorm ) );
this -> refreshSolverMonitor();
this->setResidue( ResidueGetter::getResidue( *matrix, x, b, bNorm ) );
this->refreshSolverMonitor();
return this->checkConvergence();
}
......@@ -118,9 +108,8 @@ class Jacobi : public Object,
RealType omega;
const MatrixType* matrix;
const PreconditionerType* preconditioner;
};
} // namespace Linear
} // namespace Solvers
} // namespace Linear
} // namespace Solvers
} // namespace TNL
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment