diff --git a/src/TNL/Matrices/CSR.h b/src/TNL/Matrices/CSR.h index 5ad4327a6310ed1bfaa7993771a4387e27c9dae5..7c249944bfdd7ef062078f274492cb033d183eb5 100644 --- a/src/TNL/Matrices/CSR.h +++ b/src/TNL/Matrices/CSR.h @@ -199,6 +199,38 @@ class CSR : public Sparse< Real, Device, Index > int gridIdx ) const; #endif + // The following getters allow us to interface TNL with external C-like + // libraries such as UMFPACK or SuperLU, which need the raw data. + Index* getRowPointers() + { + return this->rowPointers.getData(); + } + + const Index* getRowPointers() const + { + return this->rowPointers.getData(); + } + + Index* getColumnIndexes() + { + return this->columnIndexes.getData(); + } + + const Index* getColumnIndexes() const + { + return this->columnIndexes.getData(); + } + + Real* getValues() + { + return this->values.getData(); + } + + const Real* getValues() const + { + return this->values.getData(); + } + protected: Containers::Vector< Index, Device, Index > rowPointers; @@ -210,11 +242,6 @@ class CSR : public Sparse< Real, Device, Index > typedef CSRDeviceDependentCode< DeviceType > DeviceDependentCode; friend class CSRDeviceDependentCode< DeviceType >; friend class tnlCusparseCSR< RealType >; -#ifdef HAVE_UMFPACK - template< typename Matrix, typename Preconditioner > - friend class UmfpackWrapper; -#endif - }; } // namespace Matrices diff --git a/src/TNL/Solvers/Linear/UmfpackWrapper.h b/src/TNL/Solvers/Linear/UmfpackWrapper.h index 043e46c8d72aa14caa95cc8ebdb9d93171adb6f3..e2c59d9f8ac2771c43db903aec27da5572eecb08 100644 --- a/src/TNL/Solvers/Linear/UmfpackWrapper.h +++ b/src/TNL/Solvers/Linear/UmfpackWrapper.h @@ -9,7 +9,7 @@ #include <TNL/Object.h> #include <TNL/Config/ConfigDescription.h> #include <TNL/Matrices/CSR.h> -#include <TNL/Solvers/preconditioners/Dummy.h> +#include <TNL/Solvers/Linear/Preconditioners/Dummy.h> #include <TNL/Solvers/IterativeSolver.h> #include <TNL/Solvers/Linear/LinearResidueGetter.h> @@ -25,16 +25,16 @@ struct is_csr_matrix }; template< typename Real, typename Device, typename Index > -struct is_csr_matrix< CSR< Real, Device, Index > > +struct is_csr_matrix< Matrices::CSR< Real, Device, Index > > { static const bool value = true; }; template< typename Matrix, - typename Preconditioner = Dummy< typename Matrix :: RealType, - typename Matrix :: DeviceType, - typename Matrix :: IndexType> > + typename Preconditioner = Preconditioners::Dummy< typename Matrix :: RealType, + typename Matrix :: DeviceType, + typename Matrix :: IndexType> > class UmfpackWrapper : public Object, // just to ensure the same interface as other linear solvers @@ -88,7 +88,7 @@ public: template< typename Preconditioner > -class UmfpackWrapper< CSR< double, Devices::Host, int >, Preconditioner > +class UmfpackWrapper< Matrices::CSR< double, Devices::Host, int >, Preconditioner > : public Object, // just to ensure the same interface as other linear solvers public IterativeSolver< double, int > @@ -97,7 +97,7 @@ public: typedef double RealType; typedef int IndexType; typedef Devices::Host DeviceType; - typedef CSR< double, Devices::Host, int > MatrixType; + typedef Matrices::CSR< double, Devices::Host, int > MatrixType; typedef Preconditioner PreconditionerType; typedef SharedPointer< const MatrixType, DeviceType, true > MatrixPointer; typedef SharedPointer< const PreconditionerType, DeviceType, true > PreconditionerPointer; diff --git a/src/TNL/Solvers/Linear/UmfpackWrapper_impl.h b/src/TNL/Solvers/Linear/UmfpackWrapper_impl.h index ab05e751510117a53a94a886a6fd5b731685b30f..a0a44b6e22bc6ee84f75c4f07719610b37c3d1a0 100644 --- a/src/TNL/Solvers/Linear/UmfpackWrapper_impl.h +++ b/src/TNL/Solvers/Linear/UmfpackWrapper_impl.h @@ -11,13 +11,13 @@ namespace Solvers { namespace Linear { template< typename Preconditioner > -UmfpackWrapper< CSR< double, Devices::Host, int >, Preconditioner >:: +UmfpackWrapper< Matrices::CSR< double, Devices::Host, int >, Preconditioner >:: UmfpackWrapper() {} template< typename Preconditioner > void -UmfpackWrapper< CSR< double, Devices::Host, int >, Preconditioner >:: +UmfpackWrapper< Matrices::CSR< double, Devices::Host, int >, Preconditioner >:: configSetup( Config::ConfigDescription& config, const String& prefix ) { @@ -25,7 +25,7 @@ configSetup( Config::ConfigDescription& config, template< typename Preconditioner > bool -UmfpackWrapper< CSR< double, Devices::Host, int >, Preconditioner >:: +UmfpackWrapper< Matrices::CSR< double, Devices::Host, int >, Preconditioner >:: setup( const Config::ParameterContainer& parameters, const String& prefix ) { @@ -33,14 +33,14 @@ setup( const Config::ParameterContainer& parameters, } template< typename Preconditioner > -void UmfpackWrapper< CSR< double, Devices::Host, int >, Preconditioner >:: +void UmfpackWrapper< Matrices::CSR< double, Devices::Host, int >, Preconditioner >:: setMatrix( const MatrixPointer& matrix ) { this -> matrix = matrix; } template< typename Preconditioner > -void UmfpackWrapper< CSR< double, Devices::Host, int >, Preconditioner >:: +void UmfpackWrapper< Matrices::CSR< double, Devices::Host, int >, Preconditioner >:: setPreconditioner( const PreconditionerPointer& preconditioner ) { this -> preconditioner = preconditioner; @@ -49,7 +49,7 @@ setPreconditioner( const PreconditionerPointer& preconditioner ) template< typename Preconditioner > template< typename Vector, typename ResidueGetter > -bool UmfpackWrapper< CSR< double, Devices::Host, int >, Preconditioner >:: +bool UmfpackWrapper< Matrices::CSR< double, Devices::Host, int >, Preconditioner >:: solve( const Vector& b, Vector& x ) { @@ -77,9 +77,9 @@ solve( const Vector& b, // symbolic reordering of the sparse matrix status = umfpack_di_symbolic( size, size, - matrix->rowPointers.getData(), - matrix->columnIndexes.getData(), - matrix->values.getData(), + matrix->getRowPointers(), + matrix->getColumnIndexes(), + matrix->getValues(), &Symbolic, Control, Info ); if( status != UMFPACK_OK ) { std::cerr << "error: symbolic reordering failed" << std::endl; @@ -87,9 +87,9 @@ solve( const Vector& b, } // numeric factorization - status = umfpack_di_numeric( matrix->rowPointers.getData(), - matrix->columnIndexes.getData(), - matrix->values.getData(), + status = umfpack_di_numeric( matrix->getRowPointers(), + matrix->getColumnIndexes(), + matrix->getValues(), Symbolic, &Numeric, Control, Info ); if( status != UMFPACK_OK ) { std::cerr << "error: numeric factorization failed" << std::endl; @@ -98,9 +98,9 @@ solve( const Vector& b, // solve with specified right-hand-side status = umfpack_di_solve( system_type, - matrix->rowPointers.getData(), - matrix->columnIndexes.getData(), - matrix->values.getData(), + matrix->getRowPointers(), + matrix->getColumnIndexes(), + matrix->getValues(), x.getData(), b.getData(), Numeric, Control, Info );