(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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ }
+}
+
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ }
+
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ }
+
+}
--- /dev/null
+
+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
--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * 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;
+}
--- /dev/null
+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
+
projects / dealii.git / commitdiff