From: kanschat Date: Thu, 28 Jul 2011 17:48:13 +0000 (+0000) Subject: add function for Schur complement; needs testing X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c65c8f26847f2a944a5cef7c8a07595c6da6d6a3;p=dealii-svn.git add function for Schur complement; needs testing git-svn-id: https://svn.dealii.org/trunk@23981 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/lac/full_matrix.h b/deal.II/include/deal.II/lac/full_matrix.h index 7d5031207d..94f3f70e19 100644 --- a/deal.II/include/deal.II/lac/full_matrix.h +++ b/deal.II/include/deal.II/lac/full_matrix.h @@ -1,7 +1,7 @@ //--------------------------------------------------------------------------- // $Id$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -1156,7 +1156,23 @@ class FullMatrix : public Table<2,number> void TmTmult (FullMatrix &C, const FullMatrix &B, const bool adding=false) const; - + + /** + * Add to the current matrix the + * Schur complement B + * A-1 + * D. Optionally, use the + * transposes of the matrices + * B and D. Note + * that the argument for A + * is already the inverse. + */ + void schur_complement(const FullMatrix& Ainverse, + const FullMatrix& B, + const FullMatrix& D, + const bool transpose_B = false, + const bool transpose_D = false); + /** * Matrix-vector-multiplication. * diff --git a/deal.II/include/deal.II/lac/full_matrix.templates.h b/deal.II/include/deal.II/lac/full_matrix.templates.h index 7e8337f1c5..5b756750dd 100644 --- a/deal.II/include/deal.II/lac/full_matrix.templates.h +++ b/deal.II/include/deal.II/lac/full_matrix.templates.h @@ -880,6 +880,67 @@ void FullMatrix::TmTmult (FullMatrix &dst, } +template +void +FullMatrix::schur_complement( + const FullMatrix& Ainverse, + const FullMatrix& B, + const FullMatrix& D, + const bool transpose_B, + const bool transpose_D) +{ + AssertDimension (m(), n()); + AssertDimension (Ainverse.m(), Ainverse.n()); + + const unsigned int N = n(); + const unsigned int M = Ainverse.n(); + + if (transpose_B) + { + AssertDimension(B.m(), M); + AssertDimension(B.n(), N); + } + else + { + AssertDimension(B.n(), M); + AssertDimension(B.m(), N); + } + if (transpose_D) + { + AssertDimension(D.n(), M); + AssertDimension(D.m(), N); + } + else + { + AssertDimension(D.m(), M); + AssertDimension(D.n(), N); + } + + // For all entries of the product + // AD + for (unsigned int i=0; ioperator()(k,j) += ADij*B(i,k); + else + for (unsigned int k=0;koperator()(k,j) += ADij*B(k,i); + } +} + template template