doc/html/test_compose_action1.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 for A,P and and M
 #include <myramath/utility/Number.h>
 #include <myramath/dense/Matrix.h>
 #include <myramath/dense/MatrixRange.h>
 #include <myramath/dense/LowerMatrix.h>
 #include <myramath/dense/LowerMatrixRange.h>
 #include <myramath/dense/syrk.h>
 #include <myramath/dense/trsm.h>
 #include <myramath/dense/potrf.h>

 // Container for x/b
 #include <myramath/dense/Vector.h>
 #include <myramath/dense/VectorRange.h>

 // For examining eigenvalues of L\(A+Q*Q')/L'
 #include <myramath/dense/syev.h>
 #include <myramath/dense/DiagonalMatrix.h>
 #include <myramath/dense/DiagonalMatrixRange.h>
 #include <myramath/dense/euclidean.h>

 // Iterative solve/action stuff.
 #include <myramath/iterative/pcg.h>
 #include <myramath/iterative/gmres.h>
 #include <myramath/iterative/Action.h>
 #include <myramath/iterative/SumAction.h>
 #include <myramath/iterative/ProductAction.h>
 #include <myramath/iterative/GemmAction.h>
 #include <myramath/iterative/SymmAction.h>
 #include <myramath/iterative/TrsmAction.h>
 #include <myramath/iterative/SolveAction.h>
 #include <myramath/iterative/IdentityAction.h>

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

 using namespace myra;
 using namespace myra_stlprint;

 namespace {

 // Tests Action composition using pcg(). A little awkward, because pcg() only handles symmetric Action's
 template<class Precision> void test(int N, int K, Precision tolerance)
   {
   myra::out() << typestring<Precision>() << std::endl;
   // Form symmetric positive definite LowerMatrix A.
   auto G = Matrix<Precision>::random(N,N);
   auto A = syrk(G);
   // Compute cholesky triangle L of A.
   auto L = A;
   potrf_inplace(L);
   // Form forward action (A+PP', a rank-K perturbation of A), then precondition it from both sides.
   auto P = Matrix<Precision>::random(N,K);
   auto action1 = make_SymmAction(A) + make_GemmAction(P)*make_GemmAction(P,'T');
   auto action2 = make_TrsmAction(L) * action1 * make_TrsmAction(L,'T');
   // Examine eigenvalues of action2, should be mostly 1's with K non-1's
   auto XD = syev(action2.make_Matrix());
   myra::out() << "lambda = " << XD.second << std::endl;
   // Make random b and zero x.
   auto b = Vector<Precision>::random(N);
   auto x = Vector<Precision>::zeros(N);
   // Solve action1*x=b, by using pcg(identity,action2) and a bit of trsm() pre/post processing.
   auto b2 = b;
   auto eye = make_IdentityAction<Precision>(N);
   trsm_inplace('L','N',L,b2.column());
   // auto result = pcg(eye,action2,b2,x,tolerance); // doesn't converge as expected .. too many iterations?
   auto result = gmres(eye,action2,b2,x,tolerance); // converges as expected
   trsm_inplace('L','T',L,x.column());
   // Check residual is small.
   Precision residual = euclidean(action1*x-b) / euclidean(b);
   myra::out() << "residual = " << residual << std::endl;
   REQUIRE(residual < tolerance*10.0);
   // Since action2 is a rank-K perturbation of identity, pcg() should converge in O(K) iterations.
   myra::out() << "history = " << result.history << std::endl;
   REQUIRE(result.history.size() <= K+5);
   }

 } // namespace

 ADD_TEST("compose_action1","[iterative]")
   {
   // Test pcg() on a symmetric composition.
   // test<NumberS>(100,5,1.0e-5);
   test<NumberD>(100,10,1.0e-10);
   }

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

potrf.h
Cholesky factorization routines for positive definite matrices.

DiagonalMatrixRange.h
Interface class for representing subranges of DiagonalMatrix&#39;s.

syev.h
Computes all eigenpairs of real symmetric Matrix.

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.

myra::Matrix::random
static Matrix< Number > random(int I, int J)
Generates a random Matrix of specified size.
Definition: Matrix.cpp:353

TrsmAction.h
An Action for solving by a triangular Matrix, LowerMatrix or SparseMatrix using trsm().

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

LowerMatrixRange.h
Range construct for a lower triangular matrix stored in rectangular packed format.

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

SymmAction.h
An Action for multiplying by a symmetric dense Matrix. LowerMatrix or SparseMatrix using symm() ...

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

myra_stlprint
Definition: stlprint.h:32

trsm.h
Routines for backsolving by a triangular Matrix or LowerMatrix.

LowerMatrix.h
Specialized container for a lower triangular matrix, O(N^2/2) storage. Used by symmetry exploiting ma...

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.

SumAction.h
Composes two Action&#39;s A and B, yielding an Action that applies (A+B)*X.

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

ProductAction.h
Composes two Action&#39;s A and B, yielding an Action that cascades A*(B*X)

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

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

myra::Vector::column
CMatrixRange< Number > column() const
Returns an Nx1 column-shaped MatrixRange over all of *this.
Definition: Vector.cpp:152

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

DiagonalMatrix.h
Container for a diagonal matrix, O(n) storage. Used by SVD, row/column scaling, etc.

syrk.h
Routines for symmetric rank-k updates, a specialized form of Matrix*Matrix multiplication.