From b7162574535c9c00272bc324dbb73ed64b2c31f1 Mon Sep 17 00:00:00 2001 From: kanschat Date: Fri, 5 Aug 2011 14:46:28 +0000 Subject: [PATCH] improve and rename triple_product git-svn-id: https://svn.dealii.org/trunk@24018 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 9 +++- deal.II/include/deal.II/lac/full_matrix.h | 32 +++++++---- .../deal.II/lac/full_matrix.templates.h | 53 +++++++++---------- 3 files changed, 54 insertions(+), 40 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 73f578c4f2..f27eaee2dc 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -240,7 +240,8 @@ should be fixed now.
(Wolfgang Bangerth, 2011/01/18) -
  • Extended: Several missing instantiations of functions for triangulations and DoF handlers embedded in higher dimensional space have been added. +
  • Extended: Several missing instantiations of functions for triangulations +and DoF handlers embedded in higher dimensional space have been added.
    (Wolfgang Bangerth, 2011/01/15) @@ -251,6 +252,12 @@ should be fixed now.

    Specific improvements

      + +
    1. New: The function FullMatrix::triple_product() adds triple products +like Schur complements to existing matrices. +
      +(Guido Kanschat, 2011/08/05) +
    2. Improved: The PETScWrapper::VectorBase class was rather generous in calling the PETSc VecAssembleBegin/End functions that incur communication in the %parallel case and are therefore causes of potential diff --git a/deal.II/include/deal.II/lac/full_matrix.h b/deal.II/include/deal.II/lac/full_matrix.h index 1f61a511a8..74477a446b 100644 --- a/deal.II/include/deal.II/lac/full_matrix.h +++ b/deal.II/include/deal.II/lac/full_matrix.h @@ -1159,19 +1159,31 @@ class FullMatrix : public Table<2,number> /** * Add to the current matrix the - * Schur complement B + * triple product 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); + * B and D. The + * scaling factor scales the + * whole product, which is + * helpful when adding a multiple + * of the triple product to the + * matrix. + * + * This product was written with + * the Schur complement + * BT + * A-1 D in mind. + * Note that in this case the + * argument for A must be + * the inverse of the matrix A. + */ + void triple_product(const FullMatrix& A, + const FullMatrix& B, + const FullMatrix& D, + const bool transpose_B = false, + const bool transpose_D = false, + const number scaling = number(1.)); /** * 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 5b756750dd..a912f59811 100644 --- a/deal.II/include/deal.II/lac/full_matrix.templates.h +++ b/deal.II/include/deal.II/lac/full_matrix.templates.h @@ -882,62 +882,57 @@ void FullMatrix::TmTmult (FullMatrix &dst, template void -FullMatrix::schur_complement( - const FullMatrix& Ainverse, +FullMatrix::triple_product( + const FullMatrix& A, const FullMatrix& B, const FullMatrix& D, const bool transpose_B, - const bool transpose_D) + const bool transpose_D, + const number scaling) { - 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); + AssertDimension(B.m(), A.m()); + AssertDimension(B.n(), m()); } else { - AssertDimension(B.n(), M); - AssertDimension(B.m(), N); + AssertDimension(B.n(), A.m()); + AssertDimension(B.m(), m()); } if (transpose_D) { - AssertDimension(D.n(), M); - AssertDimension(D.m(), N); + AssertDimension(D.n(), A.n()); + AssertDimension(D.m(), n()); } else { - AssertDimension(D.m(), M); - AssertDimension(D.n(), N); + AssertDimension(D.m(), A.n()); + 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); + if (transpose_B) + for (unsigned int k=0;koperator()(k,j) += scaling * ADij * B(i,k); + else + for (unsigned int k=0;koperator()(k,j) += scaling * ADij * B(k,i); } } -- 2.39.5