MyraMath
ctrsmtrmm


Source: tests/dense/tr_roundtrip.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/dense/LowerMatrix.h>
16 #include <myramath/dense/LowerMatrixRange.h>
17 
18 // Data generators.
19 #include <myramath/utility/random.h>
20 #include <myramath/dense/potrf.h>
21 #include <myramath/dense/tril.h>
22 
23 // Algorithms.
24 #include <myramath/dense/gemm.h>
25 #include <myramath/dense/trmm.h>
26 #include <myramath/dense/trsm.h>
27 #include <myramath/dense/transpose.h>
28 #include <myramath/dense/hermitian.h>
29 #include <myramath/dense/conjugate.h>
30 #include <myramath/dense/frobenius.h>
31 
32 // Reporting.
33 #include <tests/myratest.h>
34 
35 using namespace myra;
36 
37 namespace {
38 
39 template<class Number> void test(typename ReflectPrecision<Number>::type tolerance)
40  {
41  typedef typename ReflectPrecision<Number>::type Precision;
42  myra::out() << typestring<Number>() << std::endl;
43  int I = 43;
44  int J = 26;
45  // Compute the cholesky triangle of a random SPD matrix to get a suitably conditioned L.
46  auto G = Matrix<Number>::random(I,I);
47  Matrix<Number> A = G * hermitian(G);
48  potrf_inplace('L','R',A);
49  A = tril(A);
50  LowerMatrix<Number> L(A);
51  Number alpha = random<Number>();
52  char diag = 'N';
53  Number one(1);
54  // Case 1: L\L*X
55  {
56  auto X = Matrix<Number>::random(I,J);
57  Matrix<Number> X2(X);
58  trmm_inplace('L','N',L,X2,diag,alpha);
59  trsm_inplace('L','N',L,X2,diag,one/alpha);
60  Precision error = frobenius(X-X2)/frobenius(X);
61  myra::out() << " |L\\L*X-X| = " << error << std::endl;
62  REQUIRE(error < tolerance);
63  }
64  // Case 2: tranpose(L)\transpose(L)*X
65  {
66  auto X = Matrix<Number>::random(I,J);
67  Matrix<Number> X2(X);
68  trmm_inplace('L','T',L,X2,diag,alpha);
69  trsm_inplace('L','T',L,X2,diag,one/alpha);
70  Precision error = frobenius(X-X2)/frobenius(X);
71  myra::out() << " |transpose(L)\\transpose(L)*X-X| = " << error << std::endl;
72  REQUIRE(error < tolerance);
73  }
74  // Case 3: hermitian(L)\hermitian(L)*X
75  {
76  auto X = Matrix<Number>::random(I,J);
77  Matrix<Number> X2(X);
78  trmm_inplace('L','H',L,X2,diag,alpha);
79  trsm_inplace('L','H',L,X2,diag,one/alpha);
80  Precision error = frobenius(X-X2)/frobenius(X);
81  myra::out() << " |hermitian(L)\\hermitian(L)*X-X| = " << error << std::endl;
82  REQUIRE(error < tolerance);
83  }
84  // Case 4: conjugate(L)\conjugate(L)*X
85  {
86  auto X = Matrix<Number>::random(I,J);
87  Matrix<Number> X2(X);
88  trmm_inplace('L','C',L,X2,diag,alpha);
89  trsm_inplace('L','C',L,X2,diag,one/alpha);
90  Precision error = frobenius(X-X2)/frobenius(X);
91  myra::out() << " |conjugate(L)\\conjugate(L)*X-X| = " << error << std::endl;
92  REQUIRE(error < tolerance);
93  }
94  // Case 5: X*L\L
95  {
96  auto X = Matrix<Number>::random(J,I);
97  Matrix<Number> X2(X);
98  trmm_inplace('R','N',L,X2,diag,alpha);
99  trsm_inplace('R','N',L,X2,diag,one/alpha);
100  Precision error = frobenius(X-X2)/frobenius(X);
101  myra::out() << " |X*L/L-X| = " << error << std::endl;
102  REQUIRE(error < tolerance);
103  }
104  // Case 6: X*tranpose(L)\transpose(L)
105  {
106  auto X = Matrix<Number>::random(J,I);
107  Matrix<Number> X2(X);
108  trmm_inplace('R','T',L,X2,diag,alpha);
109  trsm_inplace('R','T',L,X2,diag,one/alpha);
110  Precision error = frobenius(X-X2)/frobenius(X);
111  myra::out() << " |X*transpose(L)/transpose(L)-X| = " << error << std::endl;
112  REQUIRE(error < tolerance);
113  }
114  // Case 7: X*hermitian(L)\hermitian(L)
115  {
116  auto X = Matrix<Number>::random(J,I);
117  Matrix<Number> X2(X);
118  trmm_inplace('R','H',L,X2,diag,alpha);
119  trsm_inplace('R','H',L,X2,diag,one/alpha);
120  Precision error = frobenius(X-X2)/frobenius(X);
121  myra::out() << " |X*hermitian(L)/hermitian(L)-X| = " << error << std::endl;
122  REQUIRE(error < tolerance);
123  }
124  // Case 8: X*conjugate(L)\conjugate(L)
125  {
126  auto X = Matrix<Number>::random(J,I);
127  Matrix<Number> X2(X);
128  trmm_inplace('R','C',L,X2,diag,alpha);
129  trsm_inplace('R','C',L,X2,diag,one/alpha);
130  Precision error = frobenius(X-X2)/frobenius(X);
131  myra::out() << " |X*conjugate(L)/conjugate(L)-X| = " << error << std::endl;
132  REQUIRE(error < tolerance);
133  }
134  }
135 
136 } // namespace
137 
138 ADD_TEST("strsmtrmm","[dense][lapack]")
139  { test<NumberS>(1.0e-4f); }
140 
141 ADD_TEST("dtrsmtrmm","[dense][lapack]")
142  { test<NumberD>(1.0e-10); }
143 
144 ADD_TEST("ctrsmtrmm","[dense][lapack]")
145  { test<NumberC>(1.0e-4f); }
146 
147 ADD_TEST("ztrsmtrmm","[dense][lapack]")
148  { test<NumberZ>(1.0e-10); }
conjugate.h
Returns a conjugated copy of a Matrix or Vector. Or, conjugate one inplace.
potrf.h
Cholesky factorization routines for positive definite matrices.
MatrixRange.h
Interface class for representing subranges of dense Matrix&#39;s.
myra::Matrix
Tabulates an IxJ matrix. Allows random access, has column major layout to be compatible with BLAS/LAP...
Definition: bdsqr.h:20
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
LowerMatrixRange.h
Range construct for a lower triangular matrix stored in rectangular packed format.
myra
Definition: syntax.dox:1
trsm.h
Routines for backsolving by a triangular Matrix or LowerMatrix.
trmm.h
Routines for multiplying 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...
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.
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
hermitian.h
Returns a hermitian copy of a Matrix. The inplace version only works on a square operand.
random.h
Simplistic random number functions.
myra::LowerMatrix
Stores a lower triangular matrix in rectangular packed format.
Definition: conjugate.h:22
gemm.h
Variety of routines all for dense Matrix*Matrix multiplies. Delegates to the BLAS.


Results: [PASS]

std::complex<float>
|L\L*X-X| = 6.84072e-07
|transpose(L)\transpose(L)*X-X| = 2.9641e-07
|hermitian(L)\hermitian(L)*X-X| = 2.86997e-07
|conjugate(L)\conjugate(L)*X-X| = 8.03142e-07
|X*L/L-X| = 2.89673e-07
|X*transpose(L)/transpose(L)-X| = 6.94173e-07
|X*hermitian(L)/hermitian(L)-X| = 7.42902e-07
|X*conjugate(L)/conjugate(L)-X| = 2.85886e-07


Go back to Summary of /test programs.

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