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 86be9770aed02dc6c05ad4cbb72b0e4a6ee30e01..f3792bbc5661133d664a2e7f1a7dd5714f33f36f 100644 (file)
--- 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 567052469a78781c9413da03313ff4e511e40d41..a8c98b9724cbbe4f6960f4dc942cca2e7a7c7fe0 100644 (file)
--- 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