MyraMath
geqrf_tallskinny


Source: tests/dense/qr.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/triu.h>
20 #include <myramath/dense/vertcat.h>
21 #include <myramath/dense/horzcat.h>
22 #include <myramath/dense/gemm.h>
23 #include <myramath/dense/geqrf.h>
24 #include <myramath/dense/ormqr.h>
25 #include <myramath/dense/orgqr.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 
40  // Make random A.
41  auto A = Matrix<Number>::random(I,J);
42  // Factor A = QR inplace, then extract R and form Q explicitly.
43  Matrix<Number> QR = A;
44  Vector<Number> tau = geqrf_inplace(QR);
45  Matrix<Number> R = triu(QR);
46  Matrix<Number> QR_copy = QR;
47  CMatrixRange<Number> Q = orgqr_inplace(QR,tau);
48 
49  // Check A = Q*R
50  Precision decomposition_error = frobenius(Q*R-A);
51  myra::out() << " |Q*R-A| = " << decomposition_error << std::endl;
52  REQUIRE(decomposition_error < tolerance);
53  // Check Q'Q = I
54  Precision unitary_error = frobenius(gemm(Q,'H',Q)-Matrix<Number>::identity(Q.J));
55  myra::out() << " |Q'*Q-I| = " << unitary_error << std::endl;
56  REQUIRE(unitary_error < tolerance);
57  // Check Q*B (gemm) = Q*B (gemqr).
58  {
59  int L = 4;
60  auto B_copy = Matrix<Number>::random(Q.J,L);
61  // Note B has to be vertically padded with zero's to give it Q.I rows.
62  Matrix<Number> B = vertcat(B_copy, Matrix<Number>(Q.I-Q.J,L));
63  CMatrixRange<Number> QB = ormqr_inplace(QR_copy, tau, B, 'L', 'N');
64  Precision error = frobenius(gemm(Q,B_copy) - QB);
65  myra::out() << " |Q*B (gemm) - Q*B (gemqr) | = " << error << std::endl;
66  REQUIRE(error < tolerance);
67  }
68  // Check Q^T*B (gemm) = Q^T*B (gemqr)
69  {
70  int L = 3;
71  auto B = Matrix<Number>::random(I,L);
72  Matrix<Number> B_copy = B;
73  CMatrixRange<Number> QtB = ormqr_inplace(QR_copy, tau, B, 'L', 'T');
74  Precision error = frobenius(gemm(Q,'T',B_copy) - QtB);
75  myra::out() << " |Q^T*B (gemm) - Q^T*B (gemqr) | = " << error << std::endl;
76  REQUIRE(error < tolerance);
77  }
78  // Check Q^H*B (gemm) = Q^H*B (gemqr)
79  {
80  int L = 3;
81  auto B = Matrix<Number>::random(I,L);
82  Matrix<Number> B_copy = B;
83  CMatrixRange<Number> QhB = ormqr_inplace(QR_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 (gemqr) | = " << error << std::endl;
86  REQUIRE(error < tolerance);
87  }
88  // Check Q^C*B (gemm) = Q^C*B (gemqr).
89  {
90  int L = 4;
91  auto B_copy = Matrix<Number>::random(Q.J,L);
92  // Note B has to be vertically padded with zero's to give it Q.I rows.
93  Matrix<Number> B = vertcat(B_copy, Matrix<Number>(Q.I-Q.J,L));
94  CMatrixRange<Number> QcB = ormqr_inplace(QR_copy, tau, B, 'L', 'C');
95  Precision error = frobenius(gemm(Q,'C',B_copy) - QcB);
96  myra::out() << " |Q^C*B (gemm) - Q^C*B (gemqr) | = " << error << std::endl;
97  REQUIRE(error < tolerance);
98  }
99 
100  // Check B*Q (gemm) = B*Q (gemqr).
101  {
102  int L = 7;
103  auto B = Matrix<Number>::random(L,I);
104  Matrix<Number> B_copy = B;
105  CMatrixRange<Number> BQ = ormqr_inplace(QR_copy, tau, B, 'R', 'N');
106  Precision error = frobenius(gemm(B_copy,Q) - BQ);
107  myra::out() << " |B*Q (gemm) - B*Q (gemqr) | = " << error << std::endl;
108  REQUIRE(error < tolerance);
109  }
110  // Check B*Q^T (gemm) = B*Q^T (gemqr).
111  {
112  int L = 8;
113  auto B_copy = Matrix<Number>::random(L,Q.J);
114  // Note B has be horizontally padded with zeros to give it Q.I columns.
115  Matrix<Number> B = horzcat(B_copy, Matrix<Number>(L,Q.I-Q.J));
116  CMatrixRange<Number> BQt = ormqr_inplace(QR_copy, tau, B, 'R', 'T');
117  Precision error = frobenius(gemm(B_copy,Q,'T') - BQt);
118  myra::out() << " |B*Q^T (gemm) - B*Q^T (gemqr) | = " << error << std::endl;
119  REQUIRE(error < tolerance);
120  }
121  // Check B*Q^H (gemm) = B*Q^H (gemqr).
122  {
123  int L = 8;
124  auto B_copy = Matrix<Number>::random(L,Q.J);
125  // Note B has be horizontally padded with zeros to give it Q.I columns.
126  Matrix<Number> B = horzcat(B_copy, Matrix<Number>(L,Q.I-Q.J));
127  CMatrixRange<Number> BQh = ormqr_inplace(QR_copy, tau, B, 'R', 'H');
128  Precision error = frobenius(gemm(B_copy,Q,'H') - BQh);
129  myra::out() << " |B*Q^H (gemm) - B*Q^H (gemqr) | = " << error << std::endl;
130  REQUIRE(error < tolerance);
131  }
132  // Check B*Q^C (gemm) = B*Q^C (gemqr).
133  {
134  int L = 7;
135  auto B = Matrix<Number>::random(L,I);
136  Matrix<Number> B_copy = B;
137  CMatrixRange<Number> BQc = ormqr_inplace(QR_copy, tau, B, 'R', 'C');
138  Precision error = frobenius(gemm(B_copy,Q,'C') - BQc);
139  myra::out() << " |B*Q^C (gemm) - B*Q^C (gemqr) | = " << error << std::endl;
140  REQUIRE(error < tolerance);
141  }
142  }
143 
144 } // namespace
145 
146 ADD_TEST("geqrf_tallskinny","[dense][lapack]")
147  {
148  int N1 = 6;
149  int N2 = 5;
150  myra::out() << "Testing qr decomposition (tall/skinny).." << std::endl;
151  test<NumberS>(N1,N2,1.0e-5f);
152  test<NumberD>(N1,N2,1.0e-10);
153  test<NumberC>(N1,N2,1.0e-5f);
154  test<NumberZ>(N1,N2,1.0e-10);
155  myra::out() << std::endl;
156  }
157 
158 ADD_TEST("geqrf_shortfat","[dense][lapack]")
159  {
160  int N1 = 6;
161  int N2 = 5;
162  myra::out() << "Testing qr decomposition (short/fat).." << std::endl;
163  test<NumberS>(N2,N1,1.0e-5f);
164  test<NumberD>(N2,N1,1.0e-10);
165  test<NumberC>(N2,N1,1.0e-5f);
166  test<NumberZ>(N2,N1,1.0e-10);
167  myra::out() << std::endl;
168  }
169 
170 ADD_TEST("geqrf_square","[dense][lapack]")
171  {
172  int N = 6;
173  myra::out() << "Testing qr decomposition (square).." << std::endl;
174  test<NumberS>(N,N,1.0e-5f);
175  test<NumberD>(N,N,1.0e-10);
176  test<NumberC>(N,N,1.0e-5f);
177  test<NumberZ>(N,N,1.0e-10);
178  myra::out() << std::endl;
179  }
180 
181 
orgqr.h
Routines to reconstruct a Householder Q, as previously factored via geqrf_inplace() ...
MatrixRange.h
Interface class for representing subranges of dense Matrix&#39;s.
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
geqrf.h
Routines to compute a Householder QR decomposition.
myra
Definition: syntax.dox:1
myra::CMatrixRange
Represents a const MatrixRange.
Definition: bothcat.h:22
triu.h
Returns the upper 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
ormqr.h
Routines to apply a Householder Q, as previously factored via geqrf_inplace()
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.


Results: [PASS]

Testing qr decomposition (tall/skinny)..
float
|Q*R-A| = 2.92761e-07
|Q'*Q-I| = 3.50315e-07
|Q*B (gemm) - Q*B (gemqr) | = 4.20892e-07
|Q^T*B (gemm) - Q^T*B (gemqr) | = 3.43152e-07
|Q^H*B (gemm) - Q^H*B (gemqr) | = 5.43672e-07
|Q^C*B (gemm) - Q^C*B (gemqr) | = 3.85573e-07
|B*Q (gemm) - B*Q (gemqr) | = 5.54186e-07
|B*Q^T (gemm) - B*Q^T (gemqr) | = 6.20191e-07
|B*Q^H (gemm) - B*Q^H (gemqr) | = 5.92514e-07
|B*Q^C (gemm) - B*Q^C (gemqr) | = 6.60714e-07
double
|Q*R-A| = 6.28688e-16
|Q'*Q-I| = 7.00305e-16
|Q*B (gemm) - Q*B (gemqr) | = 6.02886e-16
|Q^T*B (gemm) - Q^T*B (gemqr) | = 1.01586e-15
|Q^H*B (gemm) - Q^H*B (gemqr) | = 6.75322e-16
|Q^C*B (gemm) - Q^C*B (gemqr) | = 5.7815e-16
|B*Q (gemm) - B*Q (gemqr) | = 1.23947e-15
|B*Q^T (gemm) - B*Q^T (gemqr) | = 7.95845e-16
|B*Q^H (gemm) - B*Q^H (gemqr) | = 9.00953e-16
|B*Q^C (gemm) - B*Q^C (gemqr) | = 6.59508e-16
std::complex<float>
|Q*R-A| = 7.18633e-07
|Q'*Q-I| = 3.57662e-07
|Q*B (gemm) - Q*B (gemqr) | = 8.2653e-07
|Q^T*B (gemm) - Q^T*B (gemqr) | = 7.84967e-07
|Q^H*B (gemm) - Q^H*B (gemqr) | = 5.1556e-07
|Q^C*B (gemm) - Q^C*B (gemqr) | = 7.74009e-07
|B*Q (gemm) - B*Q (gemqr) | = 1.22603e-06
|B*Q^T (gemm) - B*Q^T (gemqr) | = 8.92469e-07
|B*Q^H (gemm) - B*Q^H (gemqr) | = 1.08471e-06
|B*Q^C (gemm) - B*Q^C (gemqr) | = 9.12118e-07
std::complex<double>
|Q*R-A| = 1.14136e-15
|Q'*Q-I| = 6.51389e-16
|Q*B (gemm) - Q*B (gemqr) | = 1.33321e-15
|Q^T*B (gemm) - Q^T*B (gemqr) | = 1.23107e-15
|Q^H*B (gemm) - Q^H*B (gemqr) | = 7.41813e-16
|Q^C*B (gemm) - Q^C*B (gemqr) | = 1.13901e-15
|B*Q (gemm) - B*Q (gemqr) | = 1.90478e-15
|B*Q^T (gemm) - B*Q^T (gemqr) | = 1.86202e-15
|B*Q^H (gemm) - B*Q^H (gemqr) | = 1.80635e-15
|B*Q^C (gemm) - B*Q^C (gemqr) | = 1.54351e-15


Go back to Summary of /test programs.

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