From 6bcdcf4fc21c3d4e5535cdc0c299d3ba7b53f70f Mon Sep 17 00:00:00 2001 From: young Date: Mon, 10 May 2010 06:56:39 +0000 Subject: [PATCH] A class SparseDirectMUMPS that gives a interface to the sparse direct solver in MUMPS git-svn-id: https://svn.dealii.org/trunk@21110 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/sparse_direct.h | 85 +++++++++++++++++- deal.II/lac/source/sparse_direct.cc | 109 +++++++++++++++++++++++- 2 files changed, 190 insertions(+), 4 deletions(-) diff --git a/deal.II/lac/include/lac/sparse_direct.h b/deal.II/lac/include/lac/sparse_direct.h index 86be9770ae..f3792bbc56 100644 --- a/deal.II/lac/include/lac/sparse_direct.h +++ b/deal.II/lac/include/lac/sparse_direct.h @@ -23,8 +23,12 @@ #include #include -DEAL_II_NAMESPACE_OPEN +#ifdef DEAL_II_USE_MUMPS +# include +# include +#endif +DEAL_II_NAMESPACE_OPEN /** * This class provides an interface to the sparse direct solver MA27 @@ -357,6 +361,10 @@ class SparseDirectMA27 : public Subscriptor /** * Exception */ + DeclException0 (ExcInitializeNotCalled); + /** + * Exception + */ DeclException0 (ExcFactorizeNotCalled); /** * Exception @@ -1260,8 +1268,79 @@ class SparseDirectUMFPACK : public Subscriptor }; - +#ifdef DEAL_II_USE_MUMPS +/** + * This class provides an interface to the parallel sparse direct + * solver MUMPS. MUMPS is + * direct method based on a multifrontal approach, which performs a + * direct LU factorization. The matrix coming in may have either + * symmetric or nonsymmetric sparsity pattern. + * + * @note This class is useable if and only if a working installation + * of MUMPS exists on your + * system and was detected during configuration of + * deal.II. + * + *

Instantiations

+ * + * There are instantiations of this class for SparseMatrix, + * SparseMatrix, BlockSparseMatrix, and + * BlockSparseMatrix. + * + * @author Markus Buerg, 2010 + */ +class SparseDirectMUMPS +{ + private: + + DMUMPS_STRUC_C id; + + double *a; + double *rhs; + + unsigned int *irn; + unsigned int *jcn; + unsigned int n; + unsigned int nz; + + /** + * Flags storing whether the function + * initialize () has already been + * called. + */ + bool initialize_called; + + public: + + /** + * Constructor + */ + SparseDirectMUMPS (); + + /** + * Destructor + */ + ~SparseDirectMUMPS (); + + /** + * This function initializes a MUMPS instance + * and hands over the system's matrix + * matrix and right-hand side + * vector to the solver. + */ + template + void initialize (const SparseMatrix& matrix, + const Vector & vector); + + /** + * A function in which the linear system is + * solved and the solution vector is copied + * into the given vector. + */ + void solve (Vector& vector); +}; +#endif // DEAL_II_USE_MUMPS DEAL_II_NAMESPACE_CLOSE -#endif +#endif // __deal2__sparse_direct_h diff --git a/deal.II/lac/source/sparse_direct.cc b/deal.II/lac/source/sparse_direct.cc index 567052469a..a8c98b9724 100644 --- a/deal.II/lac/source/sparse_direct.cc +++ b/deal.II/lac/source/sparse_direct.cc @@ -23,7 +23,6 @@ #include #include - // this is a weird hack: on newer linux systems, some system headers // include /usr/include/linux/compiler.h which explicitly checks which // gcc is in use. in that file is also a comment that explains that @@ -1878,6 +1877,105 @@ SparseDirectUMFPACK::Tvmult_add ( Assert(false, ExcNotImplemented()); } +#ifdef DEAL_II_USE_MUMPS +SparseDirectMUMPS::SparseDirectMUMPS () +{} + +SparseDirectMUMPS::~SparseDirectMUMPS () +{} + +template +void SparseDirectMUMPS::initialize (const SparseMatrix& matrix, + const Vector & vector) +{ + + // Initialize MUMPS instance: + id.job = -1; + id.par = 1; + id.sym = 0; + + // Use MPI_COMM_WORLD as communicator + id.comm_fortran = -987654; + dmumps_c (&id); + + // Hand over matrix and right-hand side + if (Utilities::System::get_this_mpi_process (MPI_COMM_WORLD) == 0) + { + + // Objects denoting a MUMPS data structure: + // + // Set number of unknowns + n = vector.size (); + + // number of nonzero elements in matrix + nz = matrix.n_actually_nonzero_elements (); + + // matrix + a = new double[nz]; + + // vector indices pointing to the first and + // last elements of the vector respectively + irn = new unsigned int[nz]; + jcn = new unsigned int[nz]; + + unsigned int index = 0; + + for (SparseMatrix::const_iterator ptr = matrix.begin (); + ptr != matrix.end (); ++ptr) + if (std::abs (ptr->value ()) > 0.0) + { + a[index] = ptr->value (); + irn[index] = ptr->row () + 1; + jcn[index] = ptr->column () + 1; + ++index; + } + + rhs = new double[n]; + + for (unsigned int i = 0; i < n; ++i) + rhs[i] = vector (i); + + id.n = n; + id.nz = nz; + id.irn = irn; + id.jcn = jcn; + id.a = a; + id.rhs = rhs; + } + + // No outputs + id.icntl[0] = -1; + id.icntl[1] = -1; + id.icntl[2] = -1; + id.icntl[3] = 0; + + // Exit by setting this flag: + initialize_called = true; +} + +void SparseDirectMUMPS::solve (Vector& vector) +{ + + // Start solver + id.job = 6; + dmumps_c (&id); + + id.job = -2; + dmumps_c (&id); + + // Copy solution into the given vector + if (Utilities::System::get_this_mpi_process (MPI_COMM_WORLD) == 0) + { + for (unsigned int i=0; i); InstantiateUMFPACK(BlockSparseMatrix); InstantiateUMFPACK(BlockSparseMatrix); +#ifdef DEAL_II_USE_MUMPS +// explicit instantiations for SparseDirectMUMPS +template +void SparseDirectMUMPS::initialize (const SparseMatrix& matrix, + const Vector & vector); + +void SparseDirectMUMPS::solve (Vector& vector); +#endif + DEAL_II_NAMESPACE_CLOSE -- 2.39.5