From 8e434c3d01a02595f89b61dd75473c77047faf44 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Mon, 7 Dec 2015 09:01:21 +0100 Subject: [PATCH] Implementation of Schur complement linear operator and auxillary functions 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 | 14 + include/deal.II/lac/schur_complement.h | 515 ++++++++++++++++++ tests/lac/linear_operator_01b.cc | 168 ++++++ .../linear_operator_01b.with_cxx11=on.output | 33 ++ tests/lac/schur_complement_01.cc | 256 +++++++++ .../schur_complement_01.with_cxx11=on.output | 98 ++++ tests/lac/schur_complement_02.cc | 132 +++++ .../schur_complement_02.with_cxx11=on.output | 16 + tests/lac/schur_complement_03.cc | 429 +++++++++++++++ .../schur_complement_03.with_cxx11=on.output | 12 + 10 files changed, 1673 insertions(+) create mode 100644 include/deal.II/lac/schur_complement.h create mode 100644 tests/lac/linear_operator_01b.cc create mode 100644 tests/lac/linear_operator_01b.with_cxx11=on.output create mode 100644 tests/lac/schur_complement_01.cc create mode 100644 tests/lac/schur_complement_01.with_cxx11=on.output create mode 100644 tests/lac/schur_complement_02.cc create mode 100644 tests/lac/schur_complement_02.with_cxx11=on.output create mode 100644 tests/lac/schur_complement_03.cc create mode 100644 tests/lac/schur_complement_03.with_cxx11=on.output diff --git a/doc/news/changes.h b/doc/news/changes.h index c903556dc3..22407a5a05 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -215,6 +215,20 @@ inconvenience this causes. (Luca Heltai, 2015/12/13) +
  • 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 + tests/lac/schur_complement_01.cc. + The solution of a multi-component problem (namely step-22) using the + schur_complement can be found in + tests/lac/schur_complement_03.cc . +
    + (Jean-Paul Pelteret, Matthias Maier, Martin Kronbichler, 2015/12/07) +
    • +
  • 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 index 0000000000..45b2c17de4 --- /dev/null +++ b/include/deal.II/lac/schur_complement.h @@ -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 +#include +#include +#include +#include + +#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 + + // 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 + * tests/lac/schur_complement_01.cc . + * The solution of a multi-component problem (namely step-22) using the + * schur_complement can be found in + * tests/lac/schur_complement_03.cc . + * + * @see + * @ref GlossBlockLA "Block (linear algebra)" + * @author Jean-Paul Pelteret, Matthias Maier, Martin Kronbichler, 2015 + * + * @ingroup LAOperators + */ +template +LinearOperator +schur_complement(const LinearOperator &A_inv, + const LinearOperator &B, + const LinearOperator &C, + const LinearOperator &D) +{ + LinearOperator 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 vector_memory_f; + static GrowingVectorMemory vector_memory_g; + static GrowingVectorMemory 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 vector_memory_f; + static GrowingVectorMemory vector_memory_g; + static GrowingVectorMemory 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 +PackagedOperation +condense_schur_rhs (const LinearOperator &A_inv, + const LinearOperator &C, + const Range_1 &f, + const Range_2 &g) +{ + PackagedOperation 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 vector_memory_f; + static GrowingVectorMemory 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 +PackagedOperation +postprocess_schur_solution (const LinearOperator &A_inv, + const LinearOperator &B, + const Domain_2 &y, + const Range_1 &f) +{ + PackagedOperation 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 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 index 0000000000..b21bf8697c --- /dev/null +++ b/tests/lac/linear_operator_01b.cc @@ -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 +#include +#include +#include +#include + +#include +#include + +#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 A (sparsity_pattern); + Vector 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 > solver_A_1(solver_control_A_1); + PreconditionJacobi< SparseMatrix > 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 > solver_A_2(solver_control_A_2); + PreconditionJacobi< SparseMatrix > 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 x1 = lo_A*b; + const Vector x2 = lo_A_t *b; + const Vector x3 = lo_A_inv*b; + const Vector x4a = lo_A_inv_t *b; + const Vector 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 B (sparsity_pattern); + SparseMatrix C (sparsity_pattern); + SparseMatrix 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 x5 = S*b; + const Vector x6a = S_t_1*b; + const Vector 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 x5a = lo_B*b; + const Vector x5b = lo_A_inv*x5a; + const Vector x5c = lo_C*x5b; + const Vector x5d = lo_D*b - x5c; + const Vector x6a = lo_C_t *b; + const Vector x6b_1 = lo_A_inv_t *x6a; // using transpose of inverse + const Vector x6c_1 = lo_B_t *x6b_1; + const Vector x6d_1 = lo_D_t *b - x6c_1; + const Vector x6b_2 = lo_A_t_inv*x6a; // using inverse of transpose + const Vector x6c_2 = lo_B_t *x6b_2; + const Vector 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 index 0000000000..178aac01d7 --- /dev/null +++ b/tests/lac/linear_operator_01b.with_cxx11=on.output @@ -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 index 0000000000..c51c9bdd7c --- /dev/null +++ b/tests/lac/schur_complement_01.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#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 A (sparsity_pattern); + SparseMatrix B (sparsity_pattern); + SparseMatrix C (sparsity_pattern); + SparseMatrix D (sparsity_pattern); + Vector x (rc); + Vector y (rc); + Vector f (rc); + Vector 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 > solver_A (solver_control_A); + PreconditionJacobi< SparseMatrix > 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 > solver_S (solver_control_S); + PreconditionJacobi< SparseMatrix > 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 A (sparsity_pattern); + BlockVector b (blks,rc); + for (unsigned int i=0; i &f = b.block(1); + Vector &g = b.block(0); + + BlockVector s (blks,rc); + Vector &x = s.block(1); + Vector &y = s.block(0); + + SolverControl solver_control_A (100, 1.0e-10); + SolverCG< Vector > solver_A (solver_control_A); + PreconditionJacobi< SparseMatrix > 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 > solver_S (solver_control_S); + PreconditionJacobi< SparseMatrix > 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 > solver_S_approx (solver_control_S_approx); + PreconditionJacobi< SparseMatrix > 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 > 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 index 0000000000..258133e422 --- /dev/null +++ b/tests/lac/schur_complement_01.with_cxx11=on.output @@ -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 index 0000000000..958faa14fd --- /dev/null +++ b/tests/lac/schur_complement_02.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#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 A (sparsity_pattern); + SparseMatrix B (sparsity_pattern); + SparseMatrix C (sparsity_pattern); + SparseMatrix D (sparsity_pattern); + Vector y (rc); + Vector 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 > solver_A (solver_control_A); + PreconditionJacobi< SparseMatrix > 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 g1 = lo_S*y; + const Vector g2 = lo_S_t *y; + const Vector g3 = lo_S*y + g; + const Vector 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 index 0000000000..435b5ff669 --- /dev/null +++ b/tests/lac/schur_complement_02.with_cxx11=on.output @@ -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 index 0000000000..4977a610a2 --- /dev/null +++ b/tests/lac/schur_complement_03.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +namespace Step22 +{ + using namespace dealii; + template + 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 triangulation; + FESystem fe; + DoFHandler dof_handler; + ConstraintMatrix constraints; + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + BlockVector solution; + BlockVector system_rhs; + }; + template + class BoundaryValues : public Function + { + public: + BoundaryValues () : Function(dim+1) {} + virtual double value (const Point &p, + const unsigned int component = 0) const; + virtual void vector_value (const Point &p, + Vector &value) const; + }; + template + double + BoundaryValues::value (const Point &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 + void + BoundaryValues::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c=0; cn_components; ++c) + values(c) = BoundaryValues::value (p, c); + } + template + class RightHandSide : public Function + { + public: + RightHandSide () : Function(dim+1) {} + virtual double value (const Point &p, + const unsigned int component = 0) const; + virtual void vector_value (const Point &p, + Vector &value) const; + }; + template + double + RightHandSide::value (const Point &/*p*/, + const unsigned int /*component*/) const + { + return 0; + } + template + void + RightHandSide::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c=0; cn_components; ++c) + values(c) = RightHandSide::value (p, c); + } + template + StokesProblem::StokesProblem (const unsigned int degree) + : + degree (degree), + triangulation (Triangulation::maximum_smoothing), + fe (FE_Q(degree+1), dim, + FE_Q(degree), 1), + dof_handler (triangulation) + {} + template + void StokesProblem::setup_dofs () + { + system_matrix.clear (); + dof_handler.distribute_dofs (fe); + DoFRenumbering::Cuthill_McKee (dof_handler); + std::vector 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(), + constraints, + fe.component_mask(velocities)); + } + constraints.close (); + std::vector 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 + void StokesProblem::assemble_system () + { + system_matrix=0; + system_rhs=0; + QGauss quadrature_formula(degree+2); + FEValues 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 local_matrix (dofs_per_cell, dofs_per_cell); + Vector local_rhs (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + const RightHandSide right_hand_side; + std::vector > rhs_values (n_q_points, + Vector(dim+1)); + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + std::vector > symgrad_phi_u (dofs_per_cell); + std::vector div_phi_u (dofs_per_cell); + std::vector phi_p (dofs_per_cell); + typename DoFHandler::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; qget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (local_matrix, local_rhs, + local_dof_indices, + system_matrix, system_rhs); + } + } + template + void StokesProblem::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 preconditioner_A; + preconditioner_A.initialize (system_matrix.block(0,0), + SparseILU::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 preconditioner_M; + preconditioner_M.initialize (system_matrix.block(1,1), + SparseILU::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 &x = solution.block(0); + Vector &y = solution.block(1); + const Vector &f = system_rhs.block(0); + const Vector &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 + void + StokesProblem::output_results (const unsigned int refinement_cycle) const + { + std::vector solution_names (dim, "velocity"); + solution_names.push_back ("pressure"); + std::vector + data_component_interpretation + (dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation + .push_back (DataComponentInterpretation::component_is_scalar); + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, solution_names, + DataOut::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 + void + StokesProblem::refine_mesh () + { + Vector estimated_error_per_cell (triangulation.n_active_cells()); + FEValuesExtractors::Scalar pressure(dim); + KellyErrorEstimator::estimate (dof_handler, + QGauss(degree+1), + typename FunctionMap::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 + void StokesProblem::run () + { + { + std::vector subdivisions (dim, 1); +// subdivisions[0] = 4; + const Point bottom_left = (dim == 2 ? + Point(-2,-1) : + Point(-2,0,-1)); + const Point top_right = (dim == 2 ? + Point(2,0) : + Point(2,1,0)); + GridGenerator::subdivided_hyper_rectangle (triangulation, + subdivisions, + bottom_left, + top_right); + } + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + for (unsigned int f=0; f::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 index 0000000000..2b58626cbb --- /dev/null +++ b/tests/lac/schur_complement_03.with_cxx11=on.output @@ -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 + -- 2.39.5