doc/html/test_gmres.html

 // ========================================================================= //
 // 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/Action.h>
 #include <myramath/iterative/GemmAction.h>
 #include <myramath/iterative/IdentityAction.h>
 #include <myramath/iterative/SolveAction.h>
 #include <myramath/iterative/UserAction.h>

 // Reporting.
 #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;
   };

 template<class Number> void test(int X, int Y, typename ReflectPrecision<Number>::type tolerance)
   {
   typedef typename ReflectPrecision<Number>::type Precision;
   myra::out() << typestring<Number>() << std::endl;
   // Construct 2D stencil with random values.
   auto A = SparseMatrix<Number>::random(stencil2(X,Y));
   // Add some mass to it, to make it a little better conditioned.
   int N = X*Y;
   Precision sigma = 2.0;
   for (int n = 0; n < N; ++n)
     A(n,n) += sigma;
   A /= frobenius(A);
   // Construct random x/b vectors.
   auto x = Vector<Number>::random(N);
   auto b = A*x;
   // Solve using gmres() - no preconditioner.
     {
     auto m = make_IdentityAction<Number>(N);
     auto a = make_GemmAction(A);
     x = Vector<Number>::zeros(N);
     auto result = gmres(m,a,b,x,tolerance);
     Precision error = euclidean(A*x-b) / euclidean(b);
     myra::out() << "  No preconditioner:" << std::endl;
     myra::out() << "  |Ax-b| = " << error << std::endl;
     myra::out() << "  iteration count = " << result.history.size() << std::endl;
     //myra::out() << "  history = " << result.history << std::endl;
     myra::out() << std::endl;
     REQUIRE(error < tolerance*100);
     }
   // Solve using gmres() - incomplete LU preconditioner for M1.
     {
     ILUSolver<Number> M1(A);
     auto m1 = make_SolveAction(M1);
     auto m2 = make_IdentityAction<Number>(N);
     auto a = make_GemmAction(A);
     x = Vector<Number>::zeros(N);
     auto result = gmres(m1,a,m2,b,x,tolerance);
     Precision error = euclidean(A*x-b) / euclidean(b);
     myra::out() << "  ILU preconditioner (left):" << std::endl;
     myra::out() << "  |Ax-b| = " << error << std::endl;
     myra::out() << "  iteration count = " << result.history.size() << std::endl;
     //myra::out() << "  history = " << result.history << std::endl;
     myra::out() << std::endl;
     REQUIRE(error < tolerance*100);
     }

   // Solve using gmres() - incomplete LU preconditioner for M2.
     {
     ILUSolver<Number> M2(A);
     auto m1 = make_IdentityAction<Number>(N);
     auto m2 = make_SolveAction(M2);
     auto a = make_GemmAction(A);
     x = Vector<Number>::zeros(N);
     auto result = gmres(m1,a,m2,b,x,tolerance);
     Precision error = euclidean(A*x-b) / euclidean(b);
     myra::out() << "  ILU preconditioner (right):" << std::endl;
     myra::out() << "  |Ax-b| = " << error << std::endl;
     myra::out() << "  iteration count = " << result.history.size() << std::endl;
     //myra::out() << "  history = " << result.history << std::endl;
     myra::out() << std::endl;
     REQUIRE(error < tolerance*100);
     }

   // Solve using gmres() - incomplete LU preconditioner split into M1 (L) and M2 (U).
     {
     ILUSolver<Number> M(A);
     auto m1 = make_UserAction( solveL_Action<Number>(M) );
     auto m2 = make_UserAction( solveU_Action<Number>(M) );
     auto a = make_GemmAction(A);
     x = Vector<Number>::zeros(N);
     auto result = gmres(m1,a,m2,b,x,tolerance);
     Precision error = euclidean(A*x-b) / euclidean(b);
     myra::out() << "  ILU preconditioner (split):" << 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);
     }

   }

 } // namespace

 ADD_TEST("gmres","[iterative]")
   {
   test<NumberS>(20,20,1.0e-3f);
   test<NumberD>(20,20,1.0e-6);
   test<NumberC>(20,20,1.0e-3f);
   test<NumberZ>(20,20,1.0e-6);
   }

euclidean.h
Routines for computing euclidean norm of a Vector/VectorRange, or normalizing a Vector/VectorRange to...

myra::make_UserAction
Action< typename UserClass::Number > make_UserAction(const UserClass &u)
Helper function to adapt some user code (encapsulated in a class) into an Action. ...
Definition: UserAction.h:61

myra::ILUSolver
Incompletely factors A ~= L*U, presents approximate solve() method.
Definition: ILUSolver.h:28

MatrixRange.h
Interface class for representing subranges of dense Matrix&#39;s.

VectorRange.h
Interface class for representing subranges of dense Vector&#39;s.

Action.h
Applies the "Action" of a linear operator, b := A*x, used in iterative solution algorithms.

gemm.h
Variety of routines for mixed dense*sparse or dense*sparse matrix multiplies. The dense*dense case is...

myra::MatrixRange::assign
void assign(const CMatrixRange< Number > &that) const
Assigns the values in *this with the values in that.
Definition: MatrixRange.cpp:268

gmres.h
Linear system solution via gmres (for invertible action A)

myra::SparseMatrix::random
static SparseMatrix< Number > random(int I, int J, int N)
Generates a random SparseMatrix with size IxJ and (approximately) N nonzeros.
Definition: SparseMatrix.cpp:493

SparseMatrix.h
General purpose compressed-sparse-column (CSC) container.

myra::Vector::random
static Vector< Number > random(int N)
Generates a random Vector of specified size.
Definition: Vector.cpp:276

myra
Definition: syntax.dox:1

myra::CMatrixRange
Represents a const MatrixRange.
Definition: bothcat.h:22

IdentityAction.h
An Action that is just the identity operator.

myra_stlprint
Definition: stlprint.h:32

UserAction.h
Adapts user code (encapsulated in a class) into an Action.

myra::make_SolveAction
Action< typename ReflectNumber< Solver >::type > make_SolveAction(const Solver &solver, char op='N')
Free function for making SolveAction&#39;s.
Definition: SolveAction.h:67

stlprint.h
Routines for printing the contents of various std::container&#39;s to a std::ostream using operator <<...

Number.h
Various utility functions/classes related to scalar Number types.

gemv.h
Signatures for sparse matrix * dense vector multiplies. All delegate to gemm() under the hood...

myra::MatrixRange
Represents a mutable MatrixRange.
Definition: conjugate.h:26

myra::Vector::zeros
static Vector< Number > zeros(int N)
Generates a zeros Vector of specified size.
Definition: Vector.cpp:280

Matrix.h
General purpose dense matrix container, O(i*j) storage.

PatternRange.h
Range/Iterator types associated with Pattern.

Vector.h
Container for either a column vector or row vector (depends upon the usage context) ...

myra::ReflectPrecision
Reflects Precision trait for a Number, scalar Number types should specialize it.
Definition: Number.h:33

Pattern.h
Container class for a sparse nonzero pattern, used in reordering/symbolic analysis.

inverse.h
Overwrites a LowerMatrix, DiagonalMatrix, or square Matrix with its own inverse. Or, returns it as a copy.

SolveAction.h
Adapts a class with a .solve() method into an Action.

GemmAction.h
An Action for multiplying by a dense Matrix or SparseMatrix using gemm().

frobenius.h
Returns frobenius norm of a SparseMatrix.

ILUSolver.h
Incomplete LU preconditioner. Presents a solve() function, but it&#39;s only approximate.

cumsum.h
Returns cumulative sum of a std::vector<T>. Reversible by diff()

laplacian2.h
Helper routines for reordering/filling 2D structured grids. Used by many unit tests.

SparseMatrixRange.h
Range/Iterator types associated with SparseMatrix.