Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
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 
54 namespace SDMUtilsUnitTest {
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);
122  generateRandomMatrix(A, k);
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);
158  generateRandomMatrix(A, k);
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_Householder_TallSkinny ) {
243  ordinal_type m = 100;
244  ordinal_type n = 35;
245  success = test_CPQR(Stokhos::CPQR_Householder, m, n, rtol, atol, out);
246  }
247 
248  TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, CPQR_Householder_ShortFat ) {
249  ordinal_type n = 100;
250  ordinal_type m = 35;
251  success = test_CPQR(Stokhos::CPQR_Householder, m, n, rtol, atol, out);
252  }
253 
254  TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, CPQR_Householder3_TallSkinny ) {
255  ordinal_type m = 100;
256  ordinal_type n = 35;
257  success = test_CPQR(Stokhos::CPQR_Householder3, m, n, rtol, atol, out);
258  }
259 
260  TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, CPQR_Householder3_ShortFat ) {
261  ordinal_type n = 100;
262  ordinal_type m = 35;
263  success = test_CPQR(Stokhos::CPQR_Householder3, m, n, rtol, atol, out);
264  }
265 
266  TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, CPQR_Householder_threshold_TallSkinny ) {
267  ordinal_type m = 100;
268  ordinal_type n = 35;
269  ordinal_type k = 20;
271  m, n, k, rtol, atol, out);
272  }
273 
274  TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, CPQR_Householder_threshold_ShortFat ) {
275  ordinal_type n = 100;
276  ordinal_type m = 35;
277  ordinal_type k = 20;
279  m, n, k, rtol, atol, out);
280  }
281 
282  TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, CPQR_MGS_threshold_TallSkinny ) {
283  ordinal_type m = 100;
284  ordinal_type n = 35;
285  ordinal_type k = 20;
287  m, n, k, rtol, atol, out);
288  }
289 
290  TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, CPQR_MGS_threshold_ShortFat ) {
291  ordinal_type n = 100;
292  ordinal_type m = 35;
293  ordinal_type k = 20;
295  m, n, k, rtol, atol, out);
296  }
297 
298  TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, CPQR_MGS_reorthog_threshold_TallSkinny ) {
299  ordinal_type m = 100;
300  ordinal_type n = 35;
301  ordinal_type k = 20;
303  m, n, k, rtol, atol, out);
304  }
305 
306  TEUCHOS_UNIT_TEST( Stokhos_SDMUtils, CPQR_MGS_reorthog_threshold_ShortFat ) {
307  ordinal_type n = 100;
308  ordinal_type m = 35;
309  ordinal_type k = 20;
311  m, n, k, rtol, atol, out);
312  }
313 
314 }
315 
316 int main( int argc, char* argv[] ) {
317  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
318  return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
319 }
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.
ordinal_type(* wcpqr_func_type)(const scalar_type &, const SDM &, const Teuchos::Array< scalar_type > &, SDM &, SDM &, Teuchos::Array< ordinal_type > &)
void(* cpqr_func_type)(const SDM &, SDM &, SDM &, Teuchos::Array< ordinal_type > &)
int main(int argc, char *argv[])
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 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_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 generateRandomMatrix(SDM &A, ordinal_type rank)
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)
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.
bool test_CPQR(cpqr_func_type qr_func, ordinal_type m, ordinal_type n, scalar_type rtol, scalar_type atol, Teuchos::FancyOStream &out)
KOKKOS_INLINE_FUNCTION PCE< Storage > min(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
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)
expr expr expr expr j
bool test_QR(qr_func_type qr_func, ordinal_type m, ordinal_type n, scalar_type rtol, scalar_type atol, Teuchos::FancyOStream &out)
void(* qr_func_type)(ordinal_type, const SDM &, const Teuchos::Array< scalar_type > &, SDM &, SDM &)
void CPQR_Householder(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_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.
Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > SDM
TEUCHOS_UNIT_TEST(Stokhos_SDMUtils, QR_CGS_TallSkinny)
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.
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.