From: Guido Kanschat Date: Thu, 22 Mar 2001 17:42:23 +0000 (+0000) Subject: new class for Schur complements X-Git-Tag: v8.0.0~19509 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=12b352cb18f12bb597e678129fe3103382c051a3;p=dealii.git new class for Schur complements git-svn-id: https://svn.dealii.org/trunk@4265 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/schur_matrix.h b/deal.II/lac/include/lac/schur_matrix.h new file mode 100644 index 0000000000..eb4434f4d2 --- /dev/null +++ b/deal.II/lac/include/lac/schur_matrix.h @@ -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 +#include +#include +#include +#include + +template 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 SchurMatrix : + public Subscriptor +{ + public: + SchurMatrix(const MA_inverse& Ainv, + const MB& B, + const MDt& Dt, + const MC& C, + VectorMemory >& 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& dst, + const BlockVector& src) const; + + /** + * Multiplication with the Schur + * complement. + */ + void vmult (BlockVector& dst, + const BlockVector& src) const; + +// void Tmult(BlockVector& dst, const BlockVector& src) const; + + /** + * Computation of the residual of + * the Schur complement. + */ + double residual (BlockVector& dst, + const BlockVector& src, + const BlockVector& rhs) const; + + /** + * Compute the eliminated variable + * from the solution of the Schur + * complement problem. + */ + void postprocess (BlockVector& dst, + const BlockVector& src, + const BlockVector& rhs) const; + private: + const SmartPointer Ainv; + const SmartPointer B; + const SmartPointer Dt; + const SmartPointer C; + VectorMemory >& mem; +}; + +template +SchurMatrix +::SchurMatrix(const MA_inverse& Ainv, + const MB& B, + const MDt& Dt, + const MC& C, + VectorMemory >& mem) + : Ainv(&Ainv), B(&B), Dt(&Dt), C(&C), mem(mem) +{} + + +template +void SchurMatrix +::vmult(BlockVector& dst, + const BlockVector& src) const +{ + deallog.push("Schur"); + C->vmult(dst, src); + BlockVector* h1 = mem.alloc(); + h1->reinit(B->n_block_cols(), src.block(0).size()); + Dt->Tvmult(*h1,src); + BlockVector* h2 = mem.alloc(); + h2->reinit(*h1); + Ainv->vmult(*h2, *h1); + mem.free(h1); + B->vmult_add(dst, *h2); + mem.free(h2); + deallog.pop(); +} + + +template +double SchurMatrix +::residual(BlockVector& dst, + const BlockVector& src, + const BlockVector& rhs) const +{ + vmult(dst, src); + dst.scale(-1.); + dst += rhs; + return dst.l2_norm(); +} + + +template +void SchurMatrix +::prepare_rhs(BlockVector& dst, + const BlockVector& 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* 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 +void SchurMatrix +::postprocess(BlockVector& dst, + const BlockVector& src, + const BlockVector& 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* 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