Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
46 
47 int main(int argc, char* argv[])
48 {
49  int numberFailedTests = 0;
50  bool verbose = 0;
51  if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
52 
53  if (verbose)
54  std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
55 
56  using std::fabs;
57 
58  // Define some common characters
59  int info=0;
60  char char_N = 'N';
61  char char_U = 'U';
62 
63  // Create some common typedefs
65  typedef STS::magnitudeType MagnitudeType;
67 
70 
71  const int n_gesv = 4;
72  std::vector<double> Ad(n_gesv*n_gesv,0.0);
73  std::vector<double> bd(n_gesv,0.0);
74  std::vector<float> Af(n_gesv*n_gesv,0.0);
75  std::vector<float> bf(n_gesv,0.0);
76  int IPIV[n_gesv];
77 
78  Ad[0] = 1; Ad[2] = 1; Ad[5] = 1; Ad[8] = 2; Ad[9] = 1; Ad[10] = 1; Ad[14] = 2; Ad[15] = 2;
79  bd[1] = 2; bd[2] = 1; bd[3] = 2;
80  Af[0] = 1; Af[2] = 1; Af[5] = 1; Af[8] = 2; Af[9] = 1; Af[10] = 1; Af[14] = 2; Af[15] = 2;
81  bf[1] = 2; bf[2] = 1; bf[3] = 2;
82 
83  if (verbose) std::cout << "GESV test ... ";
84  L.GESV(n_gesv, 1, &Ad[0], n_gesv, IPIV, &bd[0], n_gesv, &info);
85  M.GESV(n_gesv, 1, &Af[0], n_gesv, IPIV, &bf[0], n_gesv, &info);
86  for(int i = 0; i < 4; i++)
87  {
88  if (bd[i] == bf[i]) {
89  if (verbose && i==3) std::cout << "passed!" << std::endl;
90  } else {
91  if (verbose) std::cout << "FAILED" << std::endl;
92  numberFailedTests++;
93  break;
94  }
95  }
96 
97  if (verbose) std::cout << "LAPY2 test ... ";
98  float fx = 3, fy = 4;
99  float flapy = M.LAPY2(fx, fy);
100  double dx = 3, dy = 4;
101  double dlapy = L.LAPY2(dx, dy);
102  if ( dlapy == flapy && dlapy == 5.0 && flapy == 5.0f ) {
103  if (verbose) std::cout << "passed!" << std::endl;
104  } else {
105  if (verbose) std::cout << "FAILED (" << dlapy << " != " << flapy << ")" << std::endl;
106  numberFailedTests++;
107  }
108 
109  if (verbose) std::cout << "LAMCH test ... ";
110 
111  char char_E = 'E';
112  double d_eps = L.LAMCH( char_E );
113  float f_eps = M.LAMCH( char_E );
114  if (verbose)
115  std::cout << "[ Double-precision eps = " << d_eps << ", single-precision eps = " << f_eps << " ] passed!" << std::endl;
116 
117  if (verbose) std::cout << "POTRF test ... ";
118 
119  int n_potrf = 5;
120  std::vector<double> diag_a(n_potrf*n_potrf, 0.0);
121  for (int i=0; i<n_potrf; i++)
122  diag_a[i*n_potrf + i] = (i+1)*(i+1);
123  L.POTRF(char_U, n_potrf, &diag_a[0], n_potrf, &info);
124 
125  if (info != 0)
126  {
127  if (verbose) std::cout << "FAILED" << std::endl;
128  numberFailedTests++;
129  }
130  else
131  {
132  for (int i=0; i<n_potrf; i++)
133  {
134  if ( diag_a[i*n_potrf + i] == (i+1) )
135  {
136  if (verbose && i==(n_potrf-1)) std::cout << "passed!" << std::endl;
137  }
138  else
139  {
140  if (verbose) std::cout << "FAILED" << std::endl;
141  numberFailedTests++;
142  break;
143  }
144  }
145  }
146 
147  if (verbose) std::cout << "POCON test ... ";
148 
149  double anorm = (n_potrf*n_potrf), rcond;
150  std::vector<double> work(3*n_potrf);
151  std::vector<int> iwork(n_potrf);
152 
153  L.POCON(char_U, n_potrf, &diag_a[0], n_potrf, anorm, &rcond, &work[0], &iwork[0], &info);
154  if (info != 0 || (rcond != 1.0/anorm))
155  {
156  numberFailedTests++;
157  if (verbose) std::cout << "FAILED" << std::endl;
158  }
159  else
160  {
161  if (verbose) std::cout << "passed!" << std::endl;
162  }
163 
164  if (verbose) std::cout << "POTRI test ... ";
165  std::vector<double> diag_a_trtri(diag_a); // Save a copy for TRTRI test
166 
167  L.POTRI(char_U, n_potrf, &diag_a[0], n_potrf, &info);
168 
169  if (info != 0 || (diag_a[n_potrf+1] != 1.0/4.0))
170  {
171  numberFailedTests++;
172  if (verbose) std::cout << "FAILED" << std::endl;
173  }
174  else
175  {
176  if (verbose) std::cout << "passed!" << std::endl;
177  }
178 
179  if (verbose) std::cout << "TRTRI test ... ";
180 
181  int n_trtri = n_potrf;
182  L.TRTRI( char_U, char_N, n_trtri, &diag_a_trtri[0], n_trtri, &info );
183  for (int i=0; i<n_trtri; i++)
184  {
185  if ( diag_a_trtri[i*n_trtri + i] == 1.0/(i+1) )
186  {
187  if (verbose && i==(n_trtri-1)) std::cout << "passed!" << std::endl;
188  }
189  else
190  {
191  if (verbose) std::cout << "FAILED" << std::endl;
192  numberFailedTests++;
193  break;
194  }
195  }
196 
197 
198 #if ! (defined(__INTEL_COMPILER) && defined(_WIN32) )
199 
200  // Check ILAENV with similarity transformation routine: dsytrd
201  // NOTE: Do not need to put floating point specifier [s,d,c,z] before routine name,
202  // this is handled through templating.
203  if (verbose) std::cout << "ILAENV test ... ";
204  int n1 = 100;
205  int size = L.ILAENV(1, "sytrd", "u", n1);
206  if (size > 0) {
207  if (verbose) std::cout << "passed!" << std::endl;
208  } else {
209  if (verbose) std::cout << "FAILED!" << std::endl;
210  numberFailedTests++;
211  }
212 
213 #endif
214 
215  if (verbose) std::cout << "STEQR test ... ";
216 
217 #ifndef TEUCHOSNUMERICS_DISABLE_STEQR_TEST
218 
219  const int n_steqr = 10;
220  std::vector<MagnitudeType> diagonal(n_steqr);
221  std::vector<MagnitudeType> subdiagonal(n_steqr-1);
222 
223  for (int i=0; i < n_steqr; ++i) {
224  diagonal[i] = n_steqr - i;
225  if (i < n_steqr-1)
226  subdiagonal[i] = STM::eps() * (i+1);
227  }
228 
229  std::vector<double> scalar_dummy(1,0.0);
230  std::vector<MagnitudeType> mag_dummy(4*n_steqr,0.0);
231 
232  L.STEQR (char_N, n_steqr, &diagonal[0], &subdiagonal[0],
233  &scalar_dummy[0], n_steqr, &mag_dummy[0], &info);
234 
235  if (info != 0)
236  {
237  if (verbose) std::cout << "STEQR: compute symmetric tridiagonal eigenvalues: "
238  << "LAPACK's _STEQR failed with info = "
239  << info;
240 
241  numberFailedTests++;
242  }
243 
244  MagnitudeType lambda_min = diagonal[0];
245  MagnitudeType lambda_max = diagonal[n_steqr-1];
246  MagnitudeType exp_lambda_min = STM::one();
247  MagnitudeType exp_lambda_max = STM::one()*n_steqr;
248 
249  if ((fabs(lambda_min-exp_lambda_min)<1e-12) && (fabs(lambda_max-exp_lambda_max)<1e-12))
250  {
251  if (verbose) std::cout << "passed!" << std::endl;
252  }
253  else
254  {
255  if (verbose) std::cout << "FAILED" << std::endl;
256  numberFailedTests++;
257  }
258 
259 #else // TEUCHOSNUMERICS_DISABLE_STEQR_TEST
260 
261  if (verbose) std::cout << "SKIPPED!\n";
262 
263 #endif // TEUCHOSNUMERICS_DISABLE_STEQR_TEST
264 
265  if(numberFailedTests > 0)
266  {
267  if (verbose) {
268  std::cout << "Number of failed tests: " << numberFailedTests << std::endl;
269  std::cout << "End Result: TEST FAILED" << std::endl;
270  return -1;
271  }
272  }
273  if(numberFailedTests==0)
274  std::cout << "End Result: TEST PASSED" << std::endl;
275  return 0;
276 
277 }
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 TRTRI(const char &UPLO, const char &DIAG, const OrdinalType &n, const ScalarType *A, const OrdinalType &lda, OrdinalType *info) const
Computes the inverse of an upper or lower triangular matrix A.
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 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.
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.
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 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 LAMCH(const char &CMACH) const
Determines machine parameters for floating point characteristics.