fgmres

Source: tests/iterative/fgmres.cpp

// ========================================================================= //
// This file is part of MyraMath, copyright (c) 2014-2019 by Ryan A Chilton  //
// and distributed by MyraCore, LLC. See LICENSE.txt for license terms.      //
// ========================================================================= //

// Containers.
#include <myramath/utility/Number.h>
#include <myramath/dense/Matrix.h>
#include <myramath/dense/MatrixRange.h>
#include <myramath/dense/Vector.h>
#include <myramath/dense/VectorRange.h>
#include <myramath/sparse/Pattern.h>
#include <myramath/sparse/PatternRange.h>
#include <myramath/sparse/SparseMatrix.h>
#include <myramath/sparse/SparseMatrixRange.h>
#include <myramath/iterative/ILUSolver.h>

// Algorithms.
#include <myramath/utility/cumsum.h>
#include <myramath/dense/inverse.h>
#include <myramath/dense/euclidean.h>
#include <myramath/sparse/gemm.h>
#include <myramath/sparse/gemv.h>
#include <myramath/sparse/laplacian2.h>
#include <myramath/sparse/frobenius.h>

// Iterative solve/action stuff.
#include <myramath/iterative/gmres.h>
#include <myramath/iterative/fgmres.h>
#include <myramath/iterative/Action.h>
#include <myramath/iterative/GemmAction.h>
#include <myramath/iterative/IdentityAction.h>
#include <myramath/iterative/SolveAction.h>
#include <myramath/iterative/UserAction.h>
#include <myramath/iterative/ILUSolver.h>
#include <myramath/iterative/ProductAction.h>

// Reporting.
#include <myramath/utility/Timer.h>
#include <myramath/utility/stlprint.h>
#include <tests/myratest.h>

using namespace myra;
using namespace myra_stlprint;

namespace {

// Action that solves by the L factor of an ILUSolver.
template<class InstantiatedNumber> class solveL_Action
  {
  public:
    typedef InstantiatedNumber Number;
    solveL_Action(const ILUSolver<Number>& in_solver)
      : solver(in_solver) { }
    std::pair<int,int> size() const
      { int N = solver.size(); return std::pair<int,int>(N,N); }
    void multiply(const CMatrixRange<Number>& B, const MatrixRange<Number>& X) const
      { X.assign(B); solver.solveL(X,'L','N'); }
  private:
    const ILUSolver<Number>& solver;
  };

// Action that solves by the U factor of an ILUSolver.
template<class InstantiatedNumber> class solveU_Action
  {
  public:
    typedef InstantiatedNumber Number;
    solveU_Action(const ILUSolver<Number>& in_solver)
      : solver(in_solver) { }
    std::pair<int,int> size() const
      { int N = solver.size(); return std::pair<int,int>(N,N); }
    void multiply(const CMatrixRange<Number>& B, const MatrixRange<Number>& X) const
      { X.assign(B); solver.solveU(X,'L','N'); }
  private:
    const ILUSolver<Number>& solver;
  };

// Action that applies gmres()
template<class InstantiatedNumber> class gmres_Action
  {
  public:

    // Required typedef.
    typedef InstantiatedNumber Number;
    typedef typename ReflectPrecision<Number>::type Precision;

    gmres_Action(const Action<Number>& in_A, Precision in_tol, int in_iter)
      : A(in_A), tol(in_tol), iter(in_iter) { M = make_IdentityAction<Number>(A.size().first); }

    std::pair<int,int> size() const
      { return A.size(); }

    void multiply (const CMatrixRange<Number>& B, const MatrixRange<Number>& X) const
      { X.zero(); auto output = gmres(M,A,B,X,tol,iter); }

  private:

    // Captured environment for inner gmres()
    Action<Number> M;
    Action<Number> A;
    Precision tol;
    int iter;
    
  };

template<class Number> void test(int I, int J, typename ReflectPrecision<Number>::type tolerance)
  {
  typedef typename ReflectPrecision<Number>::type Precision;
  myra::out() << typestring<Number>() << std::endl;
  int N = I*J;
  auto A = SparseMatrix<Number>::random(stencil2(I,J));
  Precision sigma = 2.0;
  for (int n = 0; n < N; ++n)
    A(n,n) += sigma;
  auto y = Vector<Number>::random(N);
  auto b = A*y;
  ILUSolver<Number> solver(A);

  // Instantiate actions for subsequent test cases.
  auto a = make_GemmAction(A);
  auto i = make_IdentityAction<Number>(N);
  auto z = make_SolveAction(solver);
  auto l = make_UserAction(solveL_Action<Number>(solver));
  auto u = make_UserAction(solveU_Action<Number>(solver));
  auto g = make_UserAction(gmres_Action<Number>(z*a,1.0e-2,20));

  // Solve using gmres(), left ILU.
    {
    auto x = Vector<Number>::zeros(N);  
    auto result = gmres(z,a,b,x,tolerance,500);
    Precision error = euclidean(A*x-b) / euclidean(b);
    myra::out() << "gmres, left ILU:" << std::endl;
    myra::out() << "  |Ax-b| = " << error << std::endl;
    myra::out() << "  iteration count = " << result.history.size() << std::endl;
    myra::out() << std::endl;
    REQUIRE(error < tolerance*100);
    }

  // Solve using fgmres(), left ILU.
    {
    auto x = Vector<Number>::zeros(N);  
    auto result = fgmres(z,a,i,b,x,tolerance,500);
    Precision error = euclidean(A*x-b) / euclidean(b);
    myra::out() << "fgmres, left ILU:" << std::endl;
    myra::out() << "  |Ax-b| = " << error << std::endl;
    myra::out() << "  iteration count = " << result.history.size() << std::endl;
    myra::out() << std::endl;
    REQUIRE(error < tolerance*100);
    }

  // Solve using fgmres(), right ILU.
    {
    auto x = Vector<Number>::zeros(N);  
    auto result = fgmres(i,a,z,b,x,tolerance,500);
    Precision error = euclidean(A*x-b) / euclidean(b);
    myra::out() << "fgmres, right ILU:" << std::endl;
    myra::out() << "  |Ax-b| = " << error << std::endl;
    myra::out() << "  iteration count = " << result.history.size() << std::endl;
    myra::out() << std::endl;
    REQUIRE(error < tolerance*100);
    }

  // Solve using fgmres(), split ILU.
    {
    auto x = Vector<Number>::zeros(N);
    auto result = fgmres(l,a,u,b,x,tolerance,500);
    Precision error = euclidean(A*x-b) / euclidean(b);
    myra::out() << "fgmres, split ILU:" << std::endl;
    myra::out() << "  |Ax-b| = " << error << std::endl;
    myra::out() << "  iteration count = " << result.history.size() << std::endl;
    myra::out() << std::endl;
    REQUIRE(error < tolerance*100);
    }

  // Solve using restarted gmres(), left ILU. Expected to fail, restarts impede convergence.
    {
    auto x = Vector<Number>::zeros(N);  
    auto result = gmres(z,a,b,x,tolerance,20,20);
    Precision error = euclidean(A*x-b) / euclidean(b);
    myra::out() << "restarted gmres, left ILU:" << std::endl;
    myra::out() << "  |Ax-b| = " << error << std::endl;
    myra::out() << "  iteration count = " << result.history.size() << std::endl;
    myra::out() << std::endl;
    // REQUIRE(error < tolerance*100); // will be big
    }

  // Solve using fgmres() with M1=ILUSolver and M2=inaccurate gmres(M1*A). Expected to work.
    {
    auto x = Vector<Number>::zeros(N);
    auto result = fgmres(z,a,g,b,x,tolerance,500);
    Precision error = euclidean(A*x-b) / euclidean(b);
    myra::out() << "fgmres, left ILU and right inner gmres:" << std::endl;
    myra::out() << "  |Ax-b| = " << error << std::endl;
    myra::out() << "  iteration count = " << result.history.size() << std::endl;
    myra::out() << std::endl;
    REQUIRE(error < tolerance*100);
    }

  // Solve using gmres() with M1=ILUSolver and M2=inaccurate gmres(M1*A). This is not 
  // expected to work, because the preconditioner is effectively changing from one outer 
  // gmres() iteration to the next, due to the limits imposed on inner gmres() tolerance 
  // and iteration count. It will converge, but the final solution will not be correct.
    {
    auto x = Vector<Number>::zeros(N);
    auto result = gmres(z,a,g,b,x,tolerance,500);
    Precision error = euclidean(A*x-b) / euclidean(b);
    myra::out() << "gmres, left ILU and right inner gmres:" << std::endl;
    myra::out() << "  |Ax-b| = " << error << std::endl;
    myra::out() << "  iteration count = " << result.history.size() << std::endl;
    myra::out() << std::endl;
    // REQUIRE(error < tolerance*100); // will be big
    }

  }

} // namespace

ADD_TEST("fgmres","[iterative]")
  {
  test<NumberD>(20,20,1.0e-8);
  test<NumberZ>(20,20,1.0e-8);
  }

ADD_TEST("fgmres_big","[iterative][.]")
  {
  test<NumberD>(100,50,1.0e-8);
  test<NumberZ>(100,50,1.0e-8);
  }
  

Results: [PASS]

double
gmres, left ILU:
  |Ax-b| = 4.37427e-09
  iteration count = 19

fgmres, left ILU:
  |Ax-b| = 4.37427e-09
  iteration count = 19

fgmres, right ILU:
  |Ax-b| = 4.0265e-09
  iteration count = 19

fgmres, split ILU:
  |Ax-b| = 4.14965e-09
  iteration count = 19

restarted gmres, left ILU:
  |Ax-b| = 4.37427e-09
  iteration count = 19

fgmres, left ILU and right inner gmres:
  |Ax-b| = 1.76921e-09
  iteration count = 5

gmres, left ILU and right inner gmres:
  |Ax-b| = 0.00393714
  iteration count = 5

std::complex<double>
gmres, left ILU:
  |Ax-b| = 7.80565e-09
  iteration count = 27

fgmres, left ILU:
  |Ax-b| = 7.80565e-09
  iteration count = 27

fgmres, right ILU:
  |Ax-b| = 7.10404e-09
  iteration count = 27

fgmres, split ILU:
  |Ax-b| = 7.75107e-09
  iteration count = 27

restarted gmres, left ILU:
  |Ax-b| = 1.02763e-08
  iteration count = 30

fgmres, left ILU and right inner gmres:
  |Ax-b| = 2.47425e-09
  iteration count = 5

gmres, left ILU and right inner gmres:
  |Ax-b| = 0.00028365
  iteration count = 5

Go back to Summary of /test programs.