From: young <young@0785d39b-7218-0410-832d-ea1e28bc413d> Date: Mon, 10 May 2010 06:56:39 +0000 (+0000) Subject: A class SparseDirectMUMPS that gives a interface to the sparse direct solver in MUMPS X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6bcdcf4fc21c3d4e5535cdc0c299d3ba7b53f70f;p=dealii-svn.git 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 --- 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 <lac/sparse_matrix.h> #include <lac/block_sparse_matrix.h> -DEAL_II_NAMESPACE_OPEN +#ifdef DEAL_II_USE_MUMPS +# include <base/utilities.h> +# include <dmumps_c.h> +#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 <a href="http://mumps.enseeiht.fr">MUMPS</a>. 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 <a href="http://mumps.enseeiht.fr">MUMPS</a> exists on your + * system and was detected during configuration of + * <code>deal.II</code>. + * + * <h4>Instantiations</h4> + * + * There are instantiations of this class for SparseMatrix<double>, + * SparseMatrix<float>, BlockSparseMatrix<double>, and + * BlockSparseMatrix<float>. + * + * @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 + * <tt>initialize ()</tt> 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 + * <tt>matrix</tt> and right-hand side + * <tt>vector</tt> to the solver. + */ + template <class Matrix> + void initialize (const SparseMatrix<double>& matrix, + const Vector<double> & vector); + + /** + * A function in which the linear system is + * solved and the solution vector is copied + * into the given <tt>vector</tt>. + */ + void solve (Vector<double>& 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 <list> #include <typeinfo> - // 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 <class Matrix> +void SparseDirectMUMPS::initialize (const SparseMatrix<double>& matrix, + const Vector<double> & 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<double>::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<double>& 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<n; ++i) + vector(i) = rhs[i]; + + delete[] a; + delete[] irn; + delete[] jcn; + delete[] rhs; + } +} +#endif // DEAL_II_USE_MUMPS // explicit instantiations for SparseMatrixMA27 template @@ -1911,4 +2009,13 @@ InstantiateUMFPACK(SparseMatrix<float>); InstantiateUMFPACK(BlockSparseMatrix<double>); InstantiateUMFPACK(BlockSparseMatrix<float>); +#ifdef DEAL_II_USE_MUMPS +// explicit instantiations for SparseDirectMUMPS +template <class Matrix> +void SparseDirectMUMPS::initialize (const SparseMatrix<double>& matrix, + const Vector<double> & vector); + +void SparseDirectMUMPS::solve (Vector<double>& vector); +#endif + DEAL_II_NAMESPACE_CLOSE