MyraMath
gelqf_tallskinny


Source: tests/dense/lq.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/Vector.h>
14 #include <myramath/dense/Matrix.h>
15 #include <myramath/dense/VectorRange.h>
16 #include <myramath/dense/MatrixRange.h>
17 
18 // Algorithms.
19 #include <myramath/dense/tril.h>
20 #include <myramath/dense/vertcat.h>
21 #include <myramath/dense/horzcat.h>
22 #include <myramath/dense/gemm.h>
23 #include <myramath/dense/gelqf.h>
24 #include <myramath/dense/ormlq.h>
25 #include <myramath/dense/orglq.h>
26 #include <myramath/dense/frobenius.h>
27 
28 // Reporting.
29 #include <tests/myratest.h>
30 
31 using namespace myra;
32 
33 namespace {
34 
35 template<class Number> void test(int I, int J, typename ReflectPrecision<Number>::type tolerance)
36  {
37  typedef typename ReflectPrecision<Number>::type Precision;
38  myra::out() << typestring<Number>() << std::endl;
39  // Make random A.
40  auto A = Matrix<Number>::random(I,J);
41  // Factor A = LQ inplace, then extract L and form Q explicitly.
42  Matrix<Number> LQ = A;
43  Vector<Number> tau = gelqf_inplace(LQ);
44  Matrix<Number> L = tril(LQ);
45  Matrix<Number> LQ_copy = LQ;
46  CMatrixRange<Number> Q = orglq_inplace(LQ,tau);
47  // Check A = L*Q
48  Precision decomposition_error = frobenius(L*Q-A);
49  myra::out() << " |L*Q-A| = " << decomposition_error << std::endl;
50  REQUIRE(decomposition_error < tolerance);
51  // Check Q'Q = I
52  Precision unitary_error = frobenius(gemm(Q,Q,'H')-Matrix<Number>::identity(Q.I));
53  myra::out() << " |Q*Q'-I| = " << unitary_error << std::endl;
54  REQUIRE(unitary_error < tolerance);
55 
56  // Check Q*B (gemm) = Q*B (gemlq).
57  {
58  int L = 4;
59  auto B = Matrix<Number>::random(J,L);
60  Matrix<Number> B_copy = B;
61  CMatrixRange<Number> QB = ormlq_inplace(LQ_copy, tau, B, 'L', 'N');
62  Precision error = frobenius(gemm(Q,B_copy) - QB);
63  myra::out() << " |Q*B (gemm) - Q*B (gemlq) | = " << error << std::endl;
64  REQUIRE(error < tolerance);
65  }
66  // Check Q^T*B (gemm) = Q^T*B (gemlq)
67  {
68  int L = 3;
69  auto B_copy = Matrix<Number>::random(Q.I,L);
70  // Note B has to be vertically padded with zero's to give it Q.J rows.
71  Matrix<Number> B = vertcat(B_copy, Matrix<Number>(Q.J-Q.I,L));
72  CMatrixRange<Number> QtB = ormlq_inplace(LQ_copy, tau, B, 'L', 'T');
73  Precision error = frobenius(gemm(Q,'T',B_copy) - QtB);
74  myra::out() << " |Q^T*B (gemm) - Q^T*B (gemlq) | = " << error << std::endl;
75  REQUIRE(error < tolerance);
76  }
77  // Check Q^H*B (gemm) = Q^H*B (gemlq)
78  {
79  int L = 3;
80  auto B_copy = Matrix<Number>::random(Q.I,L);
81  // Note B has to be vertically padded with zero's to give it Q.J rows.
82  Matrix<Number> B = vertcat(B_copy, Matrix<Number>(Q.J-Q.I,L));
83  CMatrixRange<Number> QhB = ormlq_inplace(LQ_copy, tau, B, 'L', 'H');
84  Precision error = frobenius(gemm(Q,'H',B_copy) - QhB);
85  myra::out() << " |Q^H*B (gemm) - Q^H*B (gemlq) | = " << error << std::endl;
86  REQUIRE(error < tolerance);
87  }
88  // Check Q^C*B (gemm) = Q^C*B (gemlq).
89  {
90  int L = 4;
91  auto B = Matrix<Number>::random(J,L);
92  Matrix<Number> B_copy = B;
93  CMatrixRange<Number> QcB = ormlq_inplace(LQ_copy, tau, B, 'L', 'C');
94  Precision error = frobenius(gemm(Q,'C',B_copy) - QcB);
95  myra::out() << " |Q^C*B (gemm) - Q^C*B (gemlq) | = " << error << std::endl;
96  REQUIRE(error < tolerance);
97  }
98  // Check B*Q (gemm) = B*Q (gemlq).
99  {
100  int L = 7;
101  auto B_copy = Matrix<Number>::random(L,Q.I);
102  // Note B has to be horizontally padded with zero's to give it Q.J columns.
103  Matrix<Number> B = horzcat(B_copy, Matrix<Number>(L,Q.J-Q.I));
104  CMatrixRange<Number> BQ = ormlq_inplace(LQ_copy, tau, B, 'R', 'N');
105  Precision error = frobenius(gemm(B_copy,Q) - BQ);
106  myra::out() << " |B*Q (gemm) - B*Q (gemlq) | = " << error << std::endl;
107  REQUIRE(error < tolerance);
108  }
109  // Check B*Q^T (gemm) = B*Q^T (gemlq).
110  {
111  int L = 8;
112  auto B = Matrix<Number>::random(L,J);
113  Matrix<Number> B_copy = B;
114  CMatrixRange<Number> BQt = ormlq_inplace(LQ_copy, tau, B, 'R', 'T');
115  Precision error = frobenius(gemm(B_copy,Q,'T') - BQt);
116  myra::out() << " |B*Q^T (gemm) - B*Q^T (gemlq) | = " << error << std::endl;
117  REQUIRE(error < tolerance);
118  }
119  // Check B*Q^H (gemm) = B*Q^H (gemlq).
120  {
121  int L = 8;
122  auto B = Matrix<Number>::random(L,J);
123  Matrix<Number> B_copy = B;
124  CMatrixRange<Number> BQh = ormlq_inplace(LQ_copy, tau, B, 'R', 'H');
125  Precision error = frobenius(gemm(B_copy,Q,'H') - BQh);
126  myra::out() << " |B*Q^H (gemm) - B*Q^H (gemlq) | = " << error << std::endl;
127  REQUIRE(error < tolerance);
128  }
129  // Check B*Q^C (gemm) = B*Q^C (gemlq).
130  {
131  int L = 7;
132  auto B_copy = Matrix<Number>::random(L,Q.I);
133  // Note B has to be horizontally padded with zero's to give it Q.J columns.
134  Matrix<Number> B = horzcat(B_copy, Matrix<Number>(L,Q.J-Q.I));
135  CMatrixRange<Number> BQc = ormlq_inplace(LQ_copy, tau, B, 'R', 'C');
136  Precision error = frobenius(gemm(B_copy,Q,'C') - BQc);
137  myra::out() << " |B*Q^C (gemm) - B*Q^C (gemlq) | = " << error << std::endl;
138  REQUIRE(error < tolerance);
139  }
140  }
141 
142 } // namespace
143 
144 ADD_TEST("gelqf_tallskinny","[dense][lapack]")
145  {
146  int N1 = 6;
147  int N2 = 5;
148  myra::out() << "Testing lq decomposition (tall/skinny).." << std::endl;
149  test<NumberS>(N1,N2,1.0e-5f);
150  test<NumberD>(N1,N2,1.0e-10);
151  test<NumberC>(N1,N2,1.0e-5f);
152  test<NumberZ>(N1,N2,1.0e-10);
153  myra::out() << std::endl;
154  }
155 
156 ADD_TEST("gelqf_shortfat","[dense][lapack]")
157  {
158  int N1 = 6;
159  int N2 = 5;
160  myra::out() << "Testing lq decomposition (short/fat).." << std::endl;
161  test<NumberS>(N2,N1,1.0e-5f);
162  test<NumberD>(N2,N1,1.0e-10);
163  test<NumberC>(N2,N1,1.0e-5f);
164  test<NumberZ>(N2,N1,1.0e-10);
165  myra::out() << std::endl;
166  }
167 
168 ADD_TEST("gelqf_square","[dense][lapack]")
169  {
170  int N = 6;
171  myra::out() << "Testing lq decomposition (square).." << std::endl;
172  test<NumberS>(N,N,1.0e-5f);
173  test<NumberD>(N,N,1.0e-10);
174  test<NumberC>(N,N,1.0e-5f);
175  test<NumberZ>(N,N,1.0e-10);
176  myra::out() << std::endl;
177  }
MatrixRange.h
Interface class for representing subranges of dense Matrix&#39;s.
ormlq.h
Routines to apply a Householder Q, as previously factored via gelqf_inplace()
VectorRange.h
Interface class for representing subranges of dense Vector&#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
horzcat.h
Routines to concatenate Matrix&#39;s in left/right fashion.
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
Definition: syntax.dox:1
orglq.h
Routines to reconstruct a Householder Q, as previously factored via gelqf_inplace() ...
myra::CMatrixRange
Represents a const MatrixRange.
Definition: bothcat.h:22
tril.h
Returns the lower triangle of a dense Matrix.
Number.h
Various utility functions/classes related to scalar Number types.
myra::CMatrixRange::J
int J
---------— Data members, all public ----------------—
Definition: MatrixRange.h:241
myra::CMatrixRange::I
int I
---------— Data members, all public ----------------—
Definition: MatrixRange.h:240
Matrix.h
General purpose dense matrix container, O(i*j) storage.
myra::Vector
Tabulates a vector of length N, allows random access.
Definition: conjugate.h:21
Vector.h
Container for either a column vector or row vector (depends upon the usage context) ...
myra::ReflectPrecision
Reflects Precision trait for a Number, scalar Number types should specialize it.
Definition: Number.h:33
gemm.h
Variety of routines all for dense Matrix*Matrix multiplies. Delegates to the BLAS.
vertcat.h
Routines to concatenate Matrix&#39;s in top/bottom fashion.
gelqf.h
Routines to compute a Householder LQ decomposition.


Results: [PASS]

Testing lq decomposition (tall/skinny)..
float
|L*Q-A| = 6.73072e-07
|Q*Q'-I| = 7.11805e-07
|Q*B (gemm) - Q*B (gemlq) | = 5.94006e-07
|Q^T*B (gemm) - Q^T*B (gemlq) | = 4.07677e-07
|Q^H*B (gemm) - Q^H*B (gemlq) | = 4.11075e-07
|Q^C*B (gemm) - Q^C*B (gemlq) | = 3.40314e-07
|B*Q (gemm) - B*Q (gemlq) | = 7.27474e-07
|B*Q^T (gemm) - B*Q^T (gemlq) | = 5.1378e-07
|B*Q^H (gemm) - B*Q^H (gemlq) | = 5.45857e-07
|B*Q^C (gemm) - B*Q^C (gemlq) | = 5.85929e-07
double
|L*Q-A| = 7.42332e-16
|Q*Q'-I| = 5.68639e-16
|Q*B (gemm) - Q*B (gemlq) | = 7.90098e-16
|Q^T*B (gemm) - Q^T*B (gemlq) | = 6.61574e-16
|Q^H*B (gemm) - Q^H*B (gemlq) | = 5.5706e-16
|Q^C*B (gemm) - Q^C*B (gemlq) | = 7.21822e-16
|B*Q (gemm) - B*Q (gemlq) | = 9.68257e-16
|B*Q^T (gemm) - B*Q^T (gemlq) | = 1.17306e-15
|B*Q^H (gemm) - B*Q^H (gemlq) | = 1.11006e-15
|B*Q^C (gemm) - B*Q^C (gemlq) | = 1.09735e-15
std::complex<float>
|L*Q-A| = 5.96279e-07
|Q*Q'-I| = 3.87752e-07
|Q*B (gemm) - Q*B (gemlq) | = 6.28085e-07
|Q^T*B (gemm) - Q^T*B (gemlq) | = 6.99016e-07
|Q^H*B (gemm) - Q^H*B (gemlq) | = 4.70332e-07
|Q^C*B (gemm) - Q^C*B (gemlq) | = 6.19477e-07
|B*Q (gemm) - B*Q (gemlq) | = 1.11079e-06
|B*Q^T (gemm) - B*Q^T (gemlq) | = 9.44014e-07
|B*Q^H (gemm) - B*Q^H (gemlq) | = 9.70157e-07
|B*Q^C (gemm) - B*Q^C (gemlq) | = 1.12128e-06
std::complex<double>
|L*Q-A| = 9.86202e-16
|Q*Q'-I| = 5.53216e-16
|Q*B (gemm) - Q*B (gemlq) | = 1.7541e-15
|Q^T*B (gemm) - Q^T*B (gemlq) | = 1.00391e-15
|Q^H*B (gemm) - Q^H*B (gemlq) | = 1.05188e-15
|Q^C*B (gemm) - Q^C*B (gemlq) | = 1.31897e-15
|B*Q (gemm) - B*Q (gemlq) | = 2.08321e-15
|B*Q^T (gemm) - B*Q^T (gemlq) | = 1.85926e-15
|B*Q^H (gemm) - B*Q^H (gemlq) | = 1.43506e-15
|B*Q^C (gemm) - B*Q^C (gemlq) | = 1.6724e-15


Go back to Summary of /test programs.

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