MyraMath
sytrf_Matrix


Source: tests/dense/sytrf.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 #include <myramath/dense/DiagonalMatrix.h>
18 #include <myramath/dense/DiagonalMatrixRange.h>
19 #include <myramath/dense/intRange.h>
20 
21 // Data generators.
22 #include <myramath/dense/triu.h>
23 #include <myramath/dense/tril.h>
24 
25 // Algorithms.
26 #include <myramath/dense/gemm.h>
27 #include <myramath/dense/dimm.h>
28 #include <myramath/dense/sytrf.h>
29 #include <myramath/dense/frobenius.h>
30 #include <myramath/dense/transpose.h>
31 #include <myramath/dense/swaps.h>
32 
33 // Reporting.
34 #include <tests/myratest.h>
35 
36 using namespace myra;
37 
38 namespace {
39 
40 // Tests sytrf_inplace() on 'L'ower and 'U'pper MatrixRange's
41 template<class Number> void test1(int N, typename ReflectPrecision<Number>::type tolerance)
42  {
43  typedef typename ReflectPrecision<Number>::type Precision;
44  myra::out() << typestring<Number>() << std::endl;
45  // Construct random symmetric A.
46  auto A = Matrix<Number>::random(N,N);
47  A = A + transpose(A);
48  // Check sytrf('L','R'), it factors A = P'*L*R*Q'*I*Q*R'*L'*P
49  {
50  Matrix<Number> L = A;
51  auto result = sytrf_inplace('L','R',L);
52  L = tril(L);
53  Matrix<Number> P = swaps2Matrix<Number>(result.P_swaps);
54  Matrix<Number> Q = swaps2Matrix<Number>(result.Q_swaps);
55  Matrix<Number> R = result.R.make_Matrix();
56  auto I = DiagonalMatrix<Number>::inertia(result.n_plus, result.n_minus);
57  Precision error = frobenius( transpose(P)*L*R*transpose(Q)*I*Q*transpose(R)*transpose(L)*P-A ) / frobenius(A);
58  myra::out() << " |P'*L*R'*Q'*I*Q*R*L'*P - A| = " << error << std::endl;
59  REQUIRE(error < tolerance);
60  }
61  // Check sytrf('L','L'), it factors A = P'*L'*R'*Q'*I*Q*R*L*P
62  {
63  Matrix<Number> L = A;
64  auto result = sytrf_inplace('L','L',L);
65  L = tril(L);
66  Matrix<Number> P = swaps2Matrix<Number>(result.P_swaps);
67  Matrix<Number> Q = swaps2Matrix<Number>(result.Q_swaps);
68  Matrix<Number> R = result.R.make_Matrix();
69  auto I = DiagonalMatrix<Number>::inertia(result.n_plus, result.n_minus);
70  Precision error = frobenius( transpose(P)*transpose(L)*transpose(R)*transpose(Q)*I*Q*R*L*P-A ) / frobenius(A);
71  myra::out() << " |P'*L'*R'*Q'*I*Q*R*L*P - A| = " << error << std::endl;
72  REQUIRE(error < tolerance);
73  }
74  // Check sytrf('U','R'), it factors A = P'*U*R*Q'*I*Q*R'*U'*P
75  {
76  Matrix<Number> U = A;
77  auto result = sytrf_inplace('U','R',U);
78  U = triu(U);
79  Matrix<Number> P = swaps2Matrix<Number>(result.P_swaps);
80  Matrix<Number> Q = swaps2Matrix<Number>(result.Q_swaps);
81  Matrix<Number> R = result.R.make_Matrix();
82  auto I = DiagonalMatrix<Number>::inertia(result.n_plus, result.n_minus);
83  Precision error = frobenius( transpose(P)*U*R*transpose(Q)*I*Q*transpose(R)*transpose(U)*P-A ) / frobenius(A);
84  myra::out() << " |P'*U*R*Q'*I*Q*R'*U'*P - A| = " << error << std::endl;
85  REQUIRE(error < tolerance);
86  }
87  // Check sytrf('U','L'), it factors A = P'*U'*R'*Q'*I*Q*R*U*P
88  {
89  Matrix<Number> U = A;
90  auto result = sytrf_inplace('U','L',U);
91  U = triu(U);
92  Matrix<Number> P = swaps2Matrix<Number>(result.P_swaps);
93  Matrix<Number> Q = swaps2Matrix<Number>(result.Q_swaps);
94  Matrix<Number> R = result.R.make_Matrix();
95  auto I = DiagonalMatrix<Number>::inertia(result.n_plus, result.n_minus);
96  Precision error = frobenius( transpose(P)*transpose(U)*transpose(R)*transpose(Q)*I*Q*R*U*P-A ) / frobenius(A);
97  myra::out() << " |P'*U'*R'*Q'*I*Q*R*U*P - A| = " << error << std::endl;
98  REQUIRE(error < tolerance);
99  }
100  }
101 
102 // Tests sytrf_inplace() on LowerMatrix's
103 template<class Number> void test2(int N, typename ReflectPrecision<Number>::type tolerance)
104  {
105  typedef typename ReflectPrecision<Number>::type Precision;
106  myra::out() << typestring<Number>() << std::endl;
107  // Construct random symmetric A.
108  auto A = Matrix<Number>::random(N,N);
109  A = A + transpose(A);
110  // Extract sytrf(LowerMatrix), it factors A = P'*L*R*Q'*I*Q*R'*L'*P
111  LowerMatrix<Number> L(tril(A));
112  auto result = sytrf_inplace(L);
113  Matrix<Number> P = swaps2Matrix<Number>(result.P_swaps);
114  Matrix<Number> Q = swaps2Matrix<Number>(result.Q_swaps);
115  Matrix<Number> R = result.R.make_Matrix();
116  auto I = DiagonalMatrix<Number>::inertia(result.n_plus, result.n_minus);
117  Matrix<Number> Lf = L.make_Matrix();
118  Precision error = frobenius( transpose(P)*Lf*R*transpose(Q)*I*Q*transpose(R)*transpose(Lf)*P-A ) / frobenius(A);
119  myra::out() << " |P'*L*R*Q'*I*Q*R'*L'*P - A| = " << error << std::endl;
120  REQUIRE(error < tolerance);
121  }
122 
123 } // namespace
124 
125 ADD_TEST("sytrf_Matrix","[dense][lapack]")
126  {
127  myra::out() << "Checking sytrf(MatrixRange)" << std::endl;
128  test1<NumberS>(169,1.0e-4f);
129  test1<NumberD>(169,1.0e-10);
130  test1<NumberC>(169,1.0e-4f);
131  test1<NumberZ>(169,1.0e-10);
132  myra::out() << std::endl;
133  }
134 
135 ADD_TEST("sytrf_LowerMatrix","[dense][lapack]")
136  {
137  myra::out() << "Checking sytrf(LowerMatrix)" << std::endl;
138  test2<NumberS>(192,1.0e-3f);
139  test2<NumberD>(192,1.0e-8);
140  test2<NumberC>(192,1.0e-3f);
141  test2<NumberZ>(192,1.0e-8);
142  myra::out() << std::endl;
143  }
DiagonalMatrixRange.h
Interface class for representing subranges of DiagonalMatrix&#39;s.
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
dimm.h
Routines for multiplying by a DiagonalMatrix.
LowerMatrixRange.h
Range construct for a lower triangular matrix stored in rectangular packed format.
myra
Definition: syntax.dox:1
sytrf.h
LDL&#39; decompositions for real symmetric Matrix A (indefinite is fine).
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.
triu.h
Returns the upper triangle of a dense Matrix.
Number.h
Various utility functions/classes related to scalar Number types.
swaps.h
Routines related to swap sequences, often used during pivoting.
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
myra::DiagonalMatrix::inertia
static DiagonalMatrix< Number > inertia(int N_plus, int N_minus)
Generates an inertia DiagonalMatrix, [+I -I].
Definition: DiagonalMatrix.cpp:221
myra::LowerMatrix
Stores a lower triangular matrix in rectangular packed format.
Definition: conjugate.h:22
DiagonalMatrix.h
Container for a diagonal matrix, O(n) storage. Used by SVD, row/column scaling, etc.
gemm.h
Variety of routines all for dense Matrix*Matrix multiplies. Delegates to the BLAS.
intRange.h
Interface class for representing subranges of contiguous int&#39;s.


Results: [PASS]

Checking sytrf(MatrixRange)
float
|P'*L*R'*Q'*I*Q*R*L'*P - A| = 1.54938e-06
|P'*L'*R'*Q'*I*Q*R*L*P - A| = 1.55685e-06
|P'*U*R*Q'*I*Q*R'*U'*P - A| = 1.65564e-06
|P'*U'*R'*Q'*I*Q*R*U*P - A| = 1.65183e-06
double
|P'*L*R'*Q'*I*Q*R*L'*P - A| = 2.80083e-15
|P'*L'*R'*Q'*I*Q*R*L*P - A| = 2.87291e-15
|P'*U*R*Q'*I*Q*R'*U'*P - A| = 3.17302e-15
|P'*U'*R'*Q'*I*Q*R*U*P - A| = 3.03433e-15
std::complex<float>
|P'*L*R'*Q'*I*Q*R*L'*P - A| = 3.32412e-06
|P'*L'*R'*Q'*I*Q*R*L*P - A| = 3.52411e-06
|P'*U*R*Q'*I*Q*R'*U'*P - A| = 3.42557e-06
|P'*U'*R'*Q'*I*Q*R*U*P - A| = 3.13042e-06
std::complex<double>
|P'*L*R'*Q'*I*Q*R*L'*P - A| = 5.31652e-15
|P'*L'*R'*Q'*I*Q*R*L*P - A| = 5.65931e-15
|P'*U*R*Q'*I*Q*R'*U'*P - A| = 6.94256e-15
|P'*U'*R'*Q'*I*Q*R*U*P - A| = 6.72254e-15


Go back to Summary of /test programs.

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