<br>
(Wolfgang Bangerth, 2011/01/18)
-<li> Extended: Several missing instantiations of functions for triangulations and DoF handlers embedded in higher dimensional space have been added.
+<li> Extended: Several missing instantiations of functions for triangulations
+and DoF handlers embedded in higher dimensional space have been added.
<br>
(Wolfgang Bangerth, 2011/01/15)
</ol>
<h3>Specific improvements</h3>
<ol>
+
+<li> New: The function FullMatrix::triple_product() adds triple products
+like Schur complements to existing matrices.
+<br>
+(Guido Kanschat, 2011/08/05)
+
<li> Improved: The PETScWrapper::VectorBase class was rather generous in
calling the PETSc <code>VecAssembleBegin/End</code> functions that incur
communication in the %parallel case and are therefore causes of potential
/**
* Add to the current matrix the
- * Schur complement <b>B
+ * triple product <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);
+ * <b>B</b> and <b>D</b>. 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
+ * <b>B<sup>T</sup>
+ * A<sup>-1</sup> D</b> in mind.
+ * Note that in this case the
+ * argument for <tt>A<tt> must be
+ * the inverse of the matrix <b>A</b>.
+ */
+ void triple_product(const FullMatrix<number>& A,
+ const FullMatrix<number>& B,
+ const FullMatrix<number>& D,
+ const bool transpose_B = false,
+ const bool transpose_D = false,
+ const number scaling = number(1.));
/**
* Matrix-vector-multiplication.
template <typename number>
void
-FullMatrix<number>::schur_complement(
- const FullMatrix<number>& Ainverse,
+FullMatrix<number>::triple_product(
+ const FullMatrix<number>& A,
const FullMatrix<number>& B,
const FullMatrix<number>& 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; i<M;++i)
- for (unsigned int j=0; j<N;++j)
+ for (unsigned int i=0; i<A.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);
+ for (unsigned int k=0;k<A.n();++k)
+ ADij += A(i,k)*D(j,k);
else
- for (unsigned int k=0;k<M;++k)
- ADij += Ainverse(i,k)*D(k,j);
+ for (unsigned int k=0;k<A.n();++k)
+ ADij += A(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);
+ if (transpose_B)
+ for (unsigned int k=0;k<m();++k)
+ this->operator()(k,j) += scaling * ADij * B(i,k);
+ else
+ for (unsigned int k=0;k<m();++k)
+ this->operator()(k,j) += scaling * ADij * B(k,i);
}
}