]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implementation of Schur complement linear operator and auxillary functions 1956/head
authorJean-Paul Pelteret <jppelteret@gmail.com>
Mon, 7 Dec 2015 08:01:21 +0000 (09:01 +0100)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Thu, 7 Jan 2016 21:44:28 +0000 (22:44 +0100)
A new linear operator representing the Schur complement,
namely schur_complement(), has been implemented. Some auxiliary functions
that are often used in conjunction with the Schur complement
(condense_schur_rhs() and postprocess_schur_solution()) are also provided
as a PackagedOperation. Addresses #1439 .

doc/news/changes.h
include/deal.II/lac/schur_complement.h [new file with mode: 0644]
tests/lac/linear_operator_01b.cc [new file with mode: 0644]
tests/lac/linear_operator_01b.with_cxx11=on.output [new file with mode: 0644]
tests/lac/schur_complement_01.cc [new file with mode: 0644]
tests/lac/schur_complement_01.with_cxx11=on.output [new file with mode: 0644]
tests/lac/schur_complement_02.cc [new file with mode: 0644]
tests/lac/schur_complement_02.with_cxx11=on.output [new file with mode: 0644]
tests/lac/schur_complement_03.cc [new file with mode: 0644]
tests/lac/schur_complement_03.with_cxx11=on.output [new file with mode: 0644]

index c903556dc35d1ea87f083a3a5c60c8d3b779725b..22407a5a050689fecd063f8154acee55698a53a9 100644 (file)
@@ -215,6 +215,20 @@ inconvenience this causes.
   (Luca Heltai, 2015/12/13)
   </li>
 
+  <li> New: A new linear operator representing the Schur complement,
+  namely schur_complement(), has been implemented. Some auxiliary functions
+  that are often used in conjunction with the Schur complement
+  (condense_schur_rhs() and postprocess_schur_solution()) are also provided
+  as a PackagedOperation.
+  An example of this functionality can be found in
+  <code>tests/lac/schur_complement_01.cc</code>.
+  The solution of a multi-component problem (namely step-22) using the
+  schur_complement can be found in
+  <code>tests/lac/schur_complement_03.cc</code> .
+  <br>
+  (Jean-Paul Pelteret, Matthias Maier, Martin Kronbichler, 2015/12/07)
+  </li>
+
   <li> New: There is now a function Utilities::to_string that works like
   int_to_string, but is more safe for long integers, negative integers, and
   also handles floating point numbers. The implementation of int_to_string
diff --git a/include/deal.II/lac/schur_complement.h b/include/deal.II/lac/schur_complement.h
new file mode 100644 (file)
index 0000000..45b2c17
--- /dev/null
@@ -0,0 +1,515 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__schur_complement_h
+#define dealii__schur_complement_h
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/lac/vector_memory.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/packaged_operation.h>
+
+#ifdef DEAL_II_WITH_CXX11
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * @name Creation of a LinearOperator related to the Schur Complement
+ */
+//@{
+
+/**
+ * @relates LinearOperator
+ *
+ * Returns a LinearOperator that performs the operations
+ * associated with the Schur complement. There are two additional
+ * helper functions, condense_schur_rhs() and postprocess_schur_solution(), that are likely
+ * necessary to be used in order to perform any useful tasks in linear
+ * algebra with this operator.
+ *
+ * We construct the definition of the Schur complement in the following way:
+ *
+ * Consider a general system of linear equations that can be
+ * decomposed into two major sets of equations:
+ * @f{eqnarray*}
+  \mathbf{K}\mathbf{d} = \mathbf{f}
+  \quad \Rightarrow\quad
+  \left(\begin{array}{cc}
+    A & B \\ C & D
+  \end{array}\right)
+  \left(\begin{array}{cc}
+    x \\ y
+  \end{array}\right)
+  =
+  \left(\begin{array}{cc}
+    f \\ g
+  \end{array}\right),
+ * @f}
+ * where $ A,B,C,D $  represent general subblocks of the matrix
+ * $ \mathbf{K} $ and, similarly, general subvectors of
+ * $ \mathbf{d},\mathbf{f} $ are given by $ x,y,f,g $ .
+ *
+ * This is equivalent to the following two statements:
+ * @f{eqnarray*}
+     (1) \quad Ax + By &=& f \\
+     (2) \quad Cx + Dy &=& g \quad .
+ * @f}
+ *
+ * Assuming that $ A,D $ are both square and invertible, we could
+ * then perform one of two possible substitutions,
+ * @f{eqnarray*}
+     (3) \quad x &=& A^{-1}(f - By) \quad \textnormal{from} \quad (1) \\
+     (4) \quad y &=& D^{-1}(g - Cx) \quad \textnormal{from} \quad (2) ,
+ * @f}
+ * which amount to performing block Gaussian elimination on
+ * this system of equations.
+ *
+ * For the purpose of the current implementation, we choose to
+ * substitute (3) into (2)
+ * @f{eqnarray*}
+     C \: A^{-1}(f - By) + Dy &=& g \\
+     -C \: A^{-1} \: By + Dy &=& g - C \: A^{-1} \: f \quad .
+   @f}
+ * This leads to the result
+ * @f[
+     (5) \quad (D - C\: A^{-1} \:B)y  = g - C \: A^{-1} f
+         \quad \Rightarrow \quad Sy = g'
+   @f]
+ * with $ S = (D - C\: A^{-1} \:B) $ being the Schur complement
+ * and the modified right-hand side vector $ g' = g - C \: A^{-1} f $ arising from
+ * the condensation step.
+ * Note that for this choice of $ S $, submatrix $ D $
+ * need not be invertible and may thus be the null matrix.
+ * Ideally $ A $ should be well-conditioned.
+ *
+ * So for any arbitrary vector $ a $, the Schur complement
+ * performs the following operation:
+ * @f[
+     (6) \quad Sa = (D - C \: A^{-1} \: B)a
+   @f]
+ *
+ * A typical set of steps needed the solve a linear system (1),(2)
+ * would be:
+ * 1. Define the inverse matrix @p A_inv (using inverse_operator()).
+ * 2. Define the Schur complement $ S $ (using schur_complement()).
+ * 3. Define iterative inverse matrix $ S^{-1} $ such that (6)
+ *    holds.
+ *    It is necessary to use a solver with a preconditioner
+ *    to compute the approximate inverse operation of $ S $ since
+ *    we never compute $ S $ directly, but rather the result of
+ *    its operation.
+ *    To achieve this, one may again use the inverse_operator() in
+ *    conjunction with the Schur complement that we've just
+ *    constructed.
+ *    Observe that the both $ S $ and its preconditioner operate
+ *    over the same space as $ D $.
+ * 4. Perform pre-processing step on the RHS of (5) using
+ *    condense_schur_rhs():
+ *    @f[
+        g' = g - C \: A^{-1} \: f
+      @f]
+ * 5. Solve for $ y $ in (5):
+ *    @f[
+        y =  S^{-1} g'
+      @f]
+ * 6. Perform the post-processing step from (3) using
+ *    postprocess_schur_solution():
+ *    @f[
+        x =  A^{-1} (f - By)
+      @f]
+ *
+ * An illustration of typical usage of this operator for a fully coupled
+ * system is given below.
+ * @code
+    #include<deal.II/lac/schur_complement.h>
+
+    // Given BlockMatrix K and BlockVectors d,F
+
+    // Decomposition of tangent matrix
+    const auto A = linear_operator(K.block(0,0));
+    const auto B = linear_operator(K.block(0,1));
+    const auto C = linear_operator(K.block(1,0));
+    const auto D = linear_operator(K.block(1,1));
+
+    // Decomposition of solution vector
+    auto x = d.block(0);
+    auto y = d.block(1);
+
+    // Decomposition of RHS vector
+    auto f = F.block(0);
+    auto g = F.block(1);
+
+    // Construction of inverse of Schur complement
+    const auto prec_A = PreconditionSelector<...>(A);
+    const auto A_inv = inverse_operator<...>(A,prec_A);
+    const auto S = schur_complement(A_inv,B,C,D);
+    const auto S_prec = PreconditionSelector<...>(D); // D and S operate on same space
+    const auto S_inv = inverse_operator<...>(S,...,prec_S);
+
+    // Solve reduced block system
+    auto rhs = condense_schur_rhs (A_inv,C,f,g); // PackagedOperation that represents the condensed form of g
+    y = S_inv * rhs; // Solve for y
+    x = postprocess_schur_solution (A_inv,B,y,f); // Compute x using resolved solution y
+   @endcode
+ *
+ * In the above example, the preconditioner for $ S $ was defined as the
+ * preconditioner for $ D $, which is valid since they operate on the same
+ * space.
+ * However, if $ D $ and $ S $ are too dissimilar, then this may lead to
+ * a large number of solver iterations as $ \text{prec}(D) $ is not a good
+ * approximation for $ S^{-1} $.
+ *
+ * A better preconditioner in such a case would be one that provides a more
+ * representative approximation for $ S^{-1} $.
+ * One approach is shown in step-22, where $ D $ is the null matrix and the
+ * preconditioner for $ S^{-1} $ is derived from the mass matrix over this
+ * space.
+ *
+ * From another viewpoint, a similar result can be achieved by first
+ * constructing an object that represents an approximation for $ S $ wherein
+ * expensive operation, namely $ A^{-1} $, is approximated.
+ * Thereafter we construct the approximate inverse operator $ \tilde{S}^{-1} $
+ * which is then used as the preconditioner for computing $ S^{-1} $.
+ * @code
+    // Construction of approximate inverse of Schur complement
+    const auto A_inv_approx = linear_operator(preconditioner_A);
+    const auto S_approx = schur_complement(A_inv_approx,B,C,D);
+    const auto S_approx_prec = PreconditionSelector<...>(D); // D and S_approx operate on same space
+    const auto S_inv_approx = inverse_operator(S_approx,...,S_approx_prec); // Inner solver: Typically limited to few iterations using IterationNumberControl
+
+    // Construction of exact inverse of Schur complement
+    const auto S = schur_complement(A_inv,B,C,D);
+    const auto S_inv = inverse_operator(S,...,S_inv_approx); // Outer solver
+
+    // Solve reduced block system
+    auto rhs = condense_schur_rhs (A_inv,C,f,g);
+    y = S_inv * rhs; // Solve for y
+    x = postprocess_schur_solution (A_inv,B,y,f);
+   @endcode
+ * Note that due to the construction of @c S_inv_approx and subsequently
+ * @c S_inv, there are a pair of nested iterative solvers which could
+ * collectively consume a lot of resources.
+ * Therefore care should be taken in the choices leading to the construction
+ * of the iterative inverse_operators.
+ * One might consider the use of a IterationNumberControl (or a similar
+ * mechanism) to limit the number of inner solver iterations.
+ * This controls the accuracy of the approximate inverse operation
+ * $ \tilde{S}^{-1} $ which acts only as the preconditioner for
+ * $ S^{-1} $.
+ * Furthermore, the preconditioner to $ \tilde{S}^{-1} $, which in this example is
+ * $ \text{prec}(D) $, should ideally be computationally inexpensive.
+ *
+ * However, if an iterative solver based on IterationNumberControl is used as a
+ * preconditioner then the preconditioning operation is not a linear operation.
+ * Here a flexible solver like SolverFGMRES (flexible GMRES) is best employed as an
+ * outer solver in order to deal with the variable behaviour of the preconditioner.
+ * Otherwise the iterative solver can stagnate somewhere near the tolerance of the
+ * preconditioner or generally behave erratically.
+ * Alternatively, using a ReductionControl would ensure that the preconditioner
+ * always solves to the same tolerance, thereby rendering its behaviour constant.
+ *
+ * Further examples of this functionality can be found in
+ * the test-suite, such as
+ * <code>tests/lac/schur_complement_01.cc</code> .
+ * The solution of a multi-component problem (namely step-22) using the
+ * schur_complement can be found in
+ * <code>tests/lac/schur_complement_03.cc</code> .
+ *
+ * @see
+ * @ref GlossBlockLA "Block (linear algebra)"
+ * @author Jean-Paul Pelteret, Matthias Maier, Martin Kronbichler, 2015
+ *
+ * @ingroup LAOperators
+ */
+template <typename Range_1, typename Domain_1,
+          typename Range_2, typename Domain_2>
+LinearOperator<Range_2, Domain_2>
+schur_complement(const LinearOperator<Domain_1, Range_1> &A_inv,
+                 const LinearOperator<Range_1, Domain_2> &B,
+                 const LinearOperator<Range_2, Domain_1> &C,
+                 const LinearOperator<Range_2, Domain_2> &D)
+{
+  LinearOperator<Range_2, Domain_2> return_op;
+
+  return_op.reinit_range_vector = D.reinit_range_vector;
+  return_op.reinit_domain_vector = D.reinit_domain_vector;
+
+  // ensure to have valid computation objects by catching
+  // A_inv,B,C,D by value
+
+  return_op.vmult_add = [A_inv,B,C,D](Range_2 &dst_g, const Domain_2 &src_y)
+  {
+    static GrowingVectorMemory<Range_1>  vector_memory_f;
+    static GrowingVectorMemory<Range_2>  vector_memory_g;
+    static GrowingVectorMemory<Domain_1> vector_memory_x;
+
+    Range_1  &tmp_f = *(vector_memory_f.alloc());
+    Range_2  &tmp_g = *(vector_memory_g.alloc());
+    Domain_1 &tmp_x = *(vector_memory_x.alloc());
+
+    // Reinitialise in context of how they'll be used
+    B.reinit_range_vector(tmp_f, /*bool omit_zeroing_entries =*/ true);
+    A_inv.reinit_range_vector(tmp_x, /*bool omit_zeroing_entries =*/ true);
+    C.reinit_range_vector(tmp_g, /*bool omit_zeroing_entries =*/ true);
+
+    // Need to form dst_g such that dst_g = S*src_y = (D - C*A_inv*B) src_y
+    if (D.is_null_operator == false)
+      D.vmult_add (dst_g, src_y); // dst_g += D*src_y (length y)
+
+    B.vmult (tmp_f, src_y); // tmp_f = B*src_y (length x)
+    try
+      {
+        A_inv.vmult (tmp_x, tmp_f); // tmp_x = A_inv*B*src_y (length x)
+      }
+    catch (...)
+      {
+        AssertThrow(false,
+                    ExcMessage("No convergence in A_inv vmult operation"));
+      }
+    C.vmult (tmp_g, tmp_x); // tmp_g = C*A_inv*B*src_y (length y)
+    dst_g -= tmp_g; // dst_g += D*src_y - C*A_inv*B*src_y
+
+    vector_memory_x.free(&tmp_x);
+    vector_memory_g.free(&tmp_g);
+    vector_memory_f.free(&tmp_f);
+  };
+
+  const auto vmult_add = return_op.vmult_add;
+  return_op.vmult = [vmult_add](Range_2 &dst_g, const Domain_2 &src_y)
+  {
+    dst_g = 0.;
+    vmult_add(dst_g, src_y);
+  };
+
+  return_op.Tvmult_add = [A_inv,B,C,D](Domain_2 &dst_g, const Range_2 &src_y)
+  {
+    static GrowingVectorMemory<Domain_1> vector_memory_f;
+    static GrowingVectorMemory<Domain_2> vector_memory_g;
+    static GrowingVectorMemory<Range_1>  vector_memory_x;
+
+    Domain_1 &tmp_f = *(vector_memory_f.alloc());
+    Domain_2 &tmp_g = *(vector_memory_g.alloc());
+    Range_1  &tmp_x = *(vector_memory_x.alloc());
+
+    // Reinitialise in context of how they'll be used
+    C.reinit_domain_vector(tmp_f, /*bool omit_zeroing_entries =*/ true);
+    A_inv.reinit_domain_vector(tmp_x, /*bool omit_zeroing_entries =*/ true);
+    B.reinit_domain_vector(tmp_g, /*bool omit_zeroing_entries =*/ true);
+
+    // Need to form y such that dst such that dst_g = S*src_y = (D^T - B^T*A_inv^T*C^T) src_y
+    if (D.is_null_operator == false)
+      D.Tvmult_add (dst_g, src_y); // dst_g += D^T*src_y (length y)
+
+    C.Tvmult (tmp_f, src_y); // tmp_f = C^T*src_y (length x)
+    try
+      {
+        A_inv.Tvmult (tmp_x, tmp_f); // tmp_x = A_inv^T*C^T*src_y (length x)
+      }
+    catch (...)
+      {
+        AssertThrow(false,
+                    ExcMessage("No convergence in A_inv Tvmult operation"));
+      }
+    B.Tvmult (tmp_g, tmp_x); // tmp_g = B^T*A_inv^T*C^T*src_y (length y)
+    dst_g -= tmp_g; // dst_g += D^T*src_y - B^T*A_inv^T*C^T*src_y
+
+    vector_memory_x.free(&tmp_x);
+    vector_memory_g.free(&tmp_g);
+    vector_memory_f.free(&tmp_f);
+  };
+
+  const auto Tvmult_add = return_op.Tvmult_add;
+  return_op.Tvmult = [Tvmult_add](Domain_2 &dst_g, const Range_2 &src_y)
+  {
+    dst_g = 0.;
+    Tvmult_add(dst_g, src_y);
+  };
+
+  return return_op;
+}
+
+//@}
+
+
+/**
+ * @name Creation of PackagedOperation objects related to the Schur Complement
+ */
+//@{
+
+/**
+ * @relates PackagedOperation
+ *
+ * For the system of equations
+ * @f{eqnarray*}
+     Ax + By &=& f \\
+     Cx + Dy &=& g \quad ,
+ * @f}
+ * this operation performs the pre-processing (condensation)
+ * step on the RHS subvector @p g so that the Schur complement
+ * can be used to solve this system of equations.
+ * More specifically, it produces an object that represents the
+ * condensed form of the subvector @p g, namely
+ * @f[
+     g' = g - C \: A^{-1} \: f
+   @f]
+ *
+ * @see
+ * @ref GlossBlockLA "Block (linear algebra)"
+ * @author Jean-Paul Pelteret, Matthias Maier, 2015
+ *
+ * @ingroup LAOperators
+ */
+template <typename Range_1, typename Domain_1,
+          typename Range_2>
+PackagedOperation<Range_2>
+condense_schur_rhs (const LinearOperator<Range_1, Domain_1> &A_inv,
+                    const LinearOperator<Range_2, Domain_1> &C,
+                    const Range_1                           &f,
+                    const Range_2                           &g)
+{
+  PackagedOperation<Range_2> return_comp;
+
+  return_comp.reinit_vector = C.reinit_range_vector;
+
+  // ensure to have valid computation objects by catching
+  // A_inv,C,f,g by value
+
+  return_comp.apply_add = [A_inv,C,f,g](Range_2 &g_star)
+  {
+
+    static GrowingVectorMemory<Range_1> vector_memory_f;
+    static GrowingVectorMemory<Range_2> vector_memory_g;
+
+    Range_1 &tmp_f1 = *(vector_memory_f.alloc());
+    Range_2 &tmp_g1 = *(vector_memory_g.alloc());
+    Range_2 &tmp_g2 = *(vector_memory_g.alloc());
+
+    // Reinitialise in context of how they'll be used
+    A_inv.reinit_range_vector(tmp_f1, /*bool omit_zeroing_entries =*/ true);
+    C.reinit_range_vector(tmp_g1, /*bool omit_zeroing_entries =*/ true);
+
+    // Condensation on RHS of one field
+    // Need to form g* such that g* = g - C*A_inv*f
+    try
+      {
+        A_inv.vmult(tmp_f1, f); // tmp_f1 = A_inv * f
+      }
+    catch (...)
+      {
+        AssertThrow(false,
+                    ExcMessage("No convergence in A_inv vmult operation"));
+      }
+    C.vmult(tmp_g1, tmp_f1); // tmp2 = C * A_inv * f
+
+    g_star += g;
+    g_star -= tmp_g1; // tmp_g2 = g - C * A_inv * f
+
+    vector_memory_g.free(&tmp_g2);
+    vector_memory_g.free(&tmp_g1);
+    vector_memory_f.free(&tmp_f1);
+  };
+
+  const auto apply_add = return_comp.apply_add;
+  return_comp.apply = [apply_add](Range_2 &g_star)
+  {
+    g_star = 0.;
+    apply_add(g_star);
+  };
+
+  return return_comp;
+}
+
+/**
+ * @relates PackagedOperation
+ *
+ * For the system of equations
+ * @f{eqnarray*}
+     Ax + By &=& f \\
+     Cx + Dy &=& g \quad ,
+ * @f}
+ * this operation performs the post-processing step of the
+ * Schur complement to solve for the second subvector @p x once
+ * subvector @p y is known, with the result that
+ * @f[
+     x =  A^{-1}(f - By)
+   @f]
+ *
+ * @see
+ * @ref GlossBlockLA "Block (linear algebra)"
+ * @author Jean-Paul Pelteret, Matthias Maier, 2015
+ *
+ * @ingroup LAOperators
+ */
+template <typename Range_1, typename Domain_1,
+          typename Domain_2>
+PackagedOperation<Domain_1>
+postprocess_schur_solution (const LinearOperator<Range_1, Domain_1> &A_inv,
+                            const LinearOperator<Range_1, Domain_2> &B,
+                            const Domain_2                          &y,
+                            const Range_1                           &f)
+{
+  PackagedOperation<Domain_1> return_comp;
+
+  return_comp.reinit_vector = A_inv.reinit_domain_vector;
+
+  // ensure to have valid computation objects by catching
+  // A_inv,B,y,f by value
+
+  return_comp.apply_add = [A_inv,B,y,f](Domain_1 &x)
+  {
+    static GrowingVectorMemory<Range_1> vector_memory_f;
+
+    Range_1 &tmp_f1 = *(vector_memory_f.alloc());
+    Range_1 &tmp_f2 = *(vector_memory_f.alloc());
+
+    // Reinitialise in context of how they'll be used
+    B.reinit_range_vector(tmp_f1, /*bool omit_zeroing_entries =*/ true);
+
+    // Solve for second field
+    // Need to form x such that x = A_inv*(f - B*y)
+    B.vmult(tmp_f1, y); // tmp_f1 = B*y
+    tmp_f2 = f;
+    tmp_f2 -= tmp_f1; // tmp_f2 = f - B*y
+    try
+      {
+        A_inv.vmult_add(x, tmp_f2); // x = A_inv*(f-B*y)
+      }
+    catch (...)
+      {
+        AssertThrow(false,
+                    ExcMessage("No convergence in A_inv vmult operation"));
+      }
+
+    vector_memory_f.free(&tmp_f2);
+    vector_memory_f.free(&tmp_f1);
+  };
+
+  const auto apply_add = return_comp.apply_add;
+  return_comp.apply = [apply_add](Domain_1 &x)
+  {
+    x = 0.;
+    apply_add(x);
+  };
+
+  return return_comp;
+}
+
+//@}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_CXX11
+#endif
diff --git a/tests/lac/linear_operator_01b.cc b/tests/lac/linear_operator_01b.cc
new file mode 100644 (file)
index 0000000..b21bf86
--- /dev/null
@@ -0,0 +1,168 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Test for composite operations
+
+#include "../tests.h"
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/packaged_operation.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+
+#define PRINTME(name, var) \
+  deallog \
+      << "Solution vector: " << name << ": " \
+      << var;
+
+using namespace dealii;
+
+int main()
+{
+  initlog();
+  deallog << std::setprecision(10);
+
+  // SparseMatrix:
+  {
+    const unsigned int rc=2;
+    SparsityPattern sparsity_pattern (rc, rc, 0);
+    sparsity_pattern.compress();
+
+    SparseMatrix<double> A (sparsity_pattern);
+    Vector<double> b (rc);
+    for (unsigned int i=0; i < rc; ++i)
+      {
+        A.diag_element(i) = 5.0;
+        b(i) = 1.0;
+      }
+
+    const auto lo_A = linear_operator(A);
+    const auto lo_A_t = transpose_operator(lo_A);
+
+    // build transpose of inverse
+    SolverControl solver_control_A_1 (100, 1.0e-10);
+    SolverCG< Vector<double> > solver_A_1(solver_control_A_1);
+    PreconditionJacobi< SparseMatrix<double> > preconditioner_A_1;
+    preconditioner_A_1.initialize(A);
+    const auto lo_A_inv = inverse_operator(lo_A,
+                                           solver_A_1,
+                                           preconditioner_A_1);
+    const auto lo_A_inv_t = transpose_operator(lo_A_inv);
+
+    // build inverse of transpose
+    SolverControl solver_control_A_2 (100, 1.0e-10);
+    SolverCG< Vector<double> > solver_A_2(solver_control_A_2);
+    PreconditionJacobi< SparseMatrix<double> > preconditioner_A_2;
+    preconditioner_A_2.initialize(A);
+    const auto lo_A_t_inv = inverse_operator(lo_A_t,
+                                             solver_A_2,
+                                             preconditioner_A_2);
+
+
+
+    deallog << "Normal and inverse multiplication operations" << std::endl;
+
+    const Vector<double> x1 = lo_A*b;
+    const Vector<double> x2 = lo_A_t *b;
+    const Vector<double> x3 = lo_A_inv*b;
+    const Vector<double> x4a = lo_A_inv_t *b;
+    const Vector<double> x4b = lo_A_t_inv*b;
+
+//   PRINTME("x1", x1);
+//   PRINTME("x2", x2);
+//   PRINTME("x3", x3);
+//   PRINTME("x4a", x4a);
+//   PRINTME("x4b", x4b);
+
+    deallog << "x4a==x4b : " << (x4a==x4b) << std::endl;
+
+    // Schur-type composition
+    SparseMatrix<double> B (sparsity_pattern);
+    SparseMatrix<double> C (sparsity_pattern);
+    SparseMatrix<double> D (sparsity_pattern);
+
+    for (unsigned int i=0; i < rc; ++i)
+      {
+        B.diag_element(i) = 4.0;
+        C.diag_element(i) = 4.0; // K = [A,B ; C,D] is symmetric
+        D.diag_element(i) = 3.0;
+      }
+
+    const auto lo_B = linear_operator(B);
+    const auto lo_C = linear_operator(C);
+    const auto lo_D = linear_operator(D);
+    const auto lo_B_t = transpose_operator(lo_B);
+    const auto lo_C_t = transpose_operator(lo_C);
+    const auto lo_D_t = transpose_operator(lo_D);
+
+    deallog << "Single packaged operation" << std::endl;
+    {
+      const auto S = lo_D - lo_C*lo_A_inv*lo_B;
+      const auto S_t_1 = lo_D_t - lo_B_t *lo_A_inv_t *lo_C_t; // using transpose of inverse
+      const auto S_t_2 = lo_D_t - lo_B_t *lo_A_t_inv*lo_C_t; // using inverse of transpose
+
+      const Vector<double> x5 = S*b;
+      const Vector<double> x6a = S_t_1*b;
+      const Vector<double> x6b = S_t_2*b;
+
+      deallog << "x5==x6a : " << (x5==x6a) << std::endl; // using transpose of inverse
+      deallog << "x5==x6b : " << (x5==x6b) << std::endl; // using inverse of transpose
+//      PRINTME("x5", x5);
+//      PRINTME("x6a", x6a);
+//      PRINTME("x6b", x6b);
+    }
+
+    deallog << "Manual operations" << std::endl;
+    {
+      const Vector<double> x5a = lo_B*b;
+      const Vector<double> x5b = lo_A_inv*x5a;
+      const Vector<double> x5c = lo_C*x5b;
+      const Vector<double> x5d = lo_D*b - x5c;
+      const Vector<double> x6a = lo_C_t *b;
+      const Vector<double> x6b_1 = lo_A_inv_t *x6a; // using transpose of inverse
+      const Vector<double> x6c_1 = lo_B_t *x6b_1;
+      const Vector<double> x6d_1 = lo_D_t *b - x6c_1;
+      const Vector<double> x6b_2 = lo_A_t_inv*x6a; // using inverse of transpose
+      const Vector<double> x6c_2 = lo_B_t *x6b_2;
+      const Vector<double> x6d_2 = lo_D_t *b - x6c_2;
+
+      deallog << "x5a==x6a : " << (x5a==x6a) << std::endl;
+      deallog << "x5b==x6b_1 : " << (x5b==x6b_1) << std::endl; // using transpose of inverse
+      deallog << "x5c==x6c_1 : " << (x5c==x6c_1) << std::endl;
+      deallog << "x5d==x6d_1 : " << (x5d==x6d_1) << std::endl;
+      deallog << "x5b==x6b_2 : " << (x5b==x6b_2) << std::endl; // using inverse of transpose
+      deallog << "x5c==x6c_2 : " << (x5c==x6c_2) << std::endl;
+      deallog << "x5d==x6d_2 : " << (x5d==x6d_2) << std::endl;
+//     PRINTME("x5a", x5a);
+//     PRINTME("x6a", x6a);
+//     PRINTME("x5b", x5b);
+//     PRINTME("x6b_1", x6b_1);
+//     PRINTME("x6b_2", x6b_2);
+//     PRINTME("x5c", x5c);
+//     PRINTME("x6c_1", x6c_1);
+//     PRINTME("x6c_2", x6c_2);
+//     PRINTME("x5d", x5d);
+//     PRINTME("x6d_1", x6d_1);
+//     PRINTME("x6d_2", x6d_2);
+    }
+
+    deallog << "OK" << std::endl;
+  }
+}
+
diff --git a/tests/lac/linear_operator_01b.with_cxx11=on.output b/tests/lac/linear_operator_01b.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..178aac0
--- /dev/null
@@ -0,0 +1,33 @@
+
+DEAL::Normal and inverse multiplication operations
+DEAL:cg::Starting value 1.414213562
+DEAL:cg::Convergence step 1 value 0.000000000
+DEAL:cg::Starting value 1.414213562
+DEAL:cg::Convergence step 1 value 0.000000000
+DEAL:cg::Starting value 1.414213562
+DEAL:cg::Convergence step 1 value 0.000000000
+DEAL::x4a==x4b : 1
+DEAL::Single packaged operation
+DEAL:cg::Starting value 5.656854249
+DEAL:cg::Convergence step 1 value 0.000000000
+DEAL:cg::Starting value 5.656854249
+DEAL:cg::Convergence step 1 value 0.000000000
+DEAL:cg::Starting value 5.656854249
+DEAL:cg::Convergence step 1 value 0.000000000
+DEAL::x5==x6a : 1
+DEAL::x5==x6b : 1
+DEAL::Manual operations
+DEAL:cg::Starting value 5.656854249
+DEAL:cg::Convergence step 1 value 0.000000000
+DEAL:cg::Starting value 5.656854249
+DEAL:cg::Convergence step 1 value 0.000000000
+DEAL:cg::Starting value 5.656854249
+DEAL:cg::Convergence step 1 value 0.000000000
+DEAL::x5a==x6a : 1
+DEAL::x5b==x6b_1 : 1
+DEAL::x5c==x6c_1 : 1
+DEAL::x5d==x6d_1 : 1
+DEAL::x5b==x6b_2 : 1
+DEAL::x5c==x6c_2 : 1
+DEAL::x5d==x6d_2 : 1
+DEAL::OK
diff --git a/tests/lac/schur_complement_01.cc b/tests/lac/schur_complement_01.cc
new file mode 100644 (file)
index 0000000..c51c9bd
--- /dev/null
@@ -0,0 +1,256 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Test internal preconditioner and solver options
+
+#include "../tests.h"
+
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/packaged_operation.h>
+#include <deal.II/lac/schur_complement.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+
+#define PRINTME(name, var) \
+  deallog \
+      << "Solution vector: " << name << ": " \
+      << var;
+
+using namespace dealii;
+
+
+int main()
+{
+  initlog();
+  deallog.depth_console(0);
+  deallog << std::setprecision(10);
+
+  // deal.II SparseMatrix
+  {
+
+    deallog << "Schur complement" << std::endl;
+    deallog.push("SC_SparseMatrix");
+
+    {
+      deallog << "SparseMatrix 1" << std::endl;
+
+      /* MATLAB / Gnu Octave code
+
+        clear all;
+        printf("SparseMatrix 1")
+        A = [1,2;3,4]
+        b = [5;6]
+        x = A\b
+
+       */
+
+      const unsigned int rc=1;
+      SparsityPattern sparsity_pattern (rc, rc, 0);
+      sparsity_pattern.compress();
+
+      SparseMatrix<double> A (sparsity_pattern);
+      SparseMatrix<double> B (sparsity_pattern);
+      SparseMatrix<double> C (sparsity_pattern);
+      SparseMatrix<double> D (sparsity_pattern);
+      Vector<double> x (rc);
+      Vector<double> y (rc);
+      Vector<double> f (rc);
+      Vector<double> g (rc);
+      for (unsigned int i=0; i < rc; ++i)
+        {
+          A.diag_element(i) = 1.0*(i+1);
+          B.diag_element(i) = 2.0*(i+1);
+          C.diag_element(i) = 3.0*(i+1);
+          D.diag_element(i) = 4.0*(i+1);
+          f(i) = 5.0*(i+1);
+          g(i) = 6.0*(i+1);
+        }
+
+      const auto lo_A = linear_operator(A);
+      const auto lo_B = linear_operator(B);
+      const auto lo_C = linear_operator(C);
+      const auto lo_D = linear_operator(D);
+
+      SolverControl solver_control_A (100, 1.0e-10);
+      SolverCG< Vector<double> > solver_A (solver_control_A);
+      PreconditionJacobi< SparseMatrix<double> > preconditioner_A;
+      preconditioner_A.initialize(A);
+      const auto lo_A_inv = inverse_operator(lo_A,
+                                             solver_A,
+                                             preconditioner_A);
+
+      const auto lo_S = schur_complement(lo_A_inv,lo_B,
+                                         lo_C,lo_D);
+
+      SolverControl solver_control_S (100, 1.0e-10);
+      SolverCG< Vector<double> > solver_S (solver_control_S);
+      PreconditionJacobi< SparseMatrix<double> > preconditioner_S;
+      preconditioner_S.initialize(D); // Same space as S
+      const auto lo_S_inv = inverse_operator(lo_S,
+                                             solver_S,
+                                             preconditioner_S);
+
+      auto rhs = condense_schur_rhs (lo_A_inv,lo_C,f,g);
+      y = lo_S_inv * rhs; // Solve for y
+      x = postprocess_schur_solution (lo_A_inv,lo_B,y,f);
+
+      PRINTME("x", x);
+      PRINTME("y", y);
+    }
+
+    deallog << "SparseMatrix OK" << std::endl;
+  }
+
+  // deal.II BlockSparseMatrix
+  {
+
+    deallog.push("SC_BlockSparseMatrix");
+
+    {
+      deallog << "BlockSparseMatrix 1" << std::endl;
+
+      /* MATLAB / Gnu Octave code
+
+         clear all;
+         printf("BlockSparseMatrix 1")
+         blks=2;
+         rc=10;
+         for (i=0:rc-1)
+           for (bi=0:blks-1)
+             b(bi*rc+i+1,1) = bi*rc + i;
+             for (bj=0:blks-1)
+               el_i = 1 + i + bi*rc;
+               el_j = 1 + i + bj*rc;
+               A(el_i,el_j) = 2.0*bi + 1.5*bj + (i+1);
+             endfor
+           endfor
+         endfor
+         A
+         b
+         x = A\b
+
+       */
+
+      const unsigned int blks=2;
+      const unsigned int rc=10;
+      BlockSparsityPattern sparsity_pattern;
+      {
+        BlockCompressedSimpleSparsityPattern csp(blks, blks);
+        for (unsigned int bi=0; bi<blks; ++bi)
+          for (unsigned int bj=0; bj<blks; ++bj)
+            csp.block(bi,bj).reinit(rc,rc);
+
+        csp.collect_sizes();
+        sparsity_pattern.copy_from(csp);
+      }
+
+      BlockSparseMatrix<double> A (sparsity_pattern);
+      BlockVector<double> b (blks,rc);
+      for (unsigned int i=0; i<rc; ++i)
+        {
+          for (unsigned int bi=0; bi<blks; ++bi)
+            {
+              b.block(bi)(i) = bi*rc + i;
+              for (unsigned int bj=0; bj<blks; ++bj)
+                A.block(bi,bj).diag_element(i) = 2.0*bi + 1.5*bj + (i+1);
+            }
+        }
+
+      const auto lo_A = linear_operator(A.block(1,1));
+      const auto lo_B = linear_operator(A.block(1,0));
+      const auto lo_C = linear_operator(A.block(0,1));
+      const auto lo_D = linear_operator(A.block(0,0));
+
+      Vector<double> &f = b.block(1);
+      Vector<double> &g = b.block(0);
+
+      BlockVector<double> s (blks,rc);
+      Vector<double> &x = s.block(1);
+      Vector<double> &y = s.block(0);
+
+      SolverControl solver_control_A (100, 1.0e-10);
+      SolverCG< Vector<double> > solver_A (solver_control_A);
+      PreconditionJacobi< SparseMatrix<double> > preconditioner_A;
+      preconditioner_A.initialize(A.block(1,1));
+      const auto lo_A_inv = inverse_operator(lo_A,
+                                             solver_A,
+                                             preconditioner_A);
+
+      const auto lo_S = schur_complement(lo_A_inv,lo_B,
+                                         lo_C,lo_D);
+
+      // Preconditinoed by D
+      {
+        SolverControl solver_control_S (100, 1.0e-10);
+        SolverCG< Vector<double> > solver_S (solver_control_S);
+        PreconditionJacobi< SparseMatrix<double> > preconditioner_S;
+        preconditioner_S.initialize(A.block(0,0)); // Same space as S
+        const auto lo_S_inv = inverse_operator(lo_S,
+                                               solver_S,
+                                               preconditioner_S);
+
+        auto rhs = condense_schur_rhs (lo_A_inv,lo_C,f,g);
+        y = lo_S_inv * rhs; // Solve for y
+        x = postprocess_schur_solution (lo_A_inv,lo_B,y,f);
+
+        PRINTME("x = s.block(1)", x);
+        PRINTME("y = s.block(0)", y);
+      }
+
+      // Preconditinoed by S_approx_inv
+      {
+        const auto lo_A_inv_approx = linear_operator(preconditioner_A);
+        const auto lo_S_approx = schur_complement(lo_A_inv_approx,
+                                                  lo_B,lo_C,lo_D);
+
+        // Setup inner solver: Approximation of inverse of Schur complement
+        IterationNumberControl solver_control_S_approx (1, 1.0e-10); // Perform only a limited number of sweeps
+        SolverCG< Vector<double> > solver_S_approx (solver_control_S_approx);
+        PreconditionJacobi< SparseMatrix<double> > preconditioner_S_approx;
+        preconditioner_S_approx.initialize(A.block(0,0)); // Same space as S
+        const auto lo_S_inv_approx = inverse_operator(lo_S_approx,
+                                                      solver_S_approx,
+                                                      preconditioner_S_approx);
+
+        // Setup outer solver: Exact inverse of Schur complement
+        SolverControl solver_control_S (100, 1.0e-10);
+        SolverCG< Vector<double> > solver_S (solver_control_S);
+        const auto lo_S_inv = inverse_operator(lo_S,
+                                               solver_S,
+                                               lo_S_inv_approx);
+
+        auto rhs = condense_schur_rhs (lo_A_inv,lo_C,f,g);
+        y = lo_S_inv * rhs; // Solve for y
+        x = postprocess_schur_solution (lo_A_inv,lo_B,y,f);
+
+        PRINTME("x = s.block(1)", x);
+        PRINTME("y = s.block(0)", y);
+      }
+
+//      A.print(std::cout);
+//      b.print(std::cout);
+//      s.print(std::cout);
+    }
+
+    deallog << "BlockSparseMatrix OK" << std::endl;
+  }
+
+}
diff --git a/tests/lac/schur_complement_01.with_cxx11=on.output b/tests/lac/schur_complement_01.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..258133e
--- /dev/null
@@ -0,0 +1,98 @@
+
+DEAL::Schur complement
+DEAL:SC_SparseMatrix::SparseMatrix 1
+DEAL:SC_SparseMatrix:cg::Starting value 5.000000000
+DEAL:SC_SparseMatrix:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:cg::Starting value 9.000000000
+DEAL:SC_SparseMatrix:cg:cg::Starting value 4.500000000
+DEAL:SC_SparseMatrix:cg:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:cg::Starting value 4.000000000
+DEAL:SC_SparseMatrix:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix::Solution vector: x: -4.000000000 
+DEAL:SC_SparseMatrix::Solution vector: y: 4.500000000 
+DEAL:SC_SparseMatrix::SparseMatrix OK
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix::BlockSparseMatrix 1
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg::Starting value 46.74398357
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg::Starting value 21.00919703
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 33.27448540
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 45.97490107
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 1.776356839e-15
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 28.96138274
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 1.776356839e-15
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 15.19528992
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 7.010434694
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 1.444357125e-15
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 2.860158624
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 1.146493431e-16
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 1.013942270
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 0.3017508708
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 0.07094905937
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 0.01140752780
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 2.179388120e-18
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 2.200743662e-06
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg::Convergence step 11 value 4.773887375e-11
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg::Starting value 622.4076906
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg::Convergence step 1 value 2.929642751e-14
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix::Solution vector: x = s.block(1): -3.333333333 -6.000000000 -8.666666667 -11.33333333 -14.00000000 -16.66666667 -19.33333333 -22.00000000 -24.66666667 -27.33333333 
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix::Solution vector: y = s.block(0): 8.333333333 11.00000000 13.66666667 16.33333333 19.00000000 21.66666667 24.33333333 27.00000000 29.66666667 32.33333333 
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg::Starting value 46.74398357
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg::Starting value 21.00919703
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 21.00919703
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 17.45229624
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 124.2796216
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 7.105427358e-15
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 17.45229624
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 14.35782789
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 100.0213190
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 12.34050766
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 8.991986820
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 129.0184156
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 2.587940646e-14
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 7.482680505
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 4.974308439
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 107.8383444
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 3.923140398
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 2.361041089
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 71.20654792
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 1.772116736
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 0.9501693577
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 38.95216224
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 0.6799243772
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 0.3167664073
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 17.69780734
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 3.148515104e-15
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 0.2154366229
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 0.08362844279
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 6.528344074
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 2.775557892e-17
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 0.05331266845
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 0.01584789371
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 1.854052160
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 3.565310485e-16
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 0.008945562859
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 0.001552120499
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 0.3527279632
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 6.947697991e-17
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 4.183032022e-07
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 8.661810558e-12
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Starting value 1.882364435e-06
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg:cg::Convergence step 1 value 2.908056841e-26
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg::Convergence step 11 value 4.323638387e-12
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg::Starting value 622.4076906
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix:cg::Convergence step 1 value 7.105427358e-15
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix::Solution vector: x = s.block(1): -3.333333333 -6.000000000 -8.666666667 -11.33333333 -14.00000000 -16.66666667 -19.33333333 -22.00000000 -24.66666667 -27.33333333 
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix::Solution vector: y = s.block(0): 8.333333333 11.00000000 13.66666667 16.33333333 19.00000000 21.66666667 24.33333333 27.00000000 29.66666667 32.33333333 
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix::BlockSparseMatrix OK
diff --git a/tests/lac/schur_complement_02.cc b/tests/lac/schur_complement_02.cc
new file mode 100644 (file)
index 0000000..958faa1
--- /dev/null
@@ -0,0 +1,132 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Test internal preconditioner and solver options
+
+#include "../tests.h"
+
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/packaged_operation.h>
+#include <deal.II/lac/schur_complement.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+
+#define PRINTME(name, var) \
+  deallog \
+      << "RHS vector: " << name << ": " \
+      << var;
+
+using namespace dealii;
+
+
+int main()
+{
+  initlog();
+  deallog.depth_console(0);
+  deallog << std::setprecision(10);
+
+  // deal.II SparseMatrix
+  {
+
+    deallog << "Schur complement" << std::endl;
+    deallog.push("SC_SparseMatrix");
+
+    {
+      deallog << "SparseMatrix 1" << std::endl;
+
+      /* MATLAB / Gnu Octave code
+
+        clear all;
+        printf("SparseMatrix 1")
+        A = [1];
+        B = [2];
+        C = [3];
+        D = [4];
+        y = [6];
+
+        S = D - C*inv(A)*B
+
+        % vmult
+        g1 = S*y
+        % Tvmult
+        g2 = S'*y
+
+        g = [2];
+
+        % vmult_add
+        g3 = S*y + g
+        % Tvmult_add
+        g4 = S'*y + g
+
+       */
+
+      const unsigned int rc=1;
+      SparsityPattern sparsity_pattern (rc, rc, 0);
+      sparsity_pattern.compress();
+
+      SparseMatrix<double> A (sparsity_pattern);
+      SparseMatrix<double> B (sparsity_pattern);
+      SparseMatrix<double> C (sparsity_pattern);
+      SparseMatrix<double> D (sparsity_pattern);
+      Vector<double> y (rc);
+      Vector<double> g (rc);
+      for (unsigned int i=0; i < rc; ++i)
+        {
+          A.diag_element(i) = 1.0*(i+1);
+          B.diag_element(i) = 2.0*(i+1);
+          C.diag_element(i) = 3.0*(i+1);
+          D.diag_element(i) = 4.0*(i+1);
+          y(i) = 6.0*(i+1);
+          g(i) = 2.0*(i+1);
+        }
+
+      const auto lo_A = linear_operator(A);
+      const auto lo_B = linear_operator(B);
+      const auto lo_C = linear_operator(C);
+      const auto lo_D = linear_operator(D);
+
+      SolverControl solver_control_A (100, 1.0e-10);
+      SolverCG< Vector<double> > solver_A (solver_control_A);
+      PreconditionJacobi< SparseMatrix<double> > preconditioner_A;
+      preconditioner_A.initialize(A);
+      const auto lo_A_inv = inverse_operator(lo_A,
+                                             solver_A,
+                                             preconditioner_A);
+
+      const auto lo_S = schur_complement(lo_A_inv,lo_B,
+                                         lo_C,lo_D);
+      const auto lo_S_t = transpose_operator(lo_S);
+
+      const Vector<double> g1 = lo_S*y;
+      const Vector<double> g2 = lo_S_t *y;
+      const Vector<double> g3 = lo_S*y + g;
+      const Vector<double> g4 = lo_S_t *y + g;
+
+      PRINTME("g1",g1);
+      PRINTME("g2",g2);
+      PRINTME("g3",g3);
+      PRINTME("g4",g4);
+    }
+
+    deallog << "SparseMatrix OK" << std::endl;
+  }
+
+}
diff --git a/tests/lac/schur_complement_02.with_cxx11=on.output b/tests/lac/schur_complement_02.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..435b5ff
--- /dev/null
@@ -0,0 +1,16 @@
+
+DEAL::Schur complement
+DEAL:SC_SparseMatrix::SparseMatrix 1
+DEAL:SC_SparseMatrix:cg::Starting value 12.00000000
+DEAL:SC_SparseMatrix:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:cg::Starting value 18.00000000
+DEAL:SC_SparseMatrix:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:cg::Starting value 12.00000000
+DEAL:SC_SparseMatrix:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix:cg::Starting value 18.00000000
+DEAL:SC_SparseMatrix:cg::Convergence step 1 value 0.000000000
+DEAL:SC_SparseMatrix::RHS vector: g1: -12.00000000 
+DEAL:SC_SparseMatrix::RHS vector: g2: -12.00000000 
+DEAL:SC_SparseMatrix::RHS vector: g3: -10.00000000 
+DEAL:SC_SparseMatrix::RHS vector: g4: -10.00000000 
+DEAL:SC_SparseMatrix::SparseMatrix OK
diff --git a/tests/lac/schur_complement_03.cc b/tests/lac/schur_complement_03.cc
new file mode 100644 (file)
index 0000000..4977a61
--- /dev/null
@@ -0,0 +1,429 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2008 - 2015 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE at
+ * the top level of the deal.II distribution.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * Author: Wolfgang Bangerth, Texas A&M University, 2008
+ */
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_ilu.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/packaged_operation.h>
+#include <deal.II/lac/schur_complement.h>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+namespace Step22
+{
+  using namespace dealii;
+  template <int dim>
+  class StokesProblem
+  {
+  public:
+    StokesProblem (const unsigned int degree);
+    void run ();
+  private:
+    void setup_dofs ();
+    void assemble_system ();
+    void solve ();
+    void output_results (const unsigned int refinement_cycle) const;
+    void refine_mesh ();
+    const unsigned int   degree;
+    Triangulation<dim>   triangulation;
+    FESystem<dim>        fe;
+    DoFHandler<dim>      dof_handler;
+    ConstraintMatrix     constraints;
+    BlockSparsityPattern      sparsity_pattern;
+    BlockSparseMatrix<double> system_matrix;
+    BlockVector<double> solution;
+    BlockVector<double> system_rhs;
+  };
+  template <int dim>
+  class BoundaryValues : public Function<dim>
+  {
+  public:
+    BoundaryValues () : Function<dim>(dim+1) {}
+    virtual double value (const Point<dim>   &p,
+                          const unsigned int  component = 0) const;
+    virtual void vector_value (const Point<dim> &p,
+                               Vector<double>   &value) const;
+  };
+  template <int dim>
+  double
+  BoundaryValues<dim>::value (const Point<dim>  &p,
+                              const unsigned int component) const
+  {
+    Assert (component < this->n_components,
+            ExcIndexRange (component, 0, this->n_components));
+    if (component == 0)
+      return (p[0] < 0 ? -1 : (p[0] > 0 ? 1 : 0));
+    return 0;
+  }
+  template <int dim>
+  void
+  BoundaryValues<dim>::vector_value (const Point<dim> &p,
+                                     Vector<double>   &values) const
+  {
+    for (unsigned int c=0; c<this->n_components; ++c)
+      values(c) = BoundaryValues<dim>::value (p, c);
+  }
+  template <int dim>
+  class RightHandSide : public Function<dim>
+  {
+  public:
+    RightHandSide () : Function<dim>(dim+1) {}
+    virtual double value (const Point<dim>   &p,
+                          const unsigned int  component = 0) const;
+    virtual void vector_value (const Point<dim> &p,
+                               Vector<double>   &value) const;
+  };
+  template <int dim>
+  double
+  RightHandSide<dim>::value (const Point<dim>  &/*p*/,
+                             const unsigned int /*component*/) const
+  {
+    return 0;
+  }
+  template <int dim>
+  void
+  RightHandSide<dim>::vector_value (const Point<dim> &p,
+                                    Vector<double>   &values) const
+  {
+    for (unsigned int c=0; c<this->n_components; ++c)
+      values(c) = RightHandSide<dim>::value (p, c);
+  }
+  template <int dim>
+  StokesProblem<dim>::StokesProblem (const unsigned int degree)
+    :
+    degree (degree),
+    triangulation (Triangulation<dim>::maximum_smoothing),
+    fe (FE_Q<dim>(degree+1), dim,
+        FE_Q<dim>(degree), 1),
+    dof_handler (triangulation)
+  {}
+  template <int dim>
+  void StokesProblem<dim>::setup_dofs ()
+  {
+    system_matrix.clear ();
+    dof_handler.distribute_dofs (fe);
+    DoFRenumbering::Cuthill_McKee (dof_handler);
+    std::vector<unsigned int> block_component (dim+1,0);
+    block_component[dim] = 1;
+    DoFRenumbering::component_wise (dof_handler, block_component);
+    {
+      constraints.clear ();
+      FEValuesExtractors::Vector velocities(0);
+      DoFTools::make_hanging_node_constraints (dof_handler,
+                                               constraints);
+      VectorTools::interpolate_boundary_values (dof_handler,
+                                                1,
+                                                BoundaryValues<dim>(),
+                                                constraints,
+                                                fe.component_mask(velocities));
+    }
+    constraints.close ();
+    std::vector<types::global_dof_index> dofs_per_block (2);
+    DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);
+    const unsigned int n_u = dofs_per_block[0],
+                       n_p = dofs_per_block[1];
+    std::cout << "   Number of active cells: "
+              << triangulation.n_active_cells()
+              << std::endl
+              << "   Number of degrees of freedom: "
+              << dof_handler.n_dofs()
+              << " (" << n_u << '+' << n_p << ')'
+              << std::endl;
+    {
+      BlockDynamicSparsityPattern dsp (2,2);
+      dsp.block(0,0).reinit (n_u, n_u);
+      dsp.block(1,0).reinit (n_p, n_u);
+      dsp.block(0,1).reinit (n_u, n_p);
+      dsp.block(1,1).reinit (n_p, n_p);
+      dsp.collect_sizes();
+      DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false);
+      sparsity_pattern.copy_from (dsp);
+    }
+    system_matrix.reinit (sparsity_pattern);
+    solution.reinit (2);
+    solution.block(0).reinit (n_u);
+    solution.block(1).reinit (n_p);
+    solution.collect_sizes ();
+    system_rhs.reinit (2);
+    system_rhs.block(0).reinit (n_u);
+    system_rhs.block(1).reinit (n_p);
+    system_rhs.collect_sizes ();
+  }
+  template <int dim>
+  void StokesProblem<dim>::assemble_system ()
+  {
+    system_matrix=0;
+    system_rhs=0;
+    QGauss<dim>   quadrature_formula(degree+2);
+    FEValues<dim> fe_values (fe, quadrature_formula,
+                             update_values    |
+                             update_quadrature_points  |
+                             update_JxW_values |
+                             update_gradients);
+    const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+    const unsigned int   n_q_points      = quadrature_formula.size();
+    FullMatrix<double>   local_matrix (dofs_per_cell, dofs_per_cell);
+    Vector<double>       local_rhs (dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+    const RightHandSide<dim>          right_hand_side;
+    std::vector<Vector<double> >      rhs_values (n_q_points,
+                                                  Vector<double>(dim+1));
+    const FEValuesExtractors::Vector velocities (0);
+    const FEValuesExtractors::Scalar pressure (dim);
+    std::vector<SymmetricTensor<2,dim> > symgrad_phi_u (dofs_per_cell);
+    std::vector<double>                  div_phi_u   (dofs_per_cell);
+    std::vector<double>                  phi_p       (dofs_per_cell);
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      {
+        fe_values.reinit (cell);
+        local_matrix = 0;
+        local_rhs = 0;
+        right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
+                                          rhs_values);
+        for (unsigned int q=0; q<n_q_points; ++q)
+          {
+            for (unsigned int k=0; k<dofs_per_cell; ++k)
+              {
+                symgrad_phi_u[k] = fe_values[velocities].symmetric_gradient (k, q);
+                div_phi_u[k]     = fe_values[velocities].divergence (k, q);
+                phi_p[k]         = fe_values[pressure].value (k, q);
+              }
+            for (unsigned int i=0; i<dofs_per_cell; ++i)
+              {
+                for (unsigned int j=0; j<=i; ++j)
+                  {
+                    local_matrix(i,j) += (2 * (symgrad_phi_u[i] * symgrad_phi_u[j])
+                                          - div_phi_u[i] * phi_p[j]
+                                          - phi_p[i] * div_phi_u[j]
+                                          + phi_p[i] * phi_p[j])
+                                         * fe_values.JxW(q);
+                  }
+                const unsigned int component_i =
+                  fe.system_to_component_index(i).first;
+                local_rhs(i) += fe_values.shape_value(i,q) *
+                                rhs_values[q](component_i) *
+                                fe_values.JxW(q);
+              }
+          }
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          for (unsigned int j=i+1; j<dofs_per_cell; ++j)
+            local_matrix(i,j) = local_matrix(j,i);
+        cell->get_dof_indices (local_dof_indices);
+        constraints.distribute_local_to_global (local_matrix, local_rhs,
+                                                local_dof_indices,
+                                                system_matrix, system_rhs);
+      }
+  }
+  template <int dim>
+  void StokesProblem<dim>::solve ()
+  {
+    // Linear operators
+    const auto A  = linear_operator(system_matrix.block(0,0));
+    const auto B  = linear_operator(system_matrix.block(0,1));
+    const auto C  = linear_operator(system_matrix.block(1,0));
+    const auto M  = linear_operator(system_matrix.block(1,1)); // Mass matrix stored in this block
+    const auto D0 = null_operator(M);
+
+    // Inverse of A
+    SparseILU<double> preconditioner_A;
+    preconditioner_A.initialize (system_matrix.block(0,0),
+                                 SparseILU<double>::AdditionalData());
+    ReductionControl solver_control_A (system_matrix.block(0,0).m(),
+                                       1e-10, 1e-6);
+    SolverCG<> solver_A (solver_control_A);
+    const auto A_inv = inverse_operator(A, solver_A,
+                                        preconditioner_A);
+
+    // Inverse of mass matrix stored in block "D"
+    SparseILU<double> preconditioner_M;
+    preconditioner_M.initialize (system_matrix.block(1,1),
+                                 SparseILU<double>::AdditionalData());
+    ReductionControl solver_control_M (system_matrix.block(1,1).m(),
+                                       1e-10, 1e-6);
+    SolverCG<> solver_M (solver_control_M);
+    const auto M_inv = inverse_operator(M, solver_M,
+                                        preconditioner_M);
+
+    // Schur complement
+    const auto S = schur_complement(A_inv,B,C,D0);
+
+    // Inverse of Schur complement
+    ReductionControl solver_control_S (system_matrix.block(1,1).m(),
+                                       1e-10, 1e-6);
+    SolverCG<> solver_S (solver_control_S);
+    const auto S_inv = inverse_operator(S,solver_S,M_inv);
+
+    Vector<double> &x = solution.block(0);
+    Vector<double> &y = solution.block(1);
+    const Vector<double> &f = system_rhs.block(0);
+    const Vector<double> &g = system_rhs.block(1);
+    auto rhs = condense_schur_rhs (A_inv,C,f,g);
+    y = S_inv * rhs;
+    x = postprocess_schur_solution (A_inv,B,y,f);
+
+    constraints.distribute (solution);
+    std::cout << "  "
+              << solver_control_S.last_step()
+              << " outer CG Schur complement iterations for pressure"
+              << std::endl;
+  }
+  template <int dim>
+  void
+  StokesProblem<dim>::output_results (const unsigned int refinement_cycle)  const
+  {
+    std::vector<std::string> solution_names (dim, "velocity");
+    solution_names.push_back ("pressure");
+    std::vector<DataComponentInterpretation::DataComponentInterpretation>
+    data_component_interpretation
+    (dim, DataComponentInterpretation::component_is_part_of_vector);
+    data_component_interpretation
+    .push_back (DataComponentInterpretation::component_is_scalar);
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler (dof_handler);
+    data_out.add_data_vector (solution, solution_names,
+                              DataOut<dim>::type_dof_data,
+                              data_component_interpretation);
+    data_out.build_patches ();
+    std::ostringstream filename;
+    filename << "solution-"
+             << Utilities::int_to_string (refinement_cycle, 2)
+             << ".vtk";
+    std::ofstream output (filename.str().c_str());
+    data_out.write_vtk (output);
+  }
+  template <int dim>
+  void
+  StokesProblem<dim>::refine_mesh ()
+  {
+    Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+    FEValuesExtractors::Scalar pressure(dim);
+    KellyErrorEstimator<dim>::estimate (dof_handler,
+                                        QGauss<dim-1>(degree+1),
+                                        typename FunctionMap<dim>::type(),
+                                        solution,
+                                        estimated_error_per_cell,
+                                        fe.component_mask(pressure));
+    GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                     estimated_error_per_cell,
+                                                     0.3, 0.0);
+    triangulation.execute_coarsening_and_refinement ();
+  }
+  template <int dim>
+  void StokesProblem<dim>::run ()
+  {
+    {
+      std::vector<unsigned int> subdivisions (dim, 1);
+//      subdivisions[0] = 4;
+      const Point<dim> bottom_left = (dim == 2 ?
+                                      Point<dim>(-2,-1) :
+                                      Point<dim>(-2,0,-1));
+      const Point<dim> top_right   = (dim == 2 ?
+                                      Point<dim>(2,0) :
+                                      Point<dim>(2,1,0));
+      GridGenerator::subdivided_hyper_rectangle (triangulation,
+                                                 subdivisions,
+                                                 bottom_left,
+                                                 top_right);
+    }
+    for (typename Triangulation<dim>::active_cell_iterator
+         cell = triangulation.begin_active();
+         cell != triangulation.end(); ++cell)
+      for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+        if (cell->face(f)->center()[dim-1] == 0)
+          cell->face(f)->set_all_boundary_ids(1);
+    triangulation.refine_global (4-dim);
+    for (unsigned int refinement_cycle = 0; refinement_cycle<2;
+         ++refinement_cycle)
+      {
+        std::cout << "Refinement cycle " << refinement_cycle << std::endl;
+        if (refinement_cycle > 0)
+          refine_mesh ();
+        setup_dofs ();
+        std::cout << "   Assembling..." << std::endl << std::flush;
+        assemble_system ();
+        std::cout << "   Solving..." << std::flush;
+        solve ();
+        //output_results (refinement_cycle);
+        std::cout << std::endl;
+      }
+  }
+}
+int main ()
+{
+  try
+    {
+      using namespace dealii;
+      using namespace Step22;
+      StokesProblem<2> flow_problem(1);
+      flow_problem.run ();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  return 0;
+}
diff --git a/tests/lac/schur_complement_03.with_cxx11=on.output b/tests/lac/schur_complement_03.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..2b58626
--- /dev/null
@@ -0,0 +1,12 @@
+Refinement cycle 0
+   Number of active cells: 16
+   Number of degrees of freedom: 187 (162+25)
+   Assembling...
+   Solving...  9 outer CG Schur complement iterations for pressure
+
+Refinement cycle 1
+   Number of active cells: 40
+   Number of degrees of freedom: 441 (386+55)
+   Assembling...
+   Solving...  11 outer CG Schur complement iterations for pressure
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.