Teuchos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
numerics/test/LAPACK/cxx_main.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Teuchos: Common Tools Package
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#include <iostream>
43#include <vector>
44#include "Teuchos_LAPACK.hpp"
45#include "Teuchos_Version.hpp"
48
49
50template <typename T>
51int lapackTest( bool verbose );
52
53template <typename T>
55{
56 // Specialized test for real-valued vs. complex valued types
57 static int test(bool verbose);
58
59 // Specializations for mixed complex-real operations
60 template<class R>
61 void add( const R& a, const T& b, T& result ){ result = a + b; }
62
63 template<class R>
64 void multiply( const R& a, const T& b, T& result ){ result = a * b; }
65
66 template<class R>
67 void divide( const T& a, const R& b, T& result ){ result = a / b; }
68};
69
70#ifdef HAVE_TEUCHOS_COMPLEX
71
72// Partial specialization for std::complex numbers templated on real type T
73template <typename T>
74struct specializedLAPACK< std::complex<T> >
75{
76 // Specialized test for real-valued vs. complex valued types
77 static int test(bool verbose);
78
79 // Specializations for mixed complex-real operations
80 template<class R>
81 void add( const R& a, const std::complex<T>& b, std::complex<T>& result )
82 { std::complex<T> tmp( a, 0 ); result = b + tmp; }
83
84 template<class R>
85 void multiply( const R& a, const std::complex<T>& b, std::complex<T>& result )
86 { std::complex<T> tmp( a, 0 ); result = tmp * b; }
87
88 template<class R>
89 void divide( const std::complex<T>& a, const R& b, std::complex<T>& result )
90 { std::complex<T> tmp( b, 0 ); result = a / tmp; }
91};
92
93#endif
94
95// Main test
96int main(int argc, char* argv[])
97{
98 int numberFailedTests = 0;
99 bool verbose = 0;
100 if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
101
102 if (verbose)
103 std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
104
105 using std::fabs;
106
107#ifdef HAVE_TEUCHOS_INST_FLOAT
108 if (verbose)
109 std::cout << std::endl << "LAPACK test for float" << std::endl;
110 numberFailedTests += lapackTest<float>(verbose);
111#endif
112 if (verbose)
113 std::cout << std::endl << "LAPACK test for double" << std::endl;
114 numberFailedTests += lapackTest<double>(verbose);
115
116#ifdef HAVE_TEUCHOS_COMPLEX
117#ifdef HAVE_TEUCHOS_INST_COMPLEX_FLOAT
118 if (verbose)
119 std::cout << std::endl << "LAPACK test for std::complex<float>" << std::endl;
120 numberFailedTests += lapackTest<std::complex<float> >(verbose);
121#endif
122#ifdef HAVE_TEUCHOS_INST_COMPLEX_DOUBLE
123 if (verbose)
124 std::cout << std::endl << "LAPACK test for std::complex<double>" << std::endl;
125 numberFailedTests += lapackTest<std::complex<double> >(verbose);
126#endif
127#endif
128
129 if(numberFailedTests > 0)
130 {
131 if (verbose) {
132 std::cout << "Number of failed tests: " << numberFailedTests << std::endl;
133 std::cout << "End Result: TEST FAILED" << std::endl;
134 return -1;
135 }
136 }
137 if(numberFailedTests==0)
138 std::cout << "End Result: TEST PASSED" << std::endl;
139 return 0;
140}
141
142// Common test for all four types: float, double, std::complex<float>, std::complex<double>
143// Calls the specialized test for types whose interfaces are different or undefined
144template <typename T>
145int lapackTest( bool verbose )
146{
147 int numberFailedTests = 0;
148
149 // Define some common characters
150 int info=0;
151 char char_G = 'G';
152 char char_N = 'N';
153 char char_U = 'U';
154
155 // Create some common typedefs
156 typedef Teuchos::ScalarTraits<T> STS;
157 typedef typename STS::magnitudeType MagnitudeType;
159
160 T one = STS::one();
161 MagnitudeType m_one = STM::one();
162 T zero = STS::zero();
163
166
167 const int n_gesv = 4;
168 std::vector<T> Ad(n_gesv*n_gesv,zero);
169 std::vector<T> bd(n_gesv,zero);
170 int IPIV[n_gesv];
171
172 Ad[0] = 1; Ad[2] = 1; Ad[5] = 1; Ad[8] = 2; Ad[9] = 1; Ad[10] = 1; Ad[14] = 2; Ad[15] = 2;
173 bd[1] = 2; bd[2] = 1; bd[3] = 2;
174
175 if (verbose) std::cout << "LASCL test ... ";
176 L.LASCL(char_G, 1, 1, m_one, m_one, n_gesv, n_gesv, &Ad[0], n_gesv, &info);
177 if ( !info ) {
178 if (verbose) std::cout << "passed!" << std::endl;
179 } else {
180 if (verbose) std::cout << "FAILED" << std::endl;
181 numberFailedTests++;
182 }
183
184 if (verbose) std::cout << "GESV test ... ";
185 L.GESV(n_gesv, 1, &Ad[0], n_gesv, IPIV, &bd[0], n_gesv, &info);
186 if ( !info ) {
187 if (verbose) std::cout << "passed!" << std::endl;
188 } else {
189 if (verbose) std::cout << "FAILED" << std::endl;
190 numberFailedTests++;
191 }
192
193#if ! (defined(__INTEL_COMPILER) && defined(_WIN32) )
194
195 // Check ILAENV with similarity transformation routine: dsytrd
196 // NOTE: Do not need to put floating point specifier [s,d,c,z] before routine name,
197 // this is handled through templating.
198 if (verbose) std::cout << "ILAENV test ... ";
199 int n1 = 100;
200 int size = L.ILAENV(1, "sytrd", "u", n1);
201 if (size > 0) {
202 if (verbose) std::cout << "passed!" << std::endl;
203 } else {
204 if (verbose) std::cout << "FAILED!" << std::endl;
205 numberFailedTests++;
206 }
207
208#endif
209
210 // Create a simple diagonal linear system
211 std::vector<T> Ad2_sub(n_gesv-1, zero), b2(n_gesv, one);
212 std::vector<MagnitudeType> Ad2(n_gesv, m_one);
213
214 if (verbose) std::cout << "PTTRF test ... ";
215 L.PTTRF(n_gesv, &Ad2[0], &Ad2_sub[0], &info);
216 if ( !info ) {
217 if (verbose) std::cout << "passed!" << std::endl;
218 } else {
219 if (verbose) std::cout << "FAILED" << std::endl;
220 numberFailedTests++;
221 }
222
223 int n_potrf = 5;
224 std::vector<T> diag_a(n_potrf*n_potrf, zero);
225 for (int i=0; i<n_potrf; i++)
226 {
227 T tmp = zero;
228 sL.add( i, one, tmp );
229 diag_a[i*n_potrf + i] = tmp*tmp;
230 }
231
232 if (verbose) std::cout << "POTRF test ... ";
233 L.POTRF(char_U, n_potrf, &diag_a[0], n_potrf, &info);
234
235 if (info != 0)
236 {
237 if (verbose) std::cout << "FAILED" << std::endl;
238 numberFailedTests++;
239 }
240 else
241 {
242 for (int i=0; i<n_potrf; i++)
243 {
244 T tmp = zero;
245 sL.add( i, one, tmp );
246 if ( diag_a[i*n_potrf + i] == tmp )
247 {
248 if (verbose && i==(n_potrf-1)) std::cout << "passed!" << std::endl;
249 }
250 else
251 {
252 if (verbose) std::cout << "FAILED" << std::endl;
253 numberFailedTests++;
254 break;
255 }
256 }
257 }
258
259 if (verbose) std::cout << "POTRI test ... ";
260 std::vector<T> diag_a_trtri(diag_a); // Save a copy for TRTRI test
261
262 L.POTRI(char_U, n_potrf, &diag_a[0], n_potrf, &info);
263
264 T tmp = zero;
265 sL.multiply( 1.0/4.0, one, tmp );
266 if ( info != 0 || (diag_a[n_potrf+1] != tmp) )
267 {
268 if (verbose) std::cout << "FAILED" << std::endl;
269 numberFailedTests++;
270 }
271 else
272 if (verbose) std::cout << "passed!" << std::endl;
273
274 if (verbose) std::cout << "TRTRI test ... ";
275
276 int n_trtri = n_potrf;
277 L.TRTRI( char_U, char_N, n_trtri, &diag_a_trtri[0], n_trtri, &info );
278 for (int i=0; i<n_trtri; i++)
279 {
280 tmp = zero;
281 sL.divide( one, i+1.0, tmp );
282 if ( info != 0 )
283 {
284 numberFailedTests++;
285 break;
286 }
287 else if ( diag_a_trtri[i*n_trtri + i] == tmp )
288 {
289 if (verbose && i==(n_trtri-1)) std::cout << "passed!" << std::endl;
290 }
291 else
292 {
293 if (verbose) std::cout << "FAILED" << std::endl;
294 numberFailedTests++;
295 break;
296 }
297 }
298
299#ifndef TEUCHOSNUMERICS_DISABLE_STEQR_TEST
300
301 if (verbose) std::cout << "STEQR test ... ";
302
303 const int n_steqr = 10;
304 std::vector<MagnitudeType> diagonal(n_steqr);
305 std::vector<MagnitudeType> subdiagonal(n_steqr-1);
306
307 for (int i=0; i < n_steqr; ++i) {
308 diagonal[i] = n_steqr - i;
309 if (i < n_steqr-1)
310 subdiagonal[i] = STM::eps() * (i+1);
311 }
312
313 std::vector<T> scalar_dummy(1,0.0);
314 std::vector<MagnitudeType> mag_dummy(4*n_steqr,0.0);
315
316 L.STEQR (char_N, n_steqr, &diagonal[0], &subdiagonal[0],
317 &scalar_dummy[0], n_steqr, &mag_dummy[0], &info);
318
319 if (info != 0)
320 {
321 if (verbose) std::cout << "STEQR: compute symmetric tridiagonal eigenvalues: "
322 << "LAPACK's _STEQR failed with info = "
323 << info;
324
325 numberFailedTests++;
326 }
327
328 MagnitudeType lambda_min = diagonal[0];
329 MagnitudeType lambda_max = diagonal[n_steqr-1];
330 MagnitudeType exp_lambda_min = STM::one();
331 MagnitudeType exp_lambda_max = STM::one()*n_steqr;
332
333 if ((fabs(lambda_min-exp_lambda_min)<1e-12) && (fabs(lambda_max-exp_lambda_max)<1e-12))
334 {
335 if (verbose) std::cout << "passed!" << std::endl;
336 }
337 else
338 {
339 if (verbose) std::cout << "FAILED" << std::endl;
340 numberFailedTests++;
341 }
342
343#endif // TEUCHOSNUMERICS_DISABLE_STEQR_TEST
344
345 numberFailedTests += specializedLAPACK<T>::test( verbose );
346
347 return numberFailedTests;
348}
349
350template<class T>
352{
353 // Create some common typedefs
354 typedef Teuchos::ScalarTraits<T> STS;
355 typedef typename STS::magnitudeType MagnitudeType;
357
358 T one = STS::one();
359 MagnitudeType m_one = STM::one();
360 T zero = STS::zero();
361
362 char char_E = 'E';
363 char char_U = 'U';
364
365 int info=0;
366 int numberFailedTests = 0;
367
369
370 if (verbose) std::cout << "LAPY2 test ... ";
371 T x = 3*one, y = 4*one;
372 T lapy = L.LAPY2(x, y);
373 if ( lapy == 5*one ) {
374 if (verbose) std::cout << "passed!" << std::endl;
375 } else {
376 if (verbose) std::cout << "FAILED ( " << lapy << " != 5 )" << std::endl;
377 numberFailedTests++;
378 }
379
380 if (verbose) std::cout << "LAMCH test ... ";
381
382 T st_eps = L.LAMCH( char_E );
383 if (verbose)
384 std::cout << "[ eps = " << st_eps << " ] passed!" << std::endl;
385
386 // Create a simple diagonal linear system
387 const int n = 4;
388 std::vector<T> Ad2_sub(n-1, zero), b2(n, one);
389 std::vector<MagnitudeType> Ad2(n, m_one);
390
391 if (verbose) std::cout << "PTTRS test ... ";
392 L.PTTRS(n, 1, &Ad2[0], &Ad2_sub[0], &b2[0], n, &info);
393 if ( !info ) {
394 if (verbose) std::cout << "passed!" << std::endl;
395 } else {
396 if (verbose) std::cout << "FAILED" << std::endl;
397 numberFailedTests++;
398 }
399
400 if (verbose) std::cout << "POCON test ... ";
401
402 std::vector<T> diag_a(n*n);
403 for (int i=0; i<n; i++)
404 {
405 diag_a[i*n + i] = one;
406 }
407 MagnitudeType rcond, anorm = m_one;
408 std::vector<T> work(3*n);
409 std::vector<int> iwork(n);
410
411 L.POCON(char_U, n, &diag_a[0], n, anorm, &rcond, &work[0], &iwork[0], &info);
412 if (info != 0 || (rcond != m_one))
413 {
414 if (verbose) std::cout << "FAILED" << std::endl;
415 numberFailedTests++;
416 }
417 else
418 if (verbose) std::cout << "passed!" << std::endl;
419
420
421 return numberFailedTests;
422}
423
424#ifdef HAVE_TEUCHOS_COMPLEX
425
426template<class T>
427int specializedLAPACK<std::complex<T> >::test( bool verbose )
428{
429 // Create some common typedefs
431 typedef typename STS::magnitudeType MagnitudeType;
433
434 std::complex<T> one = STS::one();
435 MagnitudeType m_one = STM::one();
436 std::complex<T> zero = STS::zero();
437
438 char char_L = 'L';
439 char char_U = 'U';
440
441 int info=0;
442 int numberFailedTests = 0;
443
445
446 // Create a simple diagonal linear system
447 const int n = 4;
448 std::vector<std::complex<T> > Ad2_sub(n-1, zero), b2(n, one);
449 std::vector<MagnitudeType> Ad2(n, m_one);
450
451 if (verbose) std::cout << "PTTRS test ... ";
452 L.PTTRS(char_L, n, 1, &Ad2[0], &Ad2_sub[0], &b2[0], n, &info);
453 if ( !info ) {
454 if (verbose) std::cout << "passed!" << std::endl;
455 } else {
456 if (verbose) std::cout << "FAILED" << std::endl;
457 numberFailedTests++;
458 }
459
460 if (verbose) std::cout << "POCON test ... ";
461
462 std::vector<std::complex<T> > diag_a(n*n);
463 for (int i=0; i<n; i++)
464 {
465 diag_a[i*n + i] = one;
466 }
467 MagnitudeType rcond, anorm = m_one;
468 std::vector<std::complex<T> > work(2*n);
469 std::vector<MagnitudeType> rwork(n);
470 std::vector<int> iwork(n);
471
472 L.POCON(char_U, n, &diag_a[0], n, anorm, &rcond, &work[0], &rwork[0], &info);
473 if (info != 0 || (rcond != m_one))
474 {
475 if (verbose) std::cout << "FAILED" << std::endl;
476 numberFailedTests++;
477 }
478 else
479 if (verbose) std::cout << "passed!" << std::endl;
480
481return numberFailedTests;
482}
483
484#endif
Templated interface class to LAPACK routines.
Templated serial dense matrix class.
Templated serial dense vector class.
The Templated LAPACK Wrapper Class.
void STEQR(const char &COMPZ, const OrdinalType &n, MagnitudeType *D, MagnitudeType *E, ScalarType *Z, const OrdinalType &ldz, MagnitudeType *WORK, OrdinalType *info) const
Computes the eigenvalues and, optionally, eigenvectors of a symmetric tridiagonal n by n matrix A usi...
ScalarType LAMCH(const char &CMACH) const
Determines machine parameters for floating point characteristics.
void PTTRS(const OrdinalType &n, const OrdinalType &nrhs, const MagnitudeType *d, const ScalarType *e, ScalarType *B, const OrdinalType &ldb, OrdinalType *info) const
Solves a tridiagonal system A*X=B using the \L*D*L' factorization of A computed by PTTRF.
void POCON(const char &UPLO, const OrdinalType &n, const ScalarType *A, const OrdinalType &lda, const ScalarType &anorm, ScalarType *rcond, ScalarType *WORK, OrdinalType *IWORK, OrdinalType *info) const
Estimates the reciprocal of the condition number (1-norm) of a real symmetric positive definite matri...
ScalarType LAPY2(const ScalarType &x, const ScalarType &y) const
Computes x^2 + y^2 safely, to avoid overflow.
void POTRI(const char &UPLO, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, OrdinalType *info) const
Computes the inverse of a real symmetric positive definite matrix A using the Cholesky factorization ...
void PTTRF(const OrdinalType &n, MagnitudeType *d, ScalarType *e, OrdinalType *info) const
Computes the L*D*L' factorization of a Hermitian/symmetric positive definite tridiagonal matrix A.
void GESV(const OrdinalType &n, const OrdinalType &nrhs, ScalarType *A, const OrdinalType &lda, OrdinalType *IPIV, ScalarType *B, const OrdinalType &ldb, OrdinalType *info) const
Computes the solution to a real system of linear equations A*X=B, where A is factored through GETRF a...
void TRTRI(const char &UPLO, const char &DIAG, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, OrdinalType *info) const
Computes the inverse of an upper or lower triangular matrix A.
OrdinalType ILAENV(const OrdinalType &ispec, const std::string &NAME, const std::string &OPTS, const OrdinalType &N1=-1, const OrdinalType &N2=-1, const OrdinalType &N3=-1, const OrdinalType &N4=-1) const
Chooses problem-dependent parameters for the local environment.
void LASCL(const char &TYPE, const OrdinalType &kl, const OrdinalType &ku, const MagnitudeType cfrom, const MagnitudeType cto, const OrdinalType &m, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, OrdinalType *info) const
Multiplies the m by n matrix A by the real scalar cto/cfrom.
void POTRF(const char &UPLO, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, OrdinalType *info) const
Computes Cholesky factorization of a real symmetric positive definite matrix A.
int main()
Definition: evilMain.cpp:75
std::string Teuchos_Version()
int lapackTest(bool verbose)
This structure defines some basic traits for a scalar field type.
static int test(bool verbose)
void multiply(const R &a, const T &b, T &result)
void add(const R &a, const T &b, T &result)
void divide(const T &a, const R &b, T &result)