//---------------------------------------------------------------------------
// $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
void TmTmult (FullMatrix<number2> &C,
const FullMatrix<number2> &B,
const bool adding=false) const;
-
+
+ /**
+ * Add to the current matrix the
+ * Schur complement <b>B
+ * A<sup>-1</sup>
+ * D</b>. Optionally, use the
+ * transposes of the matrices
+ * <b>B</b> and <b>D</b>. Note
+ * that the argument for <b>A</b>
+ * is already the inverse.
+ */
+ void schur_complement(const FullMatrix<number>& Ainverse,
+ const FullMatrix<number>& B,
+ const FullMatrix<number>& D,
+ const bool transpose_B = false,
+ const bool transpose_D = false);
+
/**
* Matrix-vector-multiplication.
*
}
+template <typename number>
+void
+FullMatrix<number>::schur_complement(
+ const FullMatrix<number>& Ainverse,
+ const FullMatrix<number>& B,
+ const FullMatrix<number>& 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; i<M;++i)
+ for (unsigned int j=0; j<N;++j)
+ {
+ // Compute the entry
+ number ADij = 0.;
+ if (transpose_D)
+ for (unsigned int k=0;k<M;++k)
+ ADij += Ainverse(i,k)*D(j,k);
+ else
+ for (unsigned int k=0;k<M;++k)
+ ADij += Ainverse(i,k)*D(k,j);
+ // And add it to this after
+ // multiplying with the right
+ // factor from B
+ if (transpose_B)
+ for (unsigned int k=0;k<N;++k)
+ this->operator()(k,j) += ADij*B(i,k);
+ else
+ for (unsigned int k=0;k<N;++k)
+ this->operator()(k,j) += ADij*B(k,i);
+ }
+}
+
template <typename number>
template <typename number2>