]> https://gitweb.dealii.org/ - dealii.git/commitdiff
new class for Schur complements
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Thu, 22 Mar 2001 17:42:23 +0000 (17:42 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Thu, 22 Mar 2001 17:42:23 +0000 (17:42 +0000)
git-svn-id: https://svn.dealii.org/trunk@4265 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/schur_matrix.h [new file with mode: 0644]

diff --git a/deal.II/lac/include/lac/schur_matrix.h b/deal.II/lac/include/lac/schur_matrix.h
new file mode 100644 (file)
index 0000000..eb4434f
--- /dev/null
@@ -0,0 +1,221 @@
+//-------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-------------------------------------------------------------------
+#ifndef __deal2__schur_matrix_h
+#define __deal2__schur_matrix_h
+
+#include <base/subscriptor.h>
+#include <base/smartpointer.h>
+#include <base/logstream.h>
+#include <lac/vector_memory.h>
+#include <vector>
+
+template <typename> class BlockVector;
+
+
+/**
+ * Schur complement of a block matrix.
+ *
+ * Given a non-singular matrix @p{A} (often positive definite) and a
+ * positive semi-definite matrix @p{C} as well as matrices @p{B} and
+ * @{Dt} of full rank, this class implements a new matrix, the Schur
+ * complement a the system of equations of the structure
+ *
+ * \begin{verbatim}
+ * /        \  /   \     /   \
+ * |  A  Dt |  | u |  -  | f |
+ * | -B  C  |  | p |  -  | g |
+ * \        /  \   /     \   /
+ * \end{verbatim}
+ *
+ * Multiplication with the Schur matrix @p{S} is the operation
+ * \begin{verbatim}
+ * S p = C p + B A-inverse Dt-transpose p,
+ * \end{verbatim}
+ * which is an operation within the space for @p{p}.
+ *
+ * The data handed to the Schur matrix are as follows:
+ *
+ * @p{A}: the inverse of @p{A} is stored, instead of @p{A}. This
+ * allows the application to use the most efficient form of inversion,
+ * iterative or direct.
+ *
+ * @p{B}, @p{C}: these matrices are stored "as is".
+ *
+ * @p{Dt}: the computation of the Schur complement involves the
+ * function @p{Tvmult} of the matrix @p{Dt}, not @p{vmult}! This way,
+ * it is sufficient to build only one matrix @p{B} for the symmetric
+ * Schur complement and use it twice.
+ *
+ * All matrices involved are of arbitrary type and vectors are
+ * @ref{BlockVector}s. This way, @p{SchurMatrix} can be coupled with
+ * any matrix classes providing @p{vmult} and @p{Tvmult} and even
+ * nested.
+ *
+ * Since the Schur complement of a matrix corresponds to a Gaussian
+ * block elimination, the right hand side of the condensed system must
+ * be preprocessed. Furthermore, the eliminated variable must be
+ * reconstructed after solving.The solution of the system above by a
+ * @p{SchurMatrix} @p{schur} is coded as follows:
+ *
+ * \begin{verbatim}
+ *   schur.prepare_rhs (g, f);
+ *   solver.solve (schur, p, g, precondition);
+ *   schur.postprocess (u, p);
+ * \end{verbatim}
+ *
+ * @author Guido Kanschat, 2000, 2001
+ */
+template <class MA_inverse, class MB, class MDt, class MC>
+class SchurMatrix :
+  public Subscriptor
+{
+  public:
+  SchurMatrix(const MA_inverse& Ainv,
+             const MB& B,
+             const MDt& Dt,
+             const MC& C,
+             VectorMemory<BlockVector<double> >& mem);
+
+                                  /**
+                                   * Do block elimination of the
+                                   * right hand side. Given right
+                                   * hand sides for both components
+                                   * of the block system, this
+                                   * function provides the right hand
+                                   * side for the Schur complement.
+                                  */
+  void prepare_rhs (BlockVector<double>& dst,
+                   const BlockVector<double>& src) const;
+
+                                  /**
+                                   * Multiplication with the Schur
+                                   * complement.
+                                   */
+  void vmult (BlockVector<double>& dst,
+             const BlockVector<double>& src) const;
+
+//  void Tmult(BlockVector<double>& dst, const BlockVector<double>& src) const;
+
+                                  /**
+                                   * Computation of the residual of
+                                   * the Schur complement.
+                                   */
+  double residual (BlockVector<double>& dst,
+                  const BlockVector<double>& src,
+                  const BlockVector<double>& rhs) const;
+
+                                  /**
+                                   * Compute the eliminated variable
+                                   * from the solution of the Schur
+                                   * complement problem.
+                                   */
+  void postprocess (BlockVector<double>& dst,
+                   const BlockVector<double>& src,
+                   const BlockVector<double>& rhs) const;
+  private:
+  const SmartPointer<const MA_inverse> Ainv;
+  const SmartPointer<const MB> B;
+  const SmartPointer<const MDt> Dt;
+  const SmartPointer<const MC> C;
+  VectorMemory<BlockVector<double> >& mem;
+};
+
+template <class MA_inverse, class MB, class MDt, class MC>
+SchurMatrix<MA_inverse, MB, MDt, MC>
+::SchurMatrix(const MA_inverse& Ainv,
+             const MB& B,
+             const MDt& Dt,
+             const MC& C,
+             VectorMemory<BlockVector<double> >& mem)
+  : Ainv(&Ainv), B(&B), Dt(&Dt), C(&C), mem(mem)
+{}
+
+
+template <class MA_inverse, class MB, class MDt, class MC>
+void SchurMatrix<MA_inverse, MB, MDt, MC>
+::vmult(BlockVector<double>& dst,
+       const BlockVector<double>& src) const
+{
+  deallog.push("Schur");
+  C->vmult(dst, src);
+  BlockVector<double>* h1 = mem.alloc();
+  h1->reinit(B->n_block_cols(), src.block(0).size());
+  Dt->Tvmult(*h1,src);
+  BlockVector<double>* h2 = mem.alloc();
+  h2->reinit(*h1);
+  Ainv->vmult(*h2, *h1);
+  mem.free(h1);
+  B->vmult_add(dst, *h2);
+  mem.free(h2);
+  deallog.pop();
+}
+
+
+template <class MA_inverse, class MB, class MDt, class MC>
+double SchurMatrix<MA_inverse, MB, MDt, MC>
+::residual(BlockVector<double>& dst,
+          const BlockVector<double>& src,
+          const BlockVector<double>& rhs) const
+{
+  vmult(dst, src);
+  dst.scale(-1.);
+  dst += rhs;
+  return dst.l2_norm();
+}
+
+
+template <class MA_inverse, class MB, class MDt, class MC>
+void SchurMatrix<MA_inverse, MB, MDt, MC>
+::prepare_rhs(BlockVector<double>& dst,
+             const BlockVector<double>& src) const
+{
+  Assert (src.n_blocks() == B->n_block_cols(),
+         ExcDimensionMismatch(src.n_blocks(), B->n_block_cols()));
+  Assert (dst.n_blocks() == B->n_block_rows(),
+         ExcDimensionMismatch(dst.n_blocks(), B->n_block_rows()));
+  
+  deallog.push("Schur-prepare");
+  BlockVector<double>* h1 = mem.alloc();
+  h1->reinit(B->n_block_cols(), src.block(0).size());
+  Ainv->vmult(*h1, src);
+  B->vmult_add(dst, *h1);
+  mem.free(h1);
+  deallog.pop();
+}
+
+
+template <class MA_inverse, class MB, class MDt, class MC>
+void SchurMatrix<MA_inverse, MB, MDt, MC>
+::postprocess(BlockVector<double>& dst,
+             const BlockVector<double>& src,
+             const BlockVector<double>& rhs) const
+{
+  Assert (dst.n_blocks() == B->n_block_cols(),
+         ExcDimensionMismatch(dst.n_blocks(), B->n_block_cols()));
+  Assert (rhs.n_blocks() == B->n_block_cols(),
+         ExcDimensionMismatch(rhs.n_blocks(), B->n_block_cols()));
+  Assert (src.n_blocks() == B->n_block_rows(),
+         ExcDimensionMismatch(src.n_blocks(), B->n_block_rows()));
+  
+  deallog.push("Schur-post");
+  BlockVector<double>* h1 = mem.alloc();
+  h1->reinit(B->n_block_cols(), src.block(0).size());
+  Dt->Tvmult(*h1, src);
+  h1->sadd(-1.,rhs);
+  Ainv->vmult(dst,*h1);
+  mem.free(h1);
+  deallog.pop();
+}
+
+
+#endif

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.