MyraMath
GemmAction


Source: tests/iterative/GemmAction.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/sparse/gemm.h>
22 #include <myramath/dense/frobenius.h>
23 
24 // Iterative solve/action stuff.
25 #include <myramath/iterative/Action.h>
26 #include <myramath/iterative/GemmAction.h>
27 
28 // Reporting.
29 #include <tests/myratest.h>
30 #include <myramath/utility/stlprint.h>
31 
32 using namespace myra;
33 using namespace myra_stlprint;
34 
35 namespace {
36 
37 // Test make_GemmAction() for a CMatrixRange
38 template<class Number> void test1(int I, int J, int K, typename ReflectPrecision<Number>::type tolerance)
39  {
40  myra::out() << typestring<Number>() << std::endl;
41  typedef typename ReflectPrecision<Number>::type Precision;
42  // Make a random Matrix A, wrap it up into a GemmAction
43  auto A = Matrix<Number>::random(I,J);
44  auto action = make_GemmAction(A);
45  // Check action.make_Matrix()
46  {
47  Precision error = frobenius(action.make_Matrix()-A);
48  myra::out() << " |A(action) - A(dense)| = " << error << std::endl;
49  REQUIRE(error < tolerance);
50  }
51  // Check A_action.multiply(X,B)
52  {
53  auto X = Matrix<Number>::random(J,K);
54  auto B = Matrix<Number>::zeros(I,K);
55  action.multiply(X,B);
56  Precision error = frobenius(A*X-B);
57  myra::out() << " |A*X - B| = " << error << std::endl;
58  REQUIRE(error < tolerance);
59  }
60  // Check A_action.multiply(X,B,alpha,beta)
61  {
62  auto X = Matrix<Number>::random(J,K);
63  auto B = Matrix<Number>::random(I,K);
64  auto alpha = random<Number>();
65  auto beta = random<Number>();
66  auto C = alpha*A*X + beta*B;
67  action.multiply(X,B,alpha,beta);
68  Precision error = frobenius(B-C);
69  myra::out() << " | (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = " << error << std::endl;
70  REQUIRE(error < tolerance);
71  }
72  }
73 
74 // Test make_GemmAction() for a CSparseMatrixRange
75 template<class Number> void test2(int I, int J, int N, int K, typename ReflectPrecision<Number>::type tolerance)
76  {
77  myra::out() << typestring<Number>() << std::endl;
78  typedef typename ReflectPrecision<Number>::type Precision;
79  // Make a random SparseMatrix A, wrap it up into a GemmAction
80  auto A = SparseMatrix<Number>::random(I,J,N);
81  auto action = make_GemmAction(A);
82  // Check action.make_Matrix()
83  {
84  Precision error = frobenius(action.make_Matrix()-A.make_Matrix());
85  myra::out() << " |A(action) - A(dense)| = " << error << std::endl;
86  REQUIRE(error < tolerance);
87  }
88  // Check A_action.multiply(X,B)
89  {
90  auto X = Matrix<Number>::random(J,K);
91  auto B = Matrix<Number>::zeros(I,K);
92  action.multiply(X,B);
93  Precision error = frobenius(A*X-B);
94  myra::out() << " |A*X - B| = " << error << std::endl;
95  REQUIRE(error < tolerance);
96  }
97  // Check A_action.multiply(X,B,alpha,beta)
98  {
99  auto X = Matrix<Number>::random(J,K);
100  auto B = Matrix<Number>::random(I,K);
101  auto alpha = random<Number>();
102  auto beta = random<Number>();
103  auto C = alpha*A*X + beta*B;
104  action.multiply(X,B,alpha,beta);
105  Precision error = frobenius(B-C);
106  myra::out() << " | (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = " << error << std::endl;
107  REQUIRE(error < tolerance);
108  }
109  }
110 
111 } // namespace
112 
113 ADD_TEST("GemmAction","[iterative]")
114  {
115  // Test dense make_GemmAction()
116  int I = 10;
117  int J = 5;
118  int K = 7;
119  test1<NumberS>(I,J,K,1.0e-4f);
120  test1<NumberD>(I,J,K,1.0e-12);
121  test1<NumberC>(I,J,K,1.0e-4f);
122  test1<NumberZ>(I,J,K,1.0e-12);
123  // Test sparse make_GemmAction()
124  int N = 20;
125  test2<NumberS>(I,J,N,K,1.0e-4f);
126  test2<NumberD>(I,J,N,K,1.0e-12);
127  test2<NumberC>(I,J,N,K,1.0e-4f);
128  test2<NumberZ>(I,J,N,K,1.0e-12);
129  }
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
myra_stlprint
Definition: stlprint.h:32
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.
Matrix.h
General purpose dense matrix container, O(i*j) storage.
myra::ReflectPrecision
Reflects Precision trait for a Number, scalar Number types should specialize it.
Definition: Number.h:33
GemmAction.h
An Action for multiplying by a dense Matrix or SparseMatrix using gemm().
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| = 0
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = 0
double
|A(action) - A(dense)| = 0
|A*X - B| = 0
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = 0
std::complex<float>
|A(action) - A(dense)| = 0
|A*X - B| = 0
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = 7.53541e-07
std::complex<double>
|A(action) - A(dense)| = 0
|A*X - B| = 0
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = 1.4281e-15
float
|A(action) - A(dense)| = 0
|A*X - B| = 0
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = 1.07195e-07
double
|A(action) - A(dense)| = 0
|A*X - B| = 0
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = 3.09384e-16
std::complex<float>
|A(action) - A(dense)| = 0
|A*X - B| = 0
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = 3.19768e-07
std::complex<double>
|A(action) - A(dense)| = 0
|A*X - B| = 0
| (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = 7.59135e-16


Go back to Summary of /test programs.

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