Implement SparseDirectUMFPACK::Tvmult

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

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index e2b2a59a3643fb705920accce4eb39765a8b64e3..6633f1bf27ceea27c114b7deb4861674fae33fb9 100644 (file)
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -147,6 +147,11 @@ this function.
 
 <ol>
 
+<li> New: The function SparseDirectUMFPACK::Tvmult is now implemented.
+<br>
+(Matthias Maier, 2013/07/03)
+</li>
+
 <li> New: In addition to the FEValuesExtractors::Scalar,
 FEValuesExtractors::Vector, and FEValuesExtractors::SymmetricTensor classes,
 there are now also fully featured FEValuesExtractors::Tensor extractors

diff --git a/deal.II/include/deal.II/lac/sparse_direct.h b/deal.II/include/deal.II/lac/sparse_direct.h
index 9146286e36974c7ec737b73a9d0db7f60d3b4dc4..fe73880c3492484796a58a1956d3b52dfd59a3af 100644 (file)
--- a/deal.II/include/deal.II/lac/sparse_direct.h
+++ b/deal.II/include/deal.II/lac/sparse_direct.h
@@ -108,7 +108,7 @@ public:
   /**
    * @{
    */
-   
+
   /**
    * This function does nothing. It is only here to provide a interface
    * consistent with other sparse direct solvers.
@@ -153,10 +153,10 @@ public:
   /**
    * @{
    */
-  
+
   /**
    * Preconditioner interface function. Usually, given the source vector,
-   * this method returns an approximated solution of <i>Ax = b</i>. As this
+   * this method returns an approximate solution of <i>Ax = b</i>. As this
    * class provides a wrapper to a direct solver, here it is actually the
    * exact solution (exact within the range of numerical accuracy of
    * course).
@@ -174,11 +174,18 @@ public:
              const BlockVector<double> &src) const;
 
   /**
-   * Not implemented but necessary for compiling certain other classes.
+   * Same as before, but uses the transpose of the matrix, i.e. this
+   * function multiplies with $A^{-T}$.
    */
   void Tvmult (Vector<double> &dst,
               const Vector<double> &src) const;
 
+  /**
+   * Same as before, but for block vectors
+   */
+  void Tvmult (BlockVector<double> &dst,
+              const BlockVector<double> &src) const;
+
   /**
    * Same as vmult(), but adding to the previous solution. Not implemented
    * yet but necessary for compiling certain other classes.
@@ -187,11 +194,12 @@ public:
                  const Vector<double> &src) const;
 
   /**
-   * Not implemented but necessary for compiling certain other classes.
+   * Same as before, but uses the transpose of the matrix, i.e. this
+   * function multiplies with $A^{-T}$.
    */
   void Tvmult_add (Vector<double> &dst,
                   const Vector<double> &src) const;
-  
+
   /**
    * @}
    */
@@ -216,14 +224,17 @@ public:
    * happen. Note that we can't actually call the factorize() function from
    * here if it has not yet been called, since we have no access to the
    * actual matrix.
+   *
+   * If @p transpose is set to true this function solves for the transpose
+   * of the matrix, i.e. $x=A^{-T}b$.
    */
-  void solve (Vector<double> &rhs_and_solution) const;
+  void solve (Vector<double> &rhs_and_solution, bool transpose = false) const;
 
   /**
    * Same as before, but for block vectors.
    */
-  void solve (BlockVector<double> &rhs_and_solution) const;
-  
+  void solve (BlockVector<double> &rhs_and_solution, bool transpose = false) const;
+
   /**
    * Call the two functions factorize() and solve() in that order, i.e. perform
    * the whole solution process for the given right hand side vector.
@@ -232,14 +243,16 @@ public:
    */
   template <class Matrix>
   void solve (const Matrix   &matrix,
-              Vector<double> &rhs_and_solution);
+              Vector<double> &rhs_and_solution,
+              bool            transpose = false);
 
   /**
    * Same as before, but for block vectors.
    */
   template <class Matrix>
   void solve (const Matrix        &matrix,
-              BlockVector<double> &rhs_and_solution);
+              BlockVector<double> &rhs_and_solution,
+              bool                 transpose = false);
 
   /**
    * @}

diff --git a/deal.II/source/lac/sparse_direct.cc b/deal.II/source/lac/sparse_direct.cc
index 3a5c543c9919f9e975afc7168d9cc483b8c90670..c8c91b628c525bd18b736e4f32640752d3fb9dcf 100644 (file)
--- a/deal.II/source/lac/sparse_direct.cc
+++ b/deal.II/source/lac/sparse_direct.cc
@@ -291,7 +291,8 @@ factorize (const Matrix &matrix)
 
 
 void
-SparseDirectUMFPACK::solve (Vector<double> &rhs_and_solution) const
+SparseDirectUMFPACK::solve (Vector<double> &rhs_and_solution,
+                            bool            transpose /*=false*/) const
 {
   // make sure that some kind of factorize() call has happened before
   Assert (Ap.size() != 0, ExcNotInitialized());
@@ -304,8 +305,11 @@ SparseDirectUMFPACK::solve (Vector<double> &rhs_and_solution) const
   // solve the system. note that since UMFPACK wants compressed column
   // storage instead of the compressed row storage format we use in
   // deal.II's SparsityPattern classes, we solve for UMFPACK's A^T instead
+
+  // Conversely, if we solve for the transpose, we have to use UMFPACK_A
+  // instead.
   const int status
-    = umfpack_dl_solve (UMFPACK_At,
+    = umfpack_dl_solve (transpose ? UMFPACK_A : UMFPACK_At,
                         &Ap[0], &Ai[0], &Ax[0],
                         rhs_and_solution.begin(), rhs.begin(),
                         numeric_decomposition,
@@ -315,14 +319,15 @@ SparseDirectUMFPACK::solve (Vector<double> &rhs_and_solution) const
 
 
 void
-SparseDirectUMFPACK::solve (BlockVector<double> &rhs_and_solution) const
+SparseDirectUMFPACK::solve (BlockVector<double> &rhs_and_solution,
+                            bool                 transpose /*=false*/) const
 {
   // the UMFPACK functions want a contiguous array of elements, so
   // there is no way around copying data around. thus, just copy the
   // data into a regular vector and back
   Vector<double> tmp (rhs_and_solution.size());
   tmp = rhs_and_solution;
-  solve (tmp);
+  solve (tmp, transpose);
   rhs_and_solution = tmp;
 }
 
@@ -331,20 +336,22 @@ SparseDirectUMFPACK::solve (BlockVector<double> &rhs_and_solution) const
 template <class Matrix>
 void
 SparseDirectUMFPACK::solve (const Matrix   &matrix,
-                            Vector<double> &rhs_and_solution)
+                            Vector<double> &rhs_and_solution,
+                            bool            transpose /*=false*/)
 {
   factorize (matrix);
-  solve (rhs_and_solution);
+  solve (rhs_and_solution, transpose);
 }
 
 
 template <class Matrix>
 void
-SparseDirectUMFPACK::solve (const Matrix   &matrix,
-                            BlockVector<double> &rhs_and_solution)
+SparseDirectUMFPACK::solve (const Matrix        &matrix,
+                            BlockVector<double> &rhs_and_solution,
+                            bool                 transpose /*=false*/)
 {
   factorize (matrix);
-  solve (rhs_and_solution);
+  solve (rhs_and_solution, transpose);
 }
 
 
@@ -439,10 +446,22 @@ SparseDirectUMFPACK::vmult (
 
 void
 SparseDirectUMFPACK::Tvmult (
-  Vector<double> &,
-  const Vector<double> &) const
+  Vector<double> &dst,
+  const Vector<double> &src) const
 {
-  Assert(false, ExcNotImplemented());
+  dst = src;
+  this->solve(dst, /*transpose=*/ true);
+}
+
+
+
+void
+SparseDirectUMFPACK::Tvmult (
+  BlockVector<double>       &dst,
+  const BlockVector<double> &src) const
+{
+  dst = src;
+  this->solve(dst, /*transpose=*/ true);
 }
 
 
@@ -465,6 +484,7 @@ SparseDirectUMFPACK::Tvmult_add (
 
 
 
+
 #ifdef DEAL_II_WITH_MUMPS
 
 SparseDirectMUMPS::SparseDirectMUMPS ()
@@ -641,17 +661,19 @@ void SparseDirectMUMPS::vmult (Vector<double>       &dst,
 
 
 // explicit instantiations for SparseMatrixUMFPACK
-#define InstantiateUMFPACK(MATRIX) \
-  template    \
-  void SparseDirectUMFPACK::factorize (const MATRIX &);    \
-  template    \
-  void SparseDirectUMFPACK::solve (const MATRIX   &,    \
-                                   Vector<double> &);    \
-  template    \
-  void SparseDirectUMFPACK::solve (const MATRIX   &,    \
-                                   BlockVector<double> &);    \
-  template    \
-  void SparseDirectUMFPACK::initialize (const MATRIX &,    \
+#define InstantiateUMFPACK(MATRIX)                        \
+  template                                                \
+  void SparseDirectUMFPACK::factorize (const MATRIX &);   \
+  template                                                \
+  void SparseDirectUMFPACK::solve (const MATRIX   &,      \
+                                   Vector<double> &,      \
+                                   bool);                 \
+  template                                                \
+  void SparseDirectUMFPACK::solve (const MATRIX   &,      \
+                                   BlockVector<double> &, \
+                                   bool);                 \
+  template                                                \
+  void SparseDirectUMFPACK::initialize (const MATRIX &,   \
                                         const AdditionalData);
 
 InstantiateUMFPACK(SparseMatrix<double>)
@@ -663,8 +685,8 @@ InstantiateUMFPACK(BlockSparseMatrix<float>)
 
 // explicit instantiations for SparseDirectMUMPS
 #ifdef DEAL_II_WITH_MUMPS
-#define InstantiateMUMPS(MATRIX) \
-  template \
+#define InstantiateMUMPS(MATRIX)                          \
+  template                                                \
   void SparseDirectMUMPS::initialize (const MATRIX &, const Vector<double> &);
 
 InstantiateMUMPS(SparseMatrix<double>)