MyraMath
TrmmAction


Source: tests/iterative/TrmmAction.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/TrmmAction.h>
32 
33 // Reporting.
34 #include <tests/myratest.h>
35 
36 using namespace myra;
37 
38 namespace {
39 
40 // Test make_TrmmAction() 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 + Matrix<Number>::identity(IJ);
48  A = tril(A);
49  auto action = make_TrmmAction(A);
50  // Check action.make_Matrix()
51  {
52  Precision error = frobenius(action.make_Matrix()-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(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*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_TrmmAction() 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_TrmmAction(A);
89  // Check action.make_Matrix()
90  {
91  Precision error = frobenius(action.make_Matrix()-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(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*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("TrmmAction","[iterative]")
121  {
122  // Test dense make_TrmmAction()
123  int IJ = 10;
124  int K = 7;
125  test1<NumberS>(IJ,K,1.0e-4f);
126  test1<NumberD>(IJ,K,1.0e-10);
127  test1<NumberC>(IJ,K,1.0e-4f);
128  test1<NumberZ>(IJ,K,1.0e-10);
129  // Test sparse make_TrmmAction()
130  int N = 20;
131  test2<NumberS>(IJ,N,K,1.0e-4f);
132  test2<NumberD>(IJ,N,K,1.0e-12);
133  test2<NumberC>(IJ,N,K,1.0e-4f);
134  test2<NumberZ>(IJ,N,K,1.0e-12);
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
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
TrmmAction.h
An Action for multiplying by a triangular Matrix, LowerMatrix or SparseMatrix using trmm()...
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)| = 0
|A*X - B| = 4.57091e-07
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = 3.49463e-07
double
|A(action) - A(dense)| = 0
|A*X - B| = 6.51702e-16
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = 7.10873e-16
std::complex<float>
|A(action) - A(dense)| = 0
|A*X - B| = 1.05538e-06
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = 1.06126e-06
std::complex<double>
|A(action) - A(dense)| = 0
|A*X - B| = 2.25379e-15
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = 2.31106e-15
float
|A(action) - A(dense)| = 0
|A*X - B| = 2.69872e-07
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = 1.60144e-07
double
|A(action) - A(dense)| = 0
|A*X - B| = 5.55382e-17
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = 2.58887e-16
std::complex<float>
|A(action) - A(dense)| = 0
|A*X - B| = 3.8289e-07
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = 4.26265e-07
std::complex<double>
|A(action) - A(dense)| = 0
|A*X - B| = 8.45065e-16
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = 1.04646e-15


Go back to Summary of /test programs.

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