From: Wolfgang Bangerth Date: Fri, 24 Nov 2023 23:46:22 +0000 (-0700) Subject: Add a test to instrument Linear Operator. X-Git-Tag: relicensing~280^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c35248141b0937024a5a7a58cd5c572ea212fa49;p=dealii.git Add a test to instrument Linear Operator. --- diff --git a/tests/lac/schur_complement_01_instrumented_vmult.cc b/tests/lac/schur_complement_01_instrumented_vmult.cc new file mode 100644 index 0000000000..19ccfc750a --- /dev/null +++ b/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 +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#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 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(1, 1.0e-10, false, false); + SolverCG> solver_A(solver_control_A); + PreconditionJacobi> 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 &dst, + const Vector &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> solver_S(solver_control_S); + PreconditionJacobi> 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 &dst, + const Vector &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 A(sparsity_pattern); + BlockVector 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 &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(1, 1.0e-10, false, false); + SolverCG> solver_A(solver_control_A); + PreconditionJacobi> 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 &dst, + const Vector &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> solver_S(solver_control_S); + PreconditionJacobi> 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 &dst, + const Vector &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> solver_S_approx(solver_control_S_approx); + PreconditionJacobi> 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> 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; + } +} diff --git a/tests/lac/schur_complement_01_instrumented_vmult.output b/tests/lac/schur_complement_01_instrumented_vmult.output new file mode 100644 index 0000000000..8f6c9a9d6f --- /dev/null +++ b/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