Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_SDMUtilsUnitTest.cpp
Go to the documentation of this file.
1// $Id$
2// $Source$
3// @HEADER
4// ***********************************************************************
5//
6// Stokhos Package
7// Copyright (2009) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
40//
41// ***********************************************************************
42// @HEADER
43
44#include "Teuchos_UnitTestHarness.hpp"
45#include "Teuchos_TestingHelpers.hpp"
46#include "Teuchos_UnitTestRepository.hpp"
47#include "Teuchos_GlobalMPISession.hpp"
48
49#include "Stokhos_SDMUtils.hpp"
51
52// Need to test pre-pivot
53
55
56 typedef int ordinal_type;
57 typedef double scalar_type;
58 typedef Teuchos::SerialDenseMatrix<ordinal_type,scalar_type> SDM;
59 typedef void (*qr_func_type)(ordinal_type, const SDM&, const Teuchos::Array<scalar_type>&, SDM&, SDM&);
60 typedef void (*cpqr_func_type)(const SDM&, SDM&, SDM&, Teuchos::Array<ordinal_type>&);
61 typedef ordinal_type (*wcpqr_func_type)(const scalar_type&, const SDM&, const Teuchos::Array<scalar_type>&, SDM&, SDM&, Teuchos::Array<ordinal_type>&);
62
63 scalar_type rtol = 1.0e-12;
64 scalar_type atol = 1.0e-12;
65
67 ordinal_type m = A.numRows();
68 ordinal_type n = A.numCols();
69 ordinal_type k = std::min(m,n);
70 SDM B(m,m), C(n,n), S(m,n), T(m,n);
71 B.random();
72 C.random();
73 S.putScalar(0.0);
74 if (rank > k)
75 rank = k;
76 for (ordinal_type i=0; i<rank; i++)
77 S(i,i) = 1.0;
78 T.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, B, S, 0.0);
79 A.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, T, C, 0.0);
80 }
81
84 Teuchos::FancyOStream& out) {
85 bool success;
86
87 SDM A(m,n);
88 A.random();
89 SDM Q, R;
90 Teuchos::Array<scalar_type> w(m, 1.0);
91 ordinal_type k = std::min(m,n);
92 qr_func(k, A, w, Q, R);
93
94 TEUCHOS_TEST_EQUALITY(Q.numRows(), m, out, success);
95 TEUCHOS_TEST_EQUALITY(Q.numCols(), k, out, success);
96 TEUCHOS_TEST_EQUALITY(R.numRows(), k, out, success);
97 TEUCHOS_TEST_EQUALITY(R.numCols(), k, out, success);
98
99 // Test A = Q*R
100 SDM QR(m,k);
101 QR.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Q, R, 0.0);
102 SDM AA(Teuchos::View, A, m, k);
103 success = Stokhos::compareSDM(AA, "A", QR, "Q*R", rtol, atol, out);
104
105 // Test Q^T*Q = I
106 SDM eye(k,k), QTQ(k,k);
107 for (ordinal_type i=0; i<k; i++)
108 eye(i,i) = 1.0;
109 QTQ.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Q, Q, 0.0);
110 success = Stokhos::compareSDM(eye, "eye", QTQ, "Q^T*Q", rtol, atol, out);
111
112 return success;
113 }
114
117 Teuchos::FancyOStream& out) {
118 bool success;
119
120 SDM A(m,n);
121 ordinal_type k = std::min(m,n);
123 SDM Q, R;
124 Teuchos::Array<ordinal_type> piv;
125 qr_func(A, Q, R, piv);
126
127 TEUCHOS_TEST_EQUALITY(Q.numRows(), m, out, success);
128 TEUCHOS_TEST_EQUALITY(Q.numCols(), k, out, success);
129 TEUCHOS_TEST_EQUALITY(R.numRows(), k, out, success);
130 TEUCHOS_TEST_EQUALITY(R.numCols(), n, out, success);
131 TEUCHOS_TEST_EQUALITY(piv.size(), n, out, success);
132
133 // Test A*P = Q*R
134 SDM AP(m,n), QR(m,n);
135 for (ordinal_type j=0; j<n; j++)
136 for (ordinal_type i=0; i<m; i++)
137 AP(i,j) = A(i,piv[j]);
138 QR.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Q, R, 0.0);
139 success = Stokhos::compareSDM(AP, "A*P", QR, "Q*R", rtol, atol, out);
140
141 // Test Q^T*Q = I
142 SDM eye(k,k), QTQ(k,k);
143 for (ordinal_type i=0; i<k; i++)
144 eye(i,i) = 1.0;
145 QTQ.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Q, Q, 0.0);
146 success = Stokhos::compareSDM(eye, "eye", QTQ, "Q^T*Q", rtol, atol, out);
147
148 return success;
149 }
150
154 Teuchos::FancyOStream& out) {
155 bool success;
156
157 SDM A(m,n);
159 SDM Q, R;
160 Teuchos::Array<ordinal_type> piv(n);
161 for (ordinal_type i=0; i<5; i++)
162 piv[i] = 1;
163 Teuchos::Array<scalar_type> w(m, 1.0);
164 ordinal_type rank = qr_func(1.0e-12, A, w, Q, R, piv);
165
166 TEUCHOS_TEST_EQUALITY(rank, k, out, success);
167 TEUCHOS_TEST_EQUALITY(Q.numRows(), m, out, success);
168 TEUCHOS_TEST_EQUALITY(Q.numCols(), k, out, success);
169 TEUCHOS_TEST_EQUALITY(R.numRows(), k, out, success);
170 TEUCHOS_TEST_EQUALITY(R.numCols(), k, out, success);
171 TEUCHOS_TEST_EQUALITY(piv.size(), n, out, success);
172
173 // Test A*P = Q*R
174 SDM AP(m,k), QR(m,k);
175 for (ordinal_type j=0; j<k; j++)
176 for (ordinal_type i=0; i<m; i++)
177 AP(i,j) = A(i,piv[j]);
178 QR.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Q, R, 0.0);
179 success = Stokhos::compareSDM(AP, "A*P", QR, "Q*R", rtol, atol, out);
180
181 // Test Q^T*Q = I
182 SDM eye(k,k), Qt(m,k), QTQ(k,k);
183 for (ordinal_type j=0; j<k; j++) {
184 eye(j,j) = 1.0;
185 for (ordinal_type i=0; i<m; i++)
186 Qt(i,j) = w[i]*Q(i,j);
187 }
188 QTQ.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Qt, Q, 0.0);
189 success = Stokhos::compareSDM(eye, "eye", QTQ, "Q^T*Q", rtol, atol, out);
190
191 return success;
192 }
193
194 TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, QR_CGS_TallSkinny ) {
195 ordinal_type m = 100;
196 ordinal_type n = 35;
197 success = test_QR(Stokhos::QR_CGS, m, n, rtol, atol, out);
198 }
199
200 TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, QR_CGS_ShortFat ) {
201 ordinal_type n = 100;
202 ordinal_type m = 35;
203 success = test_QR(Stokhos::QR_CGS, m, n, rtol, atol, out);
204 }
205
206 TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, QR_MGS_TallSkinny ) {
207 ordinal_type m = 100;
208 ordinal_type n = 35;
209 success = test_QR(Stokhos::QR_MGS, m, n, rtol, atol, out);
210 }
211
212 TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, QR_MGS_ShortFat ) {
213 ordinal_type n = 100;
214 ordinal_type m = 35;
215 success = test_QR(Stokhos::QR_MGS, m, n, rtol, atol, out);
216 }
217
218 TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, QR_MGS2_TallSkinny ) {
219 ordinal_type m = 100;
220 ordinal_type n = 35;
221 success = test_QR(Stokhos::QR_MGS2, m, n, rtol, atol, out);
222 }
223
224 TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, QR_MGS2_ShortFat ) {
225 ordinal_type n = 100;
226 ordinal_type m = 35;
227 success = test_QR(Stokhos::QR_MGS2, m, n, rtol, atol, out);
228 }
229
230 TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, QR_Householder_TallSkinny ) {
231 ordinal_type m = 100;
232 ordinal_type n = 35;
233 success = test_QR(Stokhos::QR_Householder, m, n, rtol, atol, out);
234 }
235
236 TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, QR_Householder_ShortFat ) {
237 ordinal_type n = 100;
238 ordinal_type m = 35;
239 success = test_QR(Stokhos::QR_Householder, m, n, rtol, atol, out);
240 }
241
242 TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, CPQR_Householder3_TallSkinny ) {
243 ordinal_type m = 100;
244 ordinal_type n = 35;
245 success = test_CPQR(Stokhos::CPQR_Householder3, m, n, rtol, atol, out);
246 }
247
248 TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, CPQR_Householder3_ShortFat ) {
249 ordinal_type n = 100;
250 ordinal_type m = 35;
251 success = test_CPQR(Stokhos::CPQR_Householder3, m, n, rtol, atol, out);
252 }
253
254 TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, CPQR_Householder_threshold_TallSkinny ) {
255 ordinal_type m = 100;
256 ordinal_type n = 35;
257 ordinal_type k = 20;
259 m, n, k, rtol, atol, out);
260 }
261
262 TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, CPQR_Householder_threshold_ShortFat ) {
263 ordinal_type n = 100;
264 ordinal_type m = 35;
265 ordinal_type k = 20;
267 m, n, k, rtol, atol, out);
268 }
269
270 TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, CPQR_MGS_threshold_TallSkinny ) {
271 ordinal_type m = 100;
272 ordinal_type n = 35;
273 ordinal_type k = 20;
275 m, n, k, rtol, atol, out);
276 }
277
278 TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, CPQR_MGS_threshold_ShortFat ) {
279 ordinal_type n = 100;
280 ordinal_type m = 35;
281 ordinal_type k = 20;
283 m, n, k, rtol, atol, out);
284 }
285
286 TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, CPQR_MGS_reorthog_threshold_TallSkinny ) {
287 ordinal_type m = 100;
288 ordinal_type n = 35;
289 ordinal_type k = 20;
291 m, n, k, rtol, atol, out);
292 }
293
294 TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, CPQR_MGS_reorthog_threshold_ShortFat ) {
295 ordinal_type n = 100;
296 ordinal_type m = 35;
297 ordinal_type k = 20;
299 m, n, k, rtol, atol, out);
300 }
301
302}
303
304int main( int argc, char* argv[] ) {
305 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
306 return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
307}
int main(int argc, char *argv[])
void generateRandomMatrix(SDM &A, ordinal_type rank)
bool test_QR(qr_func_type qr_func, ordinal_type m, ordinal_type n, scalar_type rtol, scalar_type atol, Teuchos::FancyOStream &out)
bool test_CPQR(cpqr_func_type qr_func, ordinal_type m, ordinal_type n, scalar_type rtol, scalar_type atol, Teuchos::FancyOStream &out)
TEUCHOS_UNIT_TEST(Stokhos_SDMUtils, QR_CGS_TallSkinny)
bool test_weighted_CPQR(wcpqr_func_type qr_func, ordinal_type m, ordinal_type n, ordinal_type k, scalar_type rtol, scalar_type atol, Teuchos::FancyOStream &out)
Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > SDM
void(* cpqr_func_type)(const SDM &, SDM &, SDM &, Teuchos::Array< ordinal_type > &)
void(* qr_func_type)(ordinal_type, const SDM &, const Teuchos::Array< scalar_type > &, SDM &, SDM &)
ordinal_type(* wcpqr_func_type)(const scalar_type &, const SDM &, const Teuchos::Array< scalar_type > &, SDM &, SDM &, Teuchos::Array< ordinal_type > &)
ordinal_type CPQR_MGS_threshold(const scalar_type &rank_threshold, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::Array< scalar_type > &w, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R, Teuchos::Array< ordinal_type > &piv)
Compute column-pivoted QR using modified Gram-Schmidt.
bool compareSDM(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &a1, const std::string &a1_name, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &a2, const std::string &a2_name, const scalar_type &rel_tol, const scalar_type &abs_tol, Teuchos::FancyOStream &out)
ordinal_type CPQR_MGS_reorthog_threshold(const scalar_type &rank_threshold, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::Array< scalar_type > &w, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R, Teuchos::Array< ordinal_type > &piv)
Compute column-pivoted QR using modified Gram-Schmidt and reorthogonalization.
ordinal_type CPQR_Householder_threshold(const scalar_type &rank_threshold, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::Array< scalar_type > &w, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R, Teuchos::Array< ordinal_type > &piv)
Compute column-pivoted QR using Householder reflections.
void CPQR_Householder3(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R, Teuchos::Array< ordinal_type > &piv)
Compute column-pivoted QR using Householder reflections.
void QR_MGS(ordinal_type k, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::Array< scalar_type > &w, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R)
Compute thin QR using modified Gram-Schmidt.
void QR_MGS2(ordinal_type k, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::Array< scalar_type > &w, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R)
Compute thin QR using modified Gram-Schmidt with reorthogonalization.
void QR_CGS(ordinal_type k, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::Array< scalar_type > &w, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R)
Compute thin QR using classical Gram-Schmidt.
void QR_Householder(ordinal_type k, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::Array< scalar_type > &w, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R)
Compute thin QR using Householder reflections.