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