From 2452e29888d927c18bbdc177aafdf07ef636cdd4 Mon Sep 17 00:00:00 2001
From: Markus Buerg <buerg@math.tamu.edu>
Date: Thu, 11 Oct 2012 09:23:54 +0000
Subject: [PATCH] Make interface similar to the one for UMFPACK.

git-svn-id: https://svn.dealii.org/trunk@27061 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/include/deal.II/lac/sparse_direct.h |  31 ++++++
 deal.II/source/lac/sparse_direct.cc         | 101 ++++++++++++++++----
 2 files changed, 112 insertions(+), 20 deletions(-)

diff --git a/deal.II/include/deal.II/lac/sparse_direct.h b/deal.II/include/deal.II/lac/sparse_direct.h
index df517a69a9..06fef86186 100644
--- a/deal.II/include/deal.II/lac/sparse_direct.h
+++ b/deal.II/include/deal.II/lac/sparse_direct.h
@@ -1310,6 +1310,20 @@ class SparseDirectMUMPS
   unsigned int  n;
   unsigned int  nz;
 
+                        /**
+                         * This function initializes a MUMPS instance
+                         * and hands over the system's matrix
+                         * <tt>matrix</tt>.
+                         */
+  template<class Matrix>
+    void initialize_matrix (const Matrix& matrix);
+
+                        /**
+                         * Copy the computed solution into the
+                         * solution vector.
+                         */
+  void copy_solution (Vector<double>& vector);
+
                         /**
                          * Flags storing whether the function
                          * <tt>initialize ()</tt> has already been
@@ -1344,12 +1358,29 @@ class SparseDirectMUMPS
     void initialize (const Matrix& matrix,
                      const Vector<double>      & vector);
 
+                        /**
+                         * This function initializes a MUMPS instance
+                         * and computes the factorization of the
+                         * system's matrix <tt>matrix</tt>.
+                         */
+  template <class Matrix>
+    void initialize (const Matrix& matrix);
+
                         /**
                          * 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);
+  
+                        /**
+                         * A function in which the inverse of the
+                         * matrix is applied to the input vector
+                         * <tt>src</tt> and the solution is
+                         * written into the output vector
+                         * <tt>dst</tt>.
+                         */
+  void vmult (Vector<double>& dst, const Vector<double>& src);
 };
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/source/lac/sparse_direct.cc b/deal.II/source/lac/sparse_direct.cc
index 157ca2c383..98ec51a194 100644
--- a/deal.II/source/lac/sparse_direct.cc
+++ b/deal.II/source/lac/sparse_direct.cc
@@ -1957,11 +1957,21 @@ SparseDirectMUMPS::SparseDirectMUMPS ()
 {}
 
 SparseDirectMUMPS::~SparseDirectMUMPS ()
-{}
+{
+  id.job = -2;
+  dmumps_c (&id);
+  
+  // Do some cleaning
+  if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0)
+  {
+    delete[] a;
+    delete[] irn;
+    delete[] jcn;
+  }
+}
 
 template <class Matrix>
-void SparseDirectMUMPS::initialize (const Matrix& matrix,
-                                    const Vector<double>      & vector)
+void SparseDirectMUMPS::initialize_matrix (const Matrix& matrix)
 {
                          // Check we haven't been here before:
   Assert (initialize_called == false, ExcInitializeAlreadyCalled());
@@ -1978,11 +1988,10 @@ void SparseDirectMUMPS::initialize (const Matrix& matrix,
                          // Hand over matrix and right-hand side
   if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0)
     {
-
                          // Objects denoting a MUMPS data structure:
                          //
                          // Set number of unknowns
-      n   = vector.size ();
+      n   = matrix.n ();
 
                          // number of nonzero elements in matrix
       nz  = matrix.n_actually_nonzero_elements ();
@@ -2009,17 +2018,11 @@ void SparseDirectMUMPS::initialize (const Matrix& matrix,
             ++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
@@ -2032,6 +2035,47 @@ void SparseDirectMUMPS::initialize (const Matrix& matrix,
   initialize_called = true;
 }
 
+template <class Matrix>
+void SparseDirectMUMPS::initialize (const Matrix& matrix,
+                                    const Vector<double>      & vector)
+{
+                        // Hand over matrix and right-hand side
+  initialize_matrix (matrix);
+
+  if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0)
+    {
+                            // Object denoting a MUMPS data structure
+      rhs = new double[n];
+
+      for (unsigned int i = 0; i < n; ++i)
+        rhs[i] = vector (i);
+
+      id.rhs = rhs;
+    }
+}
+
+void SparseDirectMUMPS::copy_solution (Vector<double>& vector)
+{
+                       // Copy solution into the given vector
+  if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0)
+    {
+      for (unsigned int i=0; i<n; ++i)
+        vector(i) = rhs[i];
+      
+      delete[] rhs;
+    }
+}
+
+template <class Matrix>
+void SparseDirectMUMPS::initialize (const Matrix& matrix)
+{
+                         // Initialize MUMPS instance:
+  initialize_matrix (matrix);
+                         // Start factorization
+  id.job = 4;
+  dmumps_c (&id);
+}
+
 void SparseDirectMUMPS::solve (Vector<double>& vector)
 {
                         // Check that the solver has been initialized
@@ -2045,21 +2089,38 @@ void SparseDirectMUMPS::solve (Vector<double>& vector)
                         // Start solver
   id.job = 6;
   dmumps_c (&id);
-  id.job = -2;
-  dmumps_c (&id);
+  copy_solution (vector);
+}
+
+void SparseDirectMUMPS::vmult (Vector<double>&       dst,
+                               const Vector<double>& src)
+{
+                        // Check that the solver has been initialized
+                        // by the routine above:
+  Assert (initialize_called == true, ExcNotInitialized());
 
-                        // Copy solution into the given vector
+                        // and that the matrix has at least one
+                        // nonzero element:
+  Assert (nz != 0, ExcNotInitialized());
+
+                       // Hand over right-hand side
   if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0)
     {
-      for (unsigned int i=0; i<n; ++i)
-        vector(i) = rhs[i];
+                           // Object denoting a MUMPS data structure:
+      rhs = new double[n];
 
-      delete[] a;
-      delete[] irn;
-      delete[] jcn;
-      delete[] rhs;
+      for (unsigned int i = 0; i < n; ++i)
+        rhs[i] = src (i);
+
+      id.rhs = rhs;
     }
+  
+                       // Start solver
+  id.job = 3;
+  dmumps_c (&id);
+  copy_solution (dst);
 }
+
 #endif // DEAL_II_USE_MUMPS
 
 // explicit instantiations for SparseMatrixMA27
-- 
2.39.5