Add a test to instrument Linear Operator.

include/deal.II/lac/linear_operator.h

@@ -94,7 +94,7 @@ identity_operator(const LinearOperator<Range, Domain, Payload> &);
  *
  * But, in contrast to a usual matrix object, the domain and range of the
  * linear operator are also bound to the LinearOperator class on the type
- * level. Because of this, <code>LinearOperator <Range, Domain></code> has two
+ * level. Because of this, `LinearOperator<Range, Domain>` has two
  * additional function objects
  * @code
  *   std::function<void(Range &, bool)> reinit_range_vector;
@@ -169,9 +169,31 @@ identity_operator(const LinearOperator<Range, Domain, Payload> &);
  * for linear operators have been provided within the respective
  * TrilinosWrappers (and, in the future, PETScWrappers) namespaces.
  *
- * @note The step-20 tutorial program has a detailed usage example of the
+ * <h3> Examples of use </h3>
+ * The step-20 tutorial program has a detailed usage example of the
  * LinearOperator class.
  *
+ * <h3> Instrumenting operations </h3>
+ * It is sometimes useful to know when functions are called, or to inject
+ * additional operations. In such cases, what one wants is to replace, for
+ * example, the `vmult` object of this class with one that does the additional
+ * operations and then calls what was originally supposed to happen. This
+ * can be done with commands such as the following:
+ * @code
+ *    auto A_inv  = inverse_operator(A, solver_A, preconditioner_A);
+ *    A_inv.vmult = [base_vmult = A_inv.vmult](Vector<double>       &dst,
+ *                                             const Vector<double> &src) {
+ *      std::cout << "Calling A_inv.vmult()" << std::endl;
+ *      base_vmult(dst, src);
+ *    };
+ * @endcode
+ * Here, we replace `A_inv.vmult` with a lambda function that first captures
+ * the previous value of `A_inv.vmult` and stores it in the `base_vmult`
+ * object. The newly installed `A_inv.vmult` function then first outputs some
+ * status information, and then calls the original functionality.
+ *
+ * This approach works for all of the other function objects mentioned above
+ * as well.
  *
  * @ingroup LAOperators
  */

tests/lac/schur_complement_01_instrumented_vmult.cc

@@ -0,0 +1,275 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 - 2020 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Test internal preconditioner and solver options
+
+#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/precondition.h>
+#include <deal.II/lac/schur_complement.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include "../tests.h"
+
+#define PRINTME(name, var) \
+  deallog << "Solution vector: " << name << ": " << var << std::endl;
+
+
+
+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(1, 1.0e-10, false, false);
+      SolverCG<Vector<double>> solver_A(solver_control_A);
+      PreconditionJacobi<SparseMatrix<double>> preconditioner_A;
+      preconditioner_A.initialize(A);
+      auto lo_A_inv  = inverse_operator(lo_A, solver_A, preconditioner_A);
+      lo_A_inv.vmult = [base_vmult =
+                          lo_A_inv.vmult](Vector<double>       &dst,
+                                          const Vector<double> &src) {
+        deallog << "Calling lo_A_inv.vmult()" << std::endl;
+        base_vmult(dst, src);
+      };
+
+      const auto lo_S = schur_complement(lo_A_inv, lo_B, lo_C, lo_D);
+
+      SolverControl            solver_control_S(1, 1.0e-10, false, false);
+      SolverCG<Vector<double>> solver_S(solver_control_S);
+      PreconditionJacobi<SparseMatrix<double>> preconditioner_S;
+      preconditioner_S.initialize(D); // Same space as S
+      auto lo_S_inv  = inverse_operator(lo_S, solver_S, preconditioner_S);
+      lo_S_inv.vmult = [base_vmult =
+                          lo_S_inv.vmult](Vector<double>       &dst,
+                                          const Vector<double> &src) {
+        deallog << "Calling lo_S_inv.vmult()" << std::endl;
+        base_vmult(dst, src);
+      };
+
+      auto rhs = condense_schur_rhs(lo_A_inv, lo_C, f, g);
+      check_solver_within_range(y = lo_S_inv * rhs, // Solve for y
+                                solver_control_S.last_step(),
+                                1,
+                                1);
+
+      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;
+      {
+        BlockDynamicSparsityPattern 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(1, 1.0e-10, false, false);
+      SolverCG<Vector<double>> solver_A(solver_control_A);
+      PreconditionJacobi<SparseMatrix<double>> preconditioner_A;
+      preconditioner_A.initialize(A.block(1, 1));
+      auto lo_A_inv  = inverse_operator(lo_A, solver_A, preconditioner_A);
+      lo_A_inv.vmult = [base_vmult =
+                          lo_A_inv.vmult](Vector<double>       &dst,
+                                          const Vector<double> &src) {
+        deallog << "Calling lo_A_inv.vmult()" << std::endl;
+        base_vmult(dst, src);
+      };
+
+      const auto lo_S = schur_complement(lo_A_inv, lo_B, lo_C, lo_D);
+
+      // Preconditinoed by D
+      {
+        SolverControl            solver_control_S(11, 1.0e-10, false, false);
+        SolverCG<Vector<double>> solver_S(solver_control_S);
+        PreconditionJacobi<SparseMatrix<double>> preconditioner_S;
+        preconditioner_S.initialize(A.block(0, 0)); // Same space as S
+        auto lo_S_inv  = inverse_operator(lo_S, solver_S, preconditioner_S);
+        lo_S_inv.vmult = [base_vmult =
+                            lo_S_inv.vmult](Vector<double>       &dst,
+                                            const Vector<double> &src) {
+          deallog << "Calling lo_S_inv.vmult()" << std::endl;
+          base_vmult(dst, src);
+        };
+
+        auto rhs = condense_schur_rhs(lo_A_inv, lo_C, f, g);
+        check_solver_within_range(y = lo_S_inv * rhs, // Solve for y
+                                  solver_control_S.last_step(),
+                                  11,
+                                  11);
+
+        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, false, false); // 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(11, 1.0e-10, false, false);
+        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);
+        check_solver_within_range(y = lo_S_inv * rhs, // Solve for y
+                                  solver_control_S.last_step(),
+                                  11,
+                                  11);
+
+        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;
+  }
+}

tests/lac/schur_complement_01_instrumented_vmult.output

@@ -0,0 +1,18 @@
+
+DEAL::Schur complement
+DEAL:SC_SparseMatrix::SparseMatrix 1
+DEAL:SC_SparseMatrix::Solver stopped within 1 - 1 iterations
+DEAL:SC_SparseMatrix::Calling lo_A_inv.vmult()
+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::Solver stopped within 11 - 11 iterations
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix::Calling lo_A_inv.vmult()
+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::Solver stopped within 11 - 11 iterations
+DEAL:SC_SparseMatrix:SC_BlockSparseMatrix::Calling lo_A_inv.vmult()
+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