MyraMath
TrsmAction


Source: tests/iterative/TrsmAction.cpp

1 // ========================================================================= //
2 // This file is part of MyraMath, copyright (c) 2014-2019 by Ryan A Chilton //
3 // and distributed by MyraCore, LLC. See LICENSE.txt for license terms. //
4 // ========================================================================= //
5 
11 // Containers.
12 #include <myramath/utility/Number.h>
13 #include <myramath/dense/Matrix.h>
14 #include <myramath/dense/MatrixRange.h>
15 #include <myramath/sparse/SparseMatrix.h>
16 #include <myramath/sparse/SparseMatrixRange.h>
17 
18 // Algorithms.
19 #include <myramath/utility/random.h>
20 #include <myramath/dense/gemm.h>
21 #include <myramath/dense/tril.h>
22 #include <myramath/dense/transpose.h>
23 #include <myramath/dense/frobenius.h>
24 #include <myramath/dense/inverse.h>
25 #include <myramath/sparse/gemm.h>
26 #include <myramath/sparse/tril.h>
27 #include <myramath/sparse/transpose.h>
28 
29 // Iterative solve/action stuff.
30 #include <myramath/iterative/Action.h>
31 #include <myramath/iterative/TrsmAction.h>
32 
33 // Reporting.
34 #include <tests/myratest.h>
35 
36 using namespace myra;
37 
38 namespace {
39 
40 // Test make_TrsmAction() for a CMatrixRange
41 template<class Number> void test1(int IJ, int K, typename ReflectPrecision<Number>::type tolerance)
42  {
43  myra::out() << typestring<Number>() << std::endl;
44  typedef typename ReflectPrecision<Number>::type Precision;
45  // Form A.
46  auto A = Matrix<Number>::random(IJ,IJ);
47  A = A + 10*Matrix<Number>::identity(IJ);
48  A = tril(A);
49  auto action = make_TrsmAction(A);
50  // Check action.make_Matrix()
51  {
52  Precision error = frobenius(action.make_Matrix()-inverse(A));
53  myra::out() << " |A(action) - A(dense)| = " << error << std::endl;
54  REQUIRE(error < tolerance);
55  }
56  // Check A_action.multiply(X,B)
57  {
58  auto X = Matrix<Number>::random(IJ,K);
59  auto B = Matrix<Number>::zeros(IJ,K);
60  action.multiply(X,B);
61  Precision error = frobenius(inverse(A)*X-B);
62  myra::out() << " |A*X - B| = " << error << std::endl;
63  REQUIRE(error < tolerance);
64  }
65  // Check A_action.multiply(X,B,alpha,beta)
66  {
67  auto X = Matrix<Number>::random(IJ,K);
68  auto B = Matrix<Number>::random(IJ,K);
69  auto alpha = random<Number>();
70  auto beta = random<Number>();
71  auto C = alpha*inverse(A)*X + beta*B;
72  action.multiply(X,B,alpha,beta);
73  Precision error = frobenius(B-C);
74  myra::out() << " | (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = " << error << std::endl;
75  REQUIRE(error < tolerance);
76  }
77  }
78 
79 // Test make_TrsmAction() for a CSparseMatrixRange
80 template<class Number> void test2(int IJ, int N, int K, typename ReflectPrecision<Number>::type tolerance)
81  {
82  myra::out() << typestring<Number>() << std::endl;
83  typedef typename ReflectPrecision<Number>::type Precision;
84  // Form A.
85  auto A = SparseMatrix<Number>::random(IJ,IJ,N);
86  A = A + SparseMatrix<Number>::identity(IJ);
87  A = tril(A);
88  auto action = make_TrsmAction(A);
89  // Check action.make_Matrix()
90  {
91  Precision error = frobenius(action.make_Matrix()-inverse(A.make_Matrix()));
92  myra::out() << " |A(action) - A(dense)| = " << error << std::endl;
93  REQUIRE(error < tolerance);
94  }
95  // Check A_action.multiply(X,B)
96  {
97  auto X = Matrix<Number>::random(IJ,K);
98  auto B = Matrix<Number>::zeros(IJ,K);
99  action.multiply(X,B);
100  Precision error = frobenius(inverse(A.make_Matrix())*X-B);
101  myra::out() << " |A*X - B| = " << error << std::endl;
102  REQUIRE(error < tolerance);
103  }
104  // Check A_action.multiply(X,B,alpha,beta)
105  {
106  auto X = Matrix<Number>::random(IJ,K);
107  auto B = Matrix<Number>::random(IJ,K);
108  auto alpha = random<Number>();
109  auto beta = random<Number>();
110  auto C = alpha*inverse(A.make_Matrix())*X + beta*B;
111  action.multiply(X,B,alpha,beta);
112  Precision error = frobenius(B-C);
113  myra::out() << " | (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = " << error << std::endl;
114  REQUIRE(error < tolerance);
115  }
116  }
117 
118 } // namespace
119 
120 ADD_TEST("TrsmAction","[iterative]")
121  {
122  // Test dense make_TrsmAction()
123  int IJ = 10;
124  int K = 7;
125  test1<NumberS>(IJ,K,1.0e-2f);
126  test1<NumberD>(IJ,K,1.0e-10);
127  test1<NumberC>(IJ,K,1.0e-2f);
128  test1<NumberZ>(IJ,K,1.0e-10);
129  // Test sparse make_TrsmAction()
130  int N = 20;
131  test2<NumberS>(IJ,N,K,1.0e-2f);
132  test2<NumberD>(IJ,N,K,1.0e-10);
133  test2<NumberC>(IJ,N,K,1.0e-2f);
134  test2<NumberZ>(IJ,N,K,1.0e-10);
135  }
transpose.h
Returns a transposed copy of a SparseMatrix.
MatrixRange.h
Interface class for representing subranges of dense Matrix&#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::Matrix::zeros
static Matrix< Number > zeros(int I, int J)
Generates a zeros Matrix of specified size.
Definition: Matrix.cpp:357
frobenius.h
Routines for computing Frobenius norms of various algebraic containers.
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().
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
Definition: syntax.dox:1
transpose.h
Returns a transposed copy of a Matrix. The inplace version only works on a square operand...
tril.h
Returns the lower triangle of a dense Matrix.
myra::SparseMatrix::identity
static SparseMatrix< Number > identity(int IJ)
Generates an identity SparseMatrix of specified size.
Definition: SparseMatrix.cpp:481
Number.h
Various utility functions/classes related to scalar Number types.
myra::Matrix::identity
static Matrix< Number > identity(int IJ)
Generates an identity Matrix of specified size.
Definition: Matrix.cpp:349
Matrix.h
General purpose dense matrix container, O(i*j) storage.
tril.h
Returns the lower triangle of a SparseMatrix.
myra::ReflectPrecision
Reflects Precision trait for a Number, scalar Number types should specialize it.
Definition: Number.h:33
inverse.h
Overwrites a LowerMatrix, DiagonalMatrix, or square Matrix with its own inverse. Or, returns it as a copy.
random.h
Simplistic random number functions.
gemm.h
Variety of routines all for dense Matrix*Matrix multiplies. Delegates to the BLAS.
SparseMatrixRange.h
Range/Iterator types associated with SparseMatrix.


Results: [PASS]

float
|A(action) - A(dense)| = 6.70731e-10
|A*X - B| = 9.31323e-09
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = 5.96337e-08
double
|A(action) - A(dense)| = 2.42435e-19
|A*X - B| = 0
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = 2.2618e-17
std::complex<float>
|A(action) - A(dense)| = 3.71502e-09
|A*X - B| = 2.69249e-08
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = 1.28792e-07
std::complex<double>
|A(action) - A(dense)| = 8.39443e-18
|A*X - B| = 4.93177e-17
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = 3.86902e-16
float
|A(action) - A(dense)| = 0
|A*X - B| = 2.25003e-07
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = 4.09888e-08
double
|A(action) - A(dense)| = 1.38778e-17
|A*X - B| = 1.89395e-16
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = 3.37447e-16
std::complex<float>
|A(action) - A(dense)| = 1.35961e-06
|A*X - B| = 3.51987e-06
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = 1.25867e-06
std::complex<double>
|A(action) - A(dense)| = 3.39284e-16
|A*X - B| = 1.33362e-15
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = 1.23501e-15


Go back to Summary of /test programs.

Generated on Sun Aug 4 2024 12:26:20 for MyraMath by   doxygen 1.8.13