doc/html/test_pcg.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/Vector.h>
 #include <myramath/dense/VectorRange.h>
 #include <myramath/sparse/SparseMatrix.h>
 #include <myramath/sparse/SparseMatrixRange.h>
 #include <myramath/iterative/ICholeskySolver.h>

 // Algorithms.
 #include <myramath/dense/euclidean.h>
 #include <myramath/sparse/gemv.h>
 #include <myramath/sparse/phasor.h>
 #include <myramath/sparse/laplacian2.h>

 // Iterative solve/action stuff.
 #include <myramath/iterative/Action.h>
 #include <myramath/iterative/GemmAction.h>
 #include <myramath/iterative/SolveAction.h>
 #include <myramath/iterative/IdentityAction.h>
 #include <myramath/iterative/pcg.h>

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

 using namespace myra;
 using namespace myra_stlprint;

 namespace {

 template<class Precision> void test_real(int X, int Y, Precision tolerance)
   {
   myra::out() << typestring<Precision>() << std::endl;
   // Construct 2D laplacian.
   int N = X*Y;
   auto A = laplacian2<Precision>(X,Y);
   // Construct random x/b vectors.
   auto x = Vector<Precision>::random(N);
   auto b = A*x;
   // Solve using pcg() - no preconditioner.
     {
     auto m = make_IdentityAction<Precision>(N);
     auto a = make_GemmAction(A);
     x = Vector<Precision>::zeros(N);
     auto result = pcg(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;
     REQUIRE(error < tolerance);
     }
   // Solve using pcg() - incomplete cholesky preconditioner.
     {
     ICholeskySolver<Precision> M(A);
     auto m = make_SolveAction(M);
     auto a = make_GemmAction(A);
     x = Vector<Precision>::zeros(N);
     auto result = pcg(m,a,b,x,tolerance);
     Precision error = euclidean(A*x-b) / euclidean(b);
     myra::out() << "  IC 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;
     REQUIRE(error < tolerance);
     }
   }

 // Tests complex pcg(), by solving a 2D graph laplacian (both with and without a preconditioner).
 template<class Precision> void test_complex(int X, int Y, Precision tolerance)
   {
   typedef std::complex<Precision> Number;
   myra::out() << typestring<Number>() << std::endl;
   // Construct 2D laplacian.
   int N = X*Y;
   auto A = laplacian2<Number>(X,Y);
   phasor_hermitian(A);
   // Construct random x/b vectors.
   auto x = Vector<Number>::random(N);
   auto b = A*x;
   // Solve using pcg() - no preconditioner.
     {
     auto m = make_IdentityAction<Number>(N);
     auto a = make_GemmAction(A);
     x = Vector<Number>::zeros(N);
     auto result = pcg(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;
     REQUIRE(error < tolerance);
     }
   // Solve using pcg() - incomplete cholesky preconditioner.
     {
     ICholeskySolver<Number> M(A);
     auto m = make_SolveAction(M);
     auto a = make_GemmAction(A);
     x = Vector<Number>::zeros(N);
     auto result = pcg(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;
     REQUIRE(error < tolerance);
     }
   }

 } // namespace

 ADD_TEST("pcg","[iterative]")
   {
   test_real<NumberS>(20,20,1.0e-4f);
   test_real<NumberD>(20,20,1.0e-8);
   test_complex<NumberS>(20,20,1.0e-4f);
   test_complex<NumberD>(20,20,1.0e-8);
   }
euclidean.h
Routines for computing euclidean norm of a Vector/VectorRange, or normalizing a Vector/VectorRange to...

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.

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

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

myra_stlprint
Definition: stlprint.h:32

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::Vector::zeros
static Vector< Number > zeros(int N)
Generates a zeros Vector of specified size.
Definition: Vector.cpp:280

pcg.h
Linear system solution via conjugate gradients (symmetric positive definite action A) ...

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

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().

myra::ICholeskySolver
Incompletely factors A ~= L*L&#39;, presents approximate solve() method.
Definition: ICholeskySolver.h:28

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

phasor.h
Applies random phase shifts to a complex square SparseMatrix.

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

SparseMatrixRange.h
Range/Iterator types associated with SparseMatrix.