From: heister <heister@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 25 Nov 2013 16:30:50 +0000 (+0000)
Subject: missing changes in r31796, clean up mumps, add asserts
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=30aba11eaa921009b8c98cec9a404af6457e268a;p=dealii-svn.git

missing changes in r31796, clean up mumps, add asserts

git-svn-id: https://svn.dealii.org/trunk@31797 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index ce205f7e4d..89ec0b712c 100644
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -94,6 +94,11 @@ inconvenience this causes.
 
 
 <ol>
+  <li> Fixed: Missing instantiations of SparseDirectMUMPS have been added.
+  <br>
+  (Timo Heister, 2013/11/25)
+  </li>
+
   <li> New: introduced "make test" that runs a minimal set of tests. We
   encourage every user to run this, especially if they run in to problems.
   The tests are automatically picked depending on the configuration and
diff --git a/deal.II/include/deal.II/lac/sparse_direct.h b/deal.II/include/deal.II/lac/sparse_direct.h
index 8ef14a7a38..c9df62991f 100644
--- a/deal.II/include/deal.II/lac/sparse_direct.h
+++ b/deal.II/include/deal.II/lac/sparse_direct.h
@@ -346,7 +346,7 @@ private:
 #endif // DEAL_II_WITH_MUMPS
 
   double   *a;
-  double   *rhs;
+  std::vector<double> rhs;
   int      *irn;
   int      *jcn;
   types::global_dof_index n;
@@ -364,6 +364,11 @@ private:
    */
   void copy_solution (Vector<double> &vector);
 
+  /**
+   *
+   */
+  void copy_rhs_to_mumps(const Vector<double> &rhs);
+
   /**
    * Flags storing whether the function <tt>initialize ()</tt> has already
    * been called.
@@ -409,7 +414,8 @@ public:
 
   /**
    * A function in which the linear system is solved and the solution
-   * vector is copied into the given <tt>vector</tt>.
+   * vector is copied into the given <tt>vector</tt>. The right-hand side
+   * need to be supplied in initialize(matrix, vector);
    */
   void solve (Vector<double> &vector);
 
diff --git a/deal.II/source/lac/sparse_direct.cc b/deal.II/source/lac/sparse_direct.cc
index ed801cea31..4333b64fa1 100644
--- a/deal.II/source/lac/sparse_direct.cc
+++ b/deal.II/source/lac/sparse_direct.cc
@@ -25,6 +25,7 @@
 #include <iostream>
 #include <list>
 #include <typeinfo>
+#include <vector>
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -513,6 +514,9 @@ SparseDirectMUMPS::~SparseDirectMUMPS ()
 template <class Matrix>
 void SparseDirectMUMPS::initialize_matrix (const Matrix &matrix)
 {
+    Assert(matrix.n() == matrix.m(), ExcMessage("Matrix needs to be square."));
+
+
   // Check we haven't been here before:
   Assert (initialize_called == false, ExcInitializeAlreadyCalled());
 
@@ -586,27 +590,35 @@ void SparseDirectMUMPS::initialize (const Matrix &matrix,
   // Hand over matrix and right-hand side
   initialize_matrix (matrix);
 
+  copy_rhs_to_mumps(vector);
+}
+
+void SparseDirectMUMPS::copy_rhs_to_mumps(const Vector<double> &new_rhs)
+{
+  Assert(n == new_rhs.size(), ExcMessage("Matrix size and rhs length must be equal."));
+
   if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0)
     {
-      // Object denoting a MUMPS data structure
-      rhs = new double[n];
-
+      rhs.resize(n);
       for (size_type i = 0; i < n; ++i)
-        rhs[i] = vector (i);
+        rhs[i] = new_rhs (i);
 
-      id.rhs = rhs;
+      id.rhs = &rhs[0];
     }
 }
 
 void SparseDirectMUMPS::copy_solution (Vector<double> &vector)
 {
+  Assert(n == vector.size(), ExcMessage("Matrix size and solution vector length must be equal."));
+  Assert(n == rhs.size(), ExcMessage("Class not initialized with a rhs vector."));
+
   // Copy solution into the given vector
   if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0)
     {
       for (size_type i=0; i<n; ++i)
         vector(i) = rhs[i];
 
-      delete[] rhs;
+      rhs.resize(0); // remove rhs again
     }
 }
 
@@ -622,6 +634,10 @@ void SparseDirectMUMPS::initialize (const Matrix &matrix)
 
 void SparseDirectMUMPS::solve (Vector<double> &vector)
 {
+  //TODO: this could be implemented similar to SparseDirectUMFPACK where
+  // the given vector will be used as the RHS. Sadly, there is no easy
+  // way to do this without breaking the interface.
+
   // Check that the solver has been initialized by the routine above:
   Assert (initialize_called == true, ExcNotInitialized());
 
@@ -629,7 +645,7 @@ void SparseDirectMUMPS::solve (Vector<double> &vector)
   Assert (nz != 0, ExcNotInitialized());
 
   // Start solver
-  id.job = 6;
+  id.job = 6; // 6 = analysis, factorization, and solve
   dmumps_c (&id);
   copy_solution (vector);
 }
@@ -643,17 +659,11 @@ void SparseDirectMUMPS::vmult (Vector<double>       &dst,
   // 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)
-    {
-      // Object denoting a MUMPS data structure:
-      rhs = new double[n];
-
-      for (size_type i = 0; i < n; ++i)
-        rhs[i] = src (i);
+  Assert(n == dst.size(), ExcMessage("Destination vector has the wrong size."));
+  Assert(n == src.size(), ExcMessage("Source vector has the wrong size."));
 
-      id.rhs = rhs;
-    }
+  // Hand over right-hand side
+  copy_rhs_to_mumps(src);
 
   // Start solver
   id.job = 3;
@@ -691,7 +701,9 @@ InstantiateUMFPACK(BlockSparseMatrix<float>)
 #ifdef DEAL_II_WITH_MUMPS
 #define InstantiateMUMPS(MATRIX)                          \
   template                                                \
-  void SparseDirectMUMPS::initialize (const MATRIX &, const Vector<double> &);
+  void SparseDirectMUMPS::initialize (const MATRIX &, const Vector<double> &);\
+  template                                                \
+  void SparseDirectMUMPS::initialize (const MATRIX &);
 
 InstantiateMUMPS(SparseMatrix<double>)
 InstantiateMUMPS(SparseMatrix<float>)