Added tnlUmfpackWrapper class (0a74e5a1) · Commits · TNL / tnl-dev

src/matrices/tnlCSRMatrix.h

+9 −0

Original line number	Diff line number	Diff line
		@@ -21,6 +21,11 @@
		#include <matrices/tnlSparseMatrix.h>
		#include <core/vectors/tnlVector.h>

		#ifdef HAVE_UMFPACK
		template< typename Matrix, typename Preconditioner >
		class tnlUmfpackWrapper;
		#endif

		template< typename Real >
		class tnlCusparseCSRMatrix;

		@@ -209,6 +214,10 @@ class tnlCSRMatrix : public tnlSparseMatrix< Real, Device, Index >
		typedef tnlCSRMatrixDeviceDependentCode< DeviceType > DeviceDependentCode;
		friend class tnlCSRMatrixDeviceDependentCode< DeviceType >;
		friend class tnlCusparseCSRMatrix< RealType >;
		#ifdef HAVE_UMFPACK
		template< typename Matrix, typename Preconditioner >
		friend class tnlUmfpackWrapper;
		#endif

		};

src/solvers/linear/tnlUmfpackWrapper.h

0 → 100644

+122 −0

Original line number	Diff line number	Diff line
		#pragma once

		#ifdef HAVE_UMFPACK

		#include <umfpack.h>

		#include <core/tnlObject.h>
		#include <config/tnlConfigDescription.h>
		#include <matrices/tnlCSRMatrix.h>
		#include <solvers/preconditioners/tnlDummyPreconditioner.h>
		#include <solvers/tnlIterativeSolver.h>
		#include <solvers/linear/tnlLinearResidueGetter.h>


		template< typename Matrix >
		struct is_csr_matrix
		{
		static const bool value = false;
		};

		template< typename Real, typename Device, typename Index >
		struct is_csr_matrix< tnlCSRMatrix< Real, Device, Index > >
		{
		static const bool value = true;
		};


		template< typename Matrix,
		typename Preconditioner = tnlDummyPreconditioner< typename Matrix :: RealType,
		typename Matrix :: DeviceType,
		typename Matrix :: IndexType> >
		class tnlUmfpackWrapper
		: public tnlObject,
		// just to ensure the same interface as other linear solvers
		public tnlIterativeSolver< typename Matrix::RealType,
		typename Matrix::IndexType >
		{
		public:
		typedef typename Matrix :: RealType RealType;
		typedef typename Matrix :: IndexType IndexType;
		typedef typename Matrix :: DeviceType DeviceType;
		typedef Matrix MatrixType;
		typedef Preconditioner PreconditionerType;

		tnlUmfpackWrapper()
		{
		if( ! is_csr_matrix< Matrix >::value )
		cerr << "The tnlUmfpackWrapper solver is available only for CSR matrices." << endl;
		if( std::is_same< typename Matrix::DeviceType, tnlCuda >::value )
		cerr << "The tnlUmfpackWrapper solver is not available on CUDA." << endl;
		if( ! std::is_same< RealType, double >::value )
		cerr << "The tnlUmfpackWrapper solver is available only for double precision." << endl;
		if( ! std::is_same< IndexType, int >::value )
		cerr << "The tnlUmfpackWrapper solver is available only for 'int' index type." << endl;
		}

		static void configSetup( tnlConfigDescription& config,
		const tnlString& prefix = "" )
		{};

		bool setup( const tnlParameterContainer& parameters,
		const tnlString& prefix = "" )
		{
		return false;
		};

		void setMatrix( const MatrixType& matrix )
		{};

		void setPreconditioner( const Preconditioner& preconditioner )
		{};

		template< typename Vector,
		typename ResidueGetter = tnlLinearResidueGetter< MatrixType, Vector > >
		bool solve( const Vector& b, Vector& x )
		{
		return false;
		};

		};


		template< typename Preconditioner >
		class tnlUmfpackWrapper< tnlCSRMatrix< double, tnlHost, int >, Preconditioner >
		: public tnlObject,
		// just to ensure the same interface as other linear solvers
		public tnlIterativeSolver< double, int >
		{
		public:
		typedef double RealType;
		typedef int IndexType;
		typedef tnlHost DeviceType;
		typedef tnlCSRMatrix< double, tnlHost, int > MatrixType;
		typedef Preconditioner PreconditionerType;

		tnlUmfpackWrapper();

		tnlString getType() const;

		static void configSetup( tnlConfigDescription& config,
		const tnlString& prefix = "" );

		bool setup( const tnlParameterContainer& parameters,
		const tnlString& prefix = "" );

		void setMatrix( const MatrixType& matrix );

		void setPreconditioner( const Preconditioner& preconditioner );

		template< typename Vector,
		typename ResidueGetter = tnlLinearResidueGetter< MatrixType, Vector > >
		bool solve( const Vector& b, Vector& x );

		protected:
		const MatrixType* matrix;

		const PreconditionerType* preconditioner;
		};

		#include "tnlUmfpackWrapper_impl.h"

		#endif

src/solvers/linear/tnlUmfpackWrapper_impl.h

0 → 100644

+129 −0

Original line number	Diff line number	Diff line
		#pragma once

		#ifdef HAVE_UMFPACK

		#include "tnlUmfpackWrapper.h"

		template< typename Preconditioner >
		tnlUmfpackWrapper< tnlCSRMatrix< double, tnlHost, int >, Preconditioner >::
		tnlUmfpackWrapper()
		{}

		template< typename Preconditioner >
		void
		tnlUmfpackWrapper< tnlCSRMatrix< double, tnlHost, int >, Preconditioner >::
		configSetup( tnlConfigDescription& config,
		const tnlString& prefix )
		{
		}

		template< typename Preconditioner >
		bool
		tnlUmfpackWrapper< tnlCSRMatrix< double, tnlHost, int >, Preconditioner >::
		setup( const tnlParameterContainer& parameters,
		const tnlString& prefix )
		{
		return true;
		}

		template< typename Preconditioner >
		void tnlUmfpackWrapper< tnlCSRMatrix< double, tnlHost, int >, Preconditioner >::
		setMatrix( const MatrixType& matrix )
		{
		this -> matrix = &matrix;
		}

		template< typename Preconditioner >
		void tnlUmfpackWrapper< tnlCSRMatrix< double, tnlHost, int >, Preconditioner >::
		setPreconditioner( const Preconditioner& preconditioner )
		{
		this -> preconditioner = &preconditioner;
		}


		template< typename Preconditioner >
		template< typename Vector, typename ResidueGetter >
		bool tnlUmfpackWrapper< tnlCSRMatrix< double, tnlHost, int >, Preconditioner >::
		solve( const Vector& b,
		Vector& x )
		{
		tnlAssert( matrix->getRows() == matrix->getColumns(), );
		tnlAssert( matrix->getColumns() == x.getSize() && matrix->getColumns() == b.getSize(), );

		const IndexType size = matrix -> getRows();

		this->resetIterations();
		this->setResidue( this -> getConvergenceResidue() + 1.0 );

		RealType bNorm = b. lpNorm( ( RealType ) 2.0 );

		// UMFPACK objects
		void* Symbolic = nullptr;
		void* Numeric = nullptr;

		int status = UMFPACK_OK;
		double Control[ UMFPACK_CONTROL ];
		double Info[ UMFPACK_INFO ];
		// Control[ UMFPACK_PRL ] = 2;

		// umfpack expects Compressed Sparse Column format, we have Compressed Sparse Row
		// so we need to solve A^T * x = rhs
		int system_type = UMFPACK_Aat;

		// symbolic reordering of the sparse matrix
		status = umfpack_di_symbolic( size, size,
		matrix->rowPointers.getData(),
		matrix->columnIndexes.getData(),
		matrix->values.getData(),
		&Symbolic, Control, Info );
		if( status != UMFPACK_OK ) {
		cerr << "error: symbolic reordering failed" << endl;
		umfpack_di_report_status( Control, status );
		// umfpack_di_report_control( Control );
		// umfpack_di_report_info( Control, Info );
		goto finished;
		}

		// numeric factorization
		status = umfpack_di_numeric( matrix->rowPointers.getData(),
		matrix->columnIndexes.getData(),
		matrix->values.getData(),
		Symbolic, &Numeric, Control, Info );
		if( status != UMFPACK_OK ) {
		cerr << "error: numeric factorization failed" << endl;
		umfpack_di_report_status( Control, status );
		// umfpack_di_report_control( Control );
		// umfpack_di_report_info( Control, Info );
		goto finished;
		}

		// solve with specified right-hand-side
		status = umfpack_di_solve( system_type,
		matrix->rowPointers.getData(),
		matrix->columnIndexes.getData(),
		matrix->values.getData(),
		x.getData(),
		b.getData(),
		Numeric, Control, Info );
		if( status != UMFPACK_OK ) {
		cerr << "error: umfpack_di_solve failed" << endl;
		umfpack_di_report_status( Control, status );
		// umfpack_di_report_control( Control );
		// umfpack_di_report_info( Control, Info );
		goto finished;
		}

		// umfpack_di_report_info( Control, Info );

		finished:
		if( Symbolic )
		umfpack_di_free_symbolic( &Symbolic );
		if( Numeric )
		umfpack_di_free_numeric( &Numeric );

		this->setResidue( ResidueGetter::getResidue( *matrix, x, b, bNorm ) );
		this->refreshSolverMonitor( true );
		return status == UMFPACK_OK;
		};

		#endif

src/solvers/tnlBuildConfigTags.h

+12 −0

Original line number	Diff line number	Diff line
		@@ -26,6 +26,7 @@
		#include <solvers/linear/krylov/tnlBICGStabSolver.h>
		#include <solvers/linear/krylov/tnlGMRESSolver.h>
		#include <solvers/linear/krylov/tnlTFQMRSolver.h>
		#include <solvers/linear/tnlUmfpackWrapper.h>
		#include <solvers/preconditioners/tnlDummyPreconditioner.h>

		class tnlDefaultBuildMeshConfig{};
		@@ -151,6 +152,17 @@ public:
		using Template = tnlTFQMRSolver< Matrix, Preconditioner >;
		};

		#ifdef HAVE_UMFPACK
		class tnlSemiImplicitUmfpackSolverTag
		{
		public:
		template< typename Matrix, typename Preconditioner = tnlDummyPreconditioner< typename Matrix::RealType,
		typename Matrix::DeviceType,
		typename Matrix::IndexType > >
		using Template = tnlUmfpackWrapper< Matrix, Preconditioner >;
		};
		#endif

		template< typename MeshConfig, typename SemiImplicitSolver > struct tnlMeshConfigSemiImplicitSolver{ enum { enabled = true }; };

		#endif /* TNLMeshConfigS_H_ */

src/solvers/tnlSolverConfig_impl.h

+4 −0

Original line number	Diff line number	Diff line
		@@ -118,6 +118,10 @@ bool tnlSolverConfig< MeshConfig, ProblemConfig >::configSetup( tnlConfigDescrip
		config.addEntryEnum( "tfqmr" );
		if( tnlMeshConfigSemiImplicitSolver< MeshConfig, tnlSemiImplicitSORSolverTag >::enabled )
		config.addEntryEnum( "sor" );
		#ifdef HAVE_UMFPACK
		if( tnlMeshConfigSemiImplicitSolver< MeshConfig, tnlSemiImplicitUmfpackSolverTag >::enabled )
		config.addEntryEnum( "umfpack" );
		#endif
		}
		config.addEntry< tnlString >( "preconditioner", "The preconditioner for the discrete solver:", "none" );
		config.addEntryEnum( "none" );