Skip to content
Snippets Groups Projects
Commit b9f067c5 authored by Tomáš Oberhuber's avatar Tomáš Oberhuber
Browse files

Linear Jacobi solver refactoring.

parent f52140d2
No related branches found
No related tags found
No related merge requests found
/***************************************************************************
tnlSORSolver.h - description
Jacobi.h - description
-------------------
begin : 2007/07/30
begin : Jul 30, 2007
copyright : (C) 2007 by Tomas Oberhuber
email : tomas.oberhuber@fjfi.cvut.cz
***************************************************************************/
......@@ -15,47 +15,52 @@
* *
***************************************************************************/
#ifndef tnlJacobiSolverH
#define tnlJacobiSolverH
#pragma once
#include <math.h>
#include <core/tnlObject.h>
#include <solvers/preconditioners/tnlDummyPreconditioner.h>
#include <solvers/tnlIterativeSolver.h>
#include <solvers/linear/tnlLinearResidueGetter.h>
#include <TNL/Object.h>
#include <TNL/Solvers/Preconditioners/Dummy.h>
#include <TNL/Solvers/IterativeSolver.h>
#include <TNL/Solvers/Linear/LinearResidueGetter.h>
namespace TNL {
namespace Solvers {
namespace Linear {
template< typename Matrix,
typename Preconditioner = tnlDummyPreconditioner< typename Matrix :: RealType,
typename Matrix :: DeviceType,
typename Matrix :: IndexType> >
class tnlJacobiSolver : public tnlObject,
public tnlIterativeSolver< typename Matrix :: RealType,
typename Matrix :: IndexType >
typename Preconditioner = DummyPreconditioner< typename Matrix :: RealType,
typename Matrix :: DeviceType,
typename Matrix :: IndexType> >
class Jacobi : public Object,
public IterativeSolver< typename Matrix :: RealType,
typename Matrix :: IndexType >
{
public:
typedef typename Matrix :: RealType RealType;
typedef typename Matrix :: IndexType IndexType;
typedef typename Matrix :: DeviceType DeviceType;
typedef typename Matrix::RealType RealType;
typedef typename Matrix::IndexType IndexType;
typedef typename Matrix::DeviceType DeviceType;
typedef Matrix MatrixType;
typedef Preconditioner PreconditionerType;
tnlJacobiSolver():omega(0){}
Jacobi():omega(0){}
tnlString getType() const
String getType() const
{
return tnlString( "tnlJacobiSolver< " ) + this -> matrix -> getType() + ", " + this -> preconditioner -> getType() + " >";
return String( "Jacobi< " ) + this -> matrix -> getType() + ", " + this -> preconditioner -> getType() + " >";
}
static void configSetup( tnlConfigDescription& config,
const tnlString& prefix = "" ) {config.addEntry< double >( prefix + "jacobi-omega", "Relaxation parameter of the weighted/damped Jacobi method.", 1.0 );}
const String& prefix = "" )
{
config.addEntry< double >( prefix + "jacobi-omega", "Relaxation parameter of the weighted/damped Jacobi method.", 1.0 );
}
bool setup( const tnlParameterContainer& parameters,
const tnlString& prefix = "" )
bool setup( const Config::ParameterContainer& parameters,
const String& prefix = "" )
{
tnlIterativeSolver< RealType, IndexType >::setup( parameters, prefix );
IterativeSolver< RealType, IndexType >::setup( parameters, prefix );
this->setOmega( parameters.getParameter< double >( prefix + "jacobi-omega" ) );
if( this->omega <= 0.0 || this->omega > 2.0 )
{
......@@ -64,49 +69,58 @@ class tnlJacobiSolver : public tnlObject,
return true;
}
void setOmega( const RealType& omega ) {this->omega = omega;}
const RealType& getOmega() const {return omega;}
void setOmega( const RealType& omega )
{
this->omega = omega;
}
const RealType& getOmega() const
{
return omega;
}
void setMatrix( const MatrixType& matrix ){this -> matrix = &matrix;}
void setPreconditioner( const Preconditioner& preconditioner ){this -> preconditioner = &preconditioner;}
void setMatrix( const MatrixType& matrix )
{
this -> matrix = &matrix;
}
void setPreconditioner( const Preconditioner& preconditioner )
{
this -> preconditioner = &preconditioner;
}
#ifdef HAVE_NOT_CXX11
template< typename Vector,
typename ResidueGetter >
bool solve( const Vector& b, Vector& x )
#else
template< typename Vector,
typename ResidueGetter = tnlLinearResidueGetter< Matrix, Vector > >
typename ResidueGetter = LinearResidueGetter< Matrix, Vector > >
bool solve( const Vector& b, Vector& x, Vector& aux)
#endif
{
const IndexType size = matrix -> getRows();
this -> resetIterations();
this -> setResidue( this -> getConvergenceResidue() + 1.0 );
RealType bNorm = b. lpNorm( ( RealType ) 2.0 );
aux = x;
while( this->nextIteration() )
{
for( IndexType row = 0; row < size; row ++ )
matrix->performJacobiIteration( b, row, x, aux, this->getOmega() );
for( IndexType row = 0; row < size; row ++ )
matrix->performJacobiIteration( b, row, aux, x, this->getOmega() );
const IndexType size = matrix -> getRows();
this -> resetIterations();
this -> setResidue( this -> getConvergenceResidue() + 1.0 );
RealType bNorm = b. lpNorm( ( RealType ) 2.0 );
aux = x;
while( this->nextIteration() )
{
for( IndexType row = 0; row < size; row ++ )
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();
return this->checkConvergence();
this -> refreshSolverMonitor();
return this->checkConvergence();
}
protected:
RealType omega;
const MatrixType* matrix;
const PreconditionerType* preconditioner;
RealType omega;
const MatrixType* matrix;
const PreconditionerType* preconditioner;
};
#endif
} // 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