From c7bfada10210045a82344797b13081ca78e78e24 Mon Sep 17 00:00:00 2001 From: maier Date: Wed, 3 Jul 2013 09:01:45 +0000 Subject: [PATCH] Implement SparseDirectUMFPACK::Tvmult git-svn-id: https://svn.dealii.org/trunk@29923 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 5 ++ deal.II/include/deal.II/lac/sparse_direct.h | 35 ++++++---- deal.II/source/lac/sparse_direct.cc | 72 ++++++++++++++------- 3 files changed, 76 insertions(+), 36 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index e2b2a59a36..6633f1bf27 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -147,6 +147,11 @@ this function.
    +
  1. New: The function SparseDirectUMFPACK::Tvmult is now implemented. +
    +(Matthias Maier, 2013/07/03) +
    2. +
  2. 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 9146286e36..fe73880c34 100644 --- 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 Ax = b. As this + * this method returns an approximate solution of Ax = b. 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 &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 &dst, const Vector &src) const; + /** + * Same as before, but for block vectors + */ + void Tvmult (BlockVector &dst, + const BlockVector &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 &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 &dst, const Vector &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 &rhs_and_solution) const; + void solve (Vector &rhs_and_solution, bool transpose = false) const; /** * Same as before, but for block vectors. */ - void solve (BlockVector &rhs_and_solution) const; - + void solve (BlockVector &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 void solve (const Matrix &matrix, - Vector &rhs_and_solution); + Vector &rhs_and_solution, + bool transpose = false); /** * Same as before, but for block vectors. */ template void solve (const Matrix &matrix, - BlockVector &rhs_and_solution); + BlockVector &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 3a5c543c99..c8c91b628c 100644 --- 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 &rhs_and_solution) const +SparseDirectUMFPACK::solve (Vector &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 &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 &rhs_and_solution) const void -SparseDirectUMFPACK::solve (BlockVector &rhs_and_solution) const +SparseDirectUMFPACK::solve (BlockVector &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 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 &rhs_and_solution) const template void SparseDirectUMFPACK::solve (const Matrix &matrix, - Vector &rhs_and_solution) + Vector &rhs_and_solution, + bool transpose /*=false*/) { factorize (matrix); - solve (rhs_and_solution); + solve (rhs_and_solution, transpose); } template void -SparseDirectUMFPACK::solve (const Matrix &matrix, - BlockVector &rhs_and_solution) +SparseDirectUMFPACK::solve (const Matrix &matrix, + BlockVector &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 &, - const Vector &) const + Vector &dst, + const Vector &src) const { - Assert(false, ExcNotImplemented()); + dst = src; + this->solve(dst, /*transpose=*/ true); +} + + + +void +SparseDirectUMFPACK::Tvmult ( + BlockVector &dst, + const BlockVector &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 &dst, // explicit instantiations for SparseMatrixUMFPACK -#define InstantiateUMFPACK(MATRIX) \ - template \ - void SparseDirectUMFPACK::factorize (const MATRIX &); \ - template \ - void SparseDirectUMFPACK::solve (const MATRIX &, \ - Vector &); \ - template \ - void SparseDirectUMFPACK::solve (const MATRIX &, \ - BlockVector &); \ - template \ - void SparseDirectUMFPACK::initialize (const MATRIX &, \ +#define InstantiateUMFPACK(MATRIX) \ + template \ + void SparseDirectUMFPACK::factorize (const MATRIX &); \ + template \ + void SparseDirectUMFPACK::solve (const MATRIX &, \ + Vector &, \ + bool); \ + template \ + void SparseDirectUMFPACK::solve (const MATRIX &, \ + BlockVector &, \ + bool); \ + template \ + void SparseDirectUMFPACK::initialize (const MATRIX &, \ const AdditionalData); InstantiateUMFPACK(SparseMatrix) @@ -663,8 +685,8 @@ InstantiateUMFPACK(BlockSparseMatrix) // explicit instantiations for SparseDirectMUMPS #ifdef DEAL_II_WITH_MUMPS -#define InstantiateMUMPS(MATRIX) \ - template \ +#define InstantiateMUMPS(MATRIX) \ + template \ void SparseDirectMUMPS::initialize (const MATRIX &, const Vector &); InstantiateMUMPS(SparseMatrix) -- 2.39.5