Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
numerics/test/LAPACK/cxx_main.cpp
Go to the documentation of this file.
2 // *****************************************************************************
3 // Teuchos: Common Tools Package
4 //
5 // Copyright 2004 NTESS and the Teuchos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
9
10 #include <iostream>
11 #include <vector>
12 #include "Teuchos_LAPACK.hpp"
13 #include "Teuchos_Version.hpp"
16
17
18 template <typename T>
19 int lapackTest( bool verbose );
20
21 template <typename T>
23 {
24  // Specialized test for real-valued vs. complex valued types
25  static int test(bool verbose);
26
27  // Specializations for mixed complex-real operations
28  template<class R>
29  void add( const R& a, const T& b, T& result ){ result = a + b; }
30
31  template<class R>
32  void multiply( const R& a, const T& b, T& result ){ result = a * b; }
33
34  template<class R>
35  void divide( const T& a, const R& b, T& result ){ result = a / b; }
36 };
37
38 #ifdef HAVE_TEUCHOS_COMPLEX
39
40 // Partial specialization for std::complex numbers templated on real type T
41 template <typename T>
42 struct specializedLAPACK< std::complex<T> >
43 {
44  // Specialized test for real-valued vs. complex valued types
45  static int test(bool verbose);
46
47  // Specializations for mixed complex-real operations
48  template<class R>
49  void add( const R& a, const std::complex<T>& b, std::complex<T>& result )
50  { std::complex<T> tmp( a, 0 ); result = b + tmp; }
51
52  template<class R>
53  void multiply( const R& a, const std::complex<T>& b, std::complex<T>& result )
54  { std::complex<T> tmp( a, 0 ); result = tmp * b; }
55
56  template<class R>
57  void divide( const std::complex<T>& a, const R& b, std::complex<T>& result )
58  { std::complex<T> tmp( b, 0 ); result = a / tmp; }
59 };
60
61 #endif
62
63 // Main test
64 int main(int argc, char* argv[])
65 {
66  int numberFailedTests = 0;
67  bool verbose = 0;
68  if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
69
70  if (verbose)
71  std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
72
73  using std::fabs;
74
75 #ifdef HAVE_TEUCHOS_INST_FLOAT
76  if (verbose)
77  std::cout << std::endl << "LAPACK test for float" << std::endl;
78  numberFailedTests += lapackTest<float>(verbose);
79 #endif
80  if (verbose)
81  std::cout << std::endl << "LAPACK test for double" << std::endl;
82  numberFailedTests += lapackTest<double>(verbose);
83
84 #ifdef HAVE_TEUCHOS_COMPLEX
85 #ifdef HAVE_TEUCHOS_INST_COMPLEX_FLOAT
86  if (verbose)
87  std::cout << std::endl << "LAPACK test for std::complex<float>" << std::endl;
88  numberFailedTests += lapackTest<std::complex<float> >(verbose);
89 #endif
90 #ifdef HAVE_TEUCHOS_INST_COMPLEX_DOUBLE
91  if (verbose)
92  std::cout << std::endl << "LAPACK test for std::complex<double>" << std::endl;
93  numberFailedTests += lapackTest<std::complex<double> >(verbose);
94 #endif
95 #endif
96
97  if(numberFailedTests > 0)
98  {
99  if (verbose) {
100  std::cout << "Number of failed tests: " << numberFailedTests << std::endl;
101  std::cout << "End Result: TEST FAILED" << std::endl;
102  return -1;
103  }
104  }
105  if(numberFailedTests==0)
106  std::cout << "End Result: TEST PASSED" << std::endl;
107  return 0;
108 }
109
110 // Common test for all four types: float, double, std::complex<float>, std::complex<double>
111 // Calls the specialized test for types whose interfaces are different or undefined
112 template <typename T>
113 int lapackTest( bool verbose )
114 {
115  int numberFailedTests = 0;
116
117  // Define some common characters
118  int info=0;
119  char char_G = 'G';
120  char char_N = 'N';
121  char char_U = 'U';
122
123  // Create some common typedefs
124  typedef Teuchos::ScalarTraits<T> STS;
125  typedef typename STS::magnitudeType MagnitudeType;
127
128  T one = STS::one();
129  MagnitudeType m_one = STM::one();
130  T zero = STS::zero();
131
134
135  const int n_gesv = 4;
137  std::vector<T> bd(n_gesv,zero);
138  int IPIV[n_gesv];
139
141  bd[1] = 2; bd[2] = 1; bd[3] = 2;
142
143  if (verbose) std::cout << "LASCL test ... ";
144  L.LASCL(char_G, 1, 1, m_one, m_one, n_gesv, n_gesv, &Ad[0], n_gesv, &info);
145  if ( !info ) {
146  if (verbose) std::cout << "passed!" << std::endl;
147  } else {
148  if (verbose) std::cout << "FAILED" << std::endl;
149  numberFailedTests++;
150  }
151
152  if (verbose) std::cout << "GESV test ... ";
153  L.GESV(n_gesv, 1, &Ad[0], n_gesv, IPIV, &bd[0], n_gesv, &info);
154  if ( !info ) {
155  if (verbose) std::cout << "passed!" << std::endl;
156  } else {
157  if (verbose) std::cout << "FAILED" << std::endl;
158  numberFailedTests++;
159  }
160
161 #if ! (defined(__INTEL_COMPILER) && defined(_WIN32) )
162
163  // Check ILAENV with similarity transformation routine: dsytrd
164  // NOTE: Do not need to put floating point specifier [s,d,c,z] before routine name,
165  // this is handled through templating.
166  if (verbose) std::cout << "ILAENV test ... ";
167  int n1 = 100;
168  int size = L.ILAENV(1, "sytrd", "u", n1);
169  if (size > 0) {
170  if (verbose) std::cout << "passed!" << std::endl;
171  } else {
172  if (verbose) std::cout << "FAILED!" << std::endl;
173  numberFailedTests++;
174  }
175
176 #endif
177
178  // Create a simple diagonal linear system
179  std::vector<T> Ad2_sub(n_gesv-1, zero), b2(n_gesv, one);
180  std::vector<MagnitudeType> Ad2(n_gesv, m_one);
181
182  if (verbose) std::cout << "PTTRF test ... ";
184  if ( !info ) {
185  if (verbose) std::cout << "passed!" << std::endl;
186  } else {
187  if (verbose) std::cout << "FAILED" << std::endl;
188  numberFailedTests++;
189  }
190
191  int n_potrf = 5;
192  std::vector<T> diag_a(n_potrf*n_potrf, zero);
193  for (int i=0; i<n_potrf; i++)
194  {
195  T tmp = zero;
196  sL.add( i, one, tmp );
197  diag_a[i*n_potrf + i] = tmp*tmp;
198  }
199
200  if (verbose) std::cout << "POTRF test ... ";
201  L.POTRF(char_U, n_potrf, &diag_a[0], n_potrf, &info);
202
203  if (info != 0)
204  {
205  if (verbose) std::cout << "FAILED" << std::endl;
206  numberFailedTests++;
207  }
208  else
209  {
210  for (int i=0; i<n_potrf; i++)
211  {
212  T tmp = zero;
213  sL.add( i, one, tmp );
214  if ( diag_a[i*n_potrf + i] == tmp )
215  {
216  if (verbose && i==(n_potrf-1)) std::cout << "passed!" << std::endl;
217  }
218  else
219  {
220  if (verbose) std::cout << "FAILED" << std::endl;
221  numberFailedTests++;
222  break;
223  }
224  }
225  }
226
227  if (verbose) std::cout << "POTRI test ... ";
228  std::vector<T> diag_a_trtri(diag_a); // Save a copy for TRTRI test
229
230  L.POTRI(char_U, n_potrf, &diag_a[0], n_potrf, &info);
231
232  T tmp = zero;
233  sL.multiply( 1.0/4.0, one, tmp );
234  if ( info != 0 || (diag_a[n_potrf+1] != tmp) )
235  {
236  if (verbose) std::cout << "FAILED" << std::endl;
237  numberFailedTests++;
238  }
239  else
240  if (verbose) std::cout << "passed!" << std::endl;
241
242  if (verbose) std::cout << "TRTRI test ... ";
243
244  int n_trtri = n_potrf;
245  L.TRTRI( char_U, char_N, n_trtri, &diag_a_trtri[0], n_trtri, &info );
246  for (int i=0; i<n_trtri; i++)
247  {
248  tmp = zero;
249  sL.divide( one, i+1.0, tmp );
250  if ( info != 0 )
251  {
252  numberFailedTests++;
253  break;
254  }
255  else if ( diag_a_trtri[i*n_trtri + i] == tmp )
256  {
257  if (verbose && i==(n_trtri-1)) std::cout << "passed!" << std::endl;
258  }
259  else
260  {
261  if (verbose) std::cout << "FAILED" << std::endl;
262  numberFailedTests++;
263  break;
264  }
265  }
266
267 #ifndef TEUCHOSNUMERICS_DISABLE_STEQR_TEST
268
269  if (verbose) std::cout << "STEQR test ... ";
270
271  const int n_steqr = 10;
272  std::vector<MagnitudeType> diagonal(n_steqr);
273  std::vector<MagnitudeType> subdiagonal(n_steqr-1);
274
275  for (int i=0; i < n_steqr; ++i) {
276  diagonal[i] = n_steqr - i;
277  if (i < n_steqr-1)
278  subdiagonal[i] = STM::eps() * (i+1);
279  }
280
281  std::vector<T> scalar_dummy(1,0.0);
282  std::vector<MagnitudeType> mag_dummy(4*n_steqr,0.0);
283
284  L.STEQR (char_N, n_steqr, &diagonal[0], &subdiagonal[0],
285  &scalar_dummy[0], n_steqr, &mag_dummy[0], &info);
286
287  if (info != 0)
288  {
289  if (verbose) std::cout << "STEQR: compute symmetric tridiagonal eigenvalues: "
290  << "LAPACK's _STEQR failed with info = "
291  << info;
292
293  numberFailedTests++;
294  }
295
296  MagnitudeType lambda_min = diagonal[0];
297  MagnitudeType lambda_max = diagonal[n_steqr-1];
298  MagnitudeType exp_lambda_min = STM::one();
299  MagnitudeType exp_lambda_max = STM::one()*n_steqr;
300
301  if ((fabs(lambda_min-exp_lambda_min)<1e-12) && (fabs(lambda_max-exp_lambda_max)<1e-12))
302  {
303  if (verbose) std::cout << "passed!" << std::endl;
304  }
305  else
306  {
307  if (verbose) std::cout << "FAILED" << std::endl;
308  numberFailedTests++;
309  }
310
311 #endif // TEUCHOSNUMERICS_DISABLE_STEQR_TEST
312
313  numberFailedTests += specializedLAPACK<T>::test( verbose );
314
315  return numberFailedTests;
316 }
317
318 template<class T>
319 int specializedLAPACK<T>::test(bool verbose)
320 {
321  // Create some common typedefs
322  typedef Teuchos::ScalarTraits<T> STS;
323  typedef typename STS::magnitudeType MagnitudeType;
325
326  T one = STS::one();
327  MagnitudeType m_one = STM::one();
328  T zero = STS::zero();
329
330  char char_E = 'E';
331  char char_U = 'U';
332
333  int info=0;
334  int numberFailedTests = 0;
335
337
338  if (verbose) std::cout << "LAPY2 test ... ";
339  T x = 3*one, y = 4*one;
340  T lapy = L.LAPY2(x, y);
341  if ( lapy == 5*one ) {
342  if (verbose) std::cout << "passed!" << std::endl;
343  } else {
344  if (verbose) std::cout << "FAILED ( " << lapy << " != 5 )" << std::endl;
345  numberFailedTests++;
346  }
347
348  if (verbose) std::cout << "LAMCH test ... ";
349
350  T st_eps = L.LAMCH( char_E );
351  if (verbose)
352  std::cout << "[ eps = " << st_eps << " ] passed!" << std::endl;
353
354  // Create a simple diagonal linear system
355  const int n = 4;
356  std::vector<T> Ad2_sub(n-1, zero), b2(n, one);
357  std::vector<MagnitudeType> Ad2(n, m_one);
358
359  if (verbose) std::cout << "PTTRS test ... ";
360  L.PTTRS(n, 1, &Ad2[0], &Ad2_sub[0], &b2[0], n, &info);
361  if ( !info ) {
362  if (verbose) std::cout << "passed!" << std::endl;
363  } else {
364  if (verbose) std::cout << "FAILED" << std::endl;
365  numberFailedTests++;
366  }
367
368  if (verbose) std::cout << "POCON test ... ";
369
370  std::vector<T> diag_a(n*n);
371  for (int i=0; i<n; i++)
372  {
373  diag_a[i*n + i] = one;
374  }
375  MagnitudeType rcond, anorm = m_one;
376  std::vector<T> work(3*n);
377  std::vector<int> iwork(n);
378
379  L.POCON(char_U, n, &diag_a[0], n, anorm, &rcond, &work[0], &iwork[0], &info);
380  if (info != 0 || (rcond != m_one))
381  {
382  if (verbose) std::cout << "FAILED" << std::endl;
383  numberFailedTests++;
384  }
385  else
386  if (verbose) std::cout << "passed!" << std::endl;
387
388
389  return numberFailedTests;
390 }
391
392 #ifdef HAVE_TEUCHOS_COMPLEX
393
394 template<class T>
395 int specializedLAPACK<std::complex<T> >::test( bool verbose )
396 {
397  // Create some common typedefs
399  typedef typename STS::magnitudeType MagnitudeType;
401
402  std::complex<T> one = STS::one();
403  MagnitudeType m_one = STM::one();
404  std::complex<T> zero = STS::zero();
405
406  char char_L = 'L';
407  char char_U = 'U';
408
409  int info=0;
410  int numberFailedTests = 0;
411
413
414  // Create a simple diagonal linear system
415  const int n = 4;
416  std::vector<std::complex<T> > Ad2_sub(n-1, zero), b2(n, one);
417  std::vector<MagnitudeType> Ad2(n, m_one);
418
419  if (verbose) std::cout << "PTTRS test ... ";
420  L.PTTRS(char_L, n, 1, &Ad2[0], &Ad2_sub[0], &b2[0], n, &info);
421  if ( !info ) {
422  if (verbose) std::cout << "passed!" << std::endl;
423  } else {
424  if (verbose) std::cout << "FAILED" << std::endl;
425  numberFailedTests++;
426  }
427
428  if (verbose) std::cout << "POCON test ... ";
429
430  std::vector<std::complex<T> > diag_a(n*n);
431  for (int i=0; i<n; i++)
432  {
433  diag_a[i*n + i] = one;
434  }
435  MagnitudeType rcond, anorm = m_one;
436  std::vector<std::complex<T> > work(2*n);
437  std::vector<MagnitudeType> rwork(n);
438  std::vector<int> iwork(n);
439
440  L.POCON(char_U, n, &diag_a[0], n, anorm, &rcond, &work[0], &rwork[0], &info);
441  if (info != 0 || (rcond != m_one))
442  {
443  if (verbose) std::cout << "FAILED" << std::endl;
444  numberFailedTests++;
445  }
446  else
447  if (verbose) std::cout << "passed!" << std::endl;
448
449 return numberFailedTests;
450 }
451
452 #endif
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.
Templated serial dense matrix class.
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&#39; factorization of a Hermitian/symmetric positive definite tridiagonal matrix A...
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...
This structure defines some basic traits for a scalar field type.
static int test(bool verbose)
Templated interface class to LAPACK routines.
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.
The Templated LAPACK Wrapper Class.
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.
std::string Teuchos_Version()
int size(const Comm< Ordinal > &comm)
Get the number of processes in the communicator.
ScalarType LAPY2(const ScalarType &x, const ScalarType &y) const
Computes x^2 + y^2 safely, to avoid overflow.
int main(int argc, char *argv[])
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 add(const R &a, const T &b, T &result)
Templated serial dense vector class.
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...
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 *D*L&#39; factorization of A computed by PTTRF.
ScalarType LAMCH(const char &CMACH) const
Determines machine parameters for floating point characteristics.
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.
void divide(const T &a, const R &b, T &result)
void multiply(const R &a, const T &b, T &result)
int lapackTest(bool verbose)