MyraMath
hetrf_Matrix


Source: tests/dense/hetrf.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/hetrf.h>
29 #include <myramath/dense/frobenius.h>
30 #include <myramath/dense/hermitian.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 hetrf_inplace() on 'L'ower and 'U'pper MatrixRange's
41 template<class Precision> void test1(int N, Precision tolerance)
42  {
43  typedef std::complex<Precision> Number;
44  myra::out() << typestring<Number>() << std::endl;
45  // Construct random hermitian A.
46  auto A = Matrix<Number>::random(N,N);
47  A = A + hermitian(A);
48 
49  // Check hetrf('L','R'), it factors A = P'*L*R*Q'*I*Q*R'*L'*P
50  {
51  Matrix<Number> L = A;
52  auto result = hetrf_inplace('L','R',L);
53  L = tril(L);
54  Matrix<Number> P = swaps2Matrix<Number>(result.P_swaps);
55  Matrix<Number> Q = swaps2Matrix<Number>(result.Q_swaps);
56  Matrix<Number> R = result.R.make_Matrix();
57  auto I = DiagonalMatrix<Number>::inertia(result.n_plus, result.n_minus);
58  Precision error = frobenius( hermitian(P)*L*R*hermitian(Q)*I*Q*hermitian(R)*hermitian(L)*P-A ) / frobenius(A);
59  myra::out() << " |P'*L*R'*Q'*I*Q*R*L'*P - A| = " << error << std::endl;
60  REQUIRE(error < tolerance);
61  }
62 
63  // Check hetrf('L','L'), it factors A = P'*L'*R'*Q'*I*Q*R*L*P
64  {
65  Matrix<Number> L = A;
66  auto result = hetrf_inplace('L','L',L);
67  L = tril(L);
68  Matrix<Number> P = swaps2Matrix<Number>(result.P_swaps);
69  Matrix<Number> Q = swaps2Matrix<Number>(result.Q_swaps);
70  Matrix<Number> R = result.R.make_Matrix();
71  auto I = DiagonalMatrix<Number>::inertia(result.n_plus, result.n_minus);
72  Precision error = frobenius( hermitian(P)*hermitian(L)*hermitian(R)*hermitian(Q)*I*Q*R*L*P-A ) / frobenius(A);
73  myra::out() << " |P'*L'*R'*Q'*I*Q*R*L*P - A| = " << error << std::endl;
74  REQUIRE(error < tolerance);
75  }
76 
77  // Check hetrf('U','R'), it factors A = P'*U*R*Q'*I*Q*R'*U'*P
78  {
79  Matrix<Number> U = A;
80  auto result = hetrf_inplace('U','R',U);
81  U = triu(U);
82  Matrix<Number> P = swaps2Matrix<Number>(result.P_swaps);
83  Matrix<Number> Q = swaps2Matrix<Number>(result.Q_swaps);
84  Matrix<Number> R = result.R.make_Matrix();
85  auto I = DiagonalMatrix<Number>::inertia(result.n_plus, result.n_minus);
86  Precision error = frobenius( hermitian(P)*U*R*hermitian(Q)*I*Q*hermitian(R)*hermitian(U)*P-A ) / frobenius(A);
87  myra::out() << " |P'*U*R*Q'*I*Q*R'*U'*P - A| = " << error << std::endl;
88  REQUIRE(error < tolerance);
89  }
90 
91  // Check hetrf('U','L'), it factors A = P'*U'*R'*Q'*I*Q*R*U*P
92  {
93  Matrix<Number> U = A;
94  auto result = hetrf_inplace('U','L',U);
95  U = triu(U);
96  Matrix<Number> P = swaps2Matrix<Number>(result.P_swaps);
97  Matrix<Number> Q = swaps2Matrix<Number>(result.Q_swaps);
98  Matrix<Number> R = result.R.make_Matrix();
99  auto I = DiagonalMatrix<Number>::inertia(result.n_plus, result.n_minus);
100  Precision error = frobenius( hermitian(P)*hermitian(U)*hermitian(R)*hermitian(Q)*I*Q*R*U*P-A ) / frobenius(A);
101  myra::out() << " |P'*U'*R'*Q'*I*Q*R*U*P - A| = " << error << std::endl;
102  REQUIRE(error < tolerance);
103  }
104 
105  }
106 
107 // Tests hetrf_inplace() on LowerMatrix's
108 template<class Precision> void test2(int N, Precision tolerance)
109  {
110  typedef std::complex<Precision> Number;
111  myra::out() << typestring<Number>() << std::endl;
112  // Construct random hermitian A.
113  auto A = Matrix<Number>::random(N,N);
114  A = A + hermitian(A);
115  // Extract hetrf(LowerMatrix), it factors A = P'*L*R*Q'*I*Q*R'*L'*P
116  LowerMatrix<Number> L(A);
117  auto result = hetrf_inplace(L);
118  Matrix<Number> P = swaps2Matrix<Number>(result.P_swaps);
119  Matrix<Number> Q = swaps2Matrix<Number>(result.Q_swaps);
120  Matrix<Number> R = result.R.make_Matrix();
121  auto I = DiagonalMatrix<Number>::inertia(result.n_plus, result.n_minus);
122  Matrix<Number> Lf = L.make_Matrix();
123  Precision error = frobenius( hermitian(P)*Lf*R*hermitian(Q)*I*Q*hermitian(R)*hermitian(Lf)*P-A ) / frobenius(A);
124  myra::out() << " |P'*L*R*Q'*I*Q*R'*L'*P-A| = " << error << std::endl;
125  REQUIRE(error < tolerance);
126  }
127 
128 } // namespace
129 
130 ADD_TEST("hetrf_Matrix","[dense][lapack]")
131  {
132  myra::out() << "Checking hetrf(MatrixRange)" << std::endl;
133  test1<float > (169, 1.0e-4f);
134  test1<double> (169, 1.0e-10);
135  }
136 
137 ADD_TEST("hetrf_LowerMatrix","[dense][lapack]")
138  {
139  myra::out() << "Checking hetrf(LowerMatrix)" << std::endl;
140  test2<float > (192, 1.0e-3f);
141  test2<double> (192, 1.0e-8);
142  }
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
hetrf.h
LDL&#39; decompositions for real hermitian Matrix A (indefinite is fine).
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
LowerMatrix.h
Specialized container for a lower triangular matrix, O(N^2/2) storage. Used by symmetry exploiting ma...
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::DiagonalMatrix::inertia
static DiagonalMatrix< Number > inertia(int N_plus, int N_minus)
Generates an inertia DiagonalMatrix, [+I -I].
Definition: DiagonalMatrix.cpp:221
hermitian.h
Returns a hermitian copy of a Matrix. The inplace version only works on a square operand.
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 hetrf(MatrixRange)
std::complex<float>
|P'*L*R'*Q'*I*Q*R*L'*P - A| = 3.1803e-06
|P'*L'*R'*Q'*I*Q*R*L*P - A| = 3.29095e-06
|P'*U*R*Q'*I*Q*R'*U'*P - A| = 3.80268e-06
|P'*U'*R'*Q'*I*Q*R*U*P - A| = 3.62216e-06
std::complex<double>
|P'*L*R'*Q'*I*Q*R*L'*P - A| = 5.94817e-15
|P'*L'*R'*Q'*I*Q*R*L*P - A| = 6.05741e-15
|P'*U*R*Q'*I*Q*R'*U'*P - A| = 5.69709e-15
|P'*U'*R'*Q'*I*Q*R*U*P - A| = 5.42728e-15


Go back to Summary of /test programs.

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