57 static int test(
bool verbose);
61 void add(
const R& a,
const T& b, T& result ){ result = a + b; }
64 void multiply(
const R& a,
const T& b, T& result ){ result = a * b; }
67 void divide(
const T& a,
const R& b, T& result ){ result = a / b; }
70 #ifdef HAVE_TEUCHOS_COMPLEX
77 static int test(
bool verbose);
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; }
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; }
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; }
96 int main(
int argc,
char* argv[])
98 int numberFailedTests = 0;
100 if (argc>1)
if (argv[1][0]==
'-' && argv[1][1]==
'v') verbose =
true;
107 #ifdef HAVE_TEUCHOS_INST_FLOAT
109 std::cout << std::endl <<
"LAPACK test for float" << std::endl;
110 numberFailedTests += lapackTest<float>(verbose);
113 std::cout << std::endl <<
"LAPACK test for double" << std::endl;
114 numberFailedTests += lapackTest<double>(verbose);
116 #ifdef HAVE_TEUCHOS_COMPLEX
117 #ifdef HAVE_TEUCHOS_INST_COMPLEX_FLOAT
119 std::cout << std::endl <<
"LAPACK test for std::complex<float>" << std::endl;
120 numberFailedTests += lapackTest<std::complex<float> >(verbose);
122 #ifdef HAVE_TEUCHOS_INST_COMPLEX_DOUBLE
124 std::cout << std::endl <<
"LAPACK test for std::complex<double>" << std::endl;
125 numberFailedTests += lapackTest<std::complex<double> >(verbose);
129 if(numberFailedTests > 0)
132 std::cout <<
"Number of failed tests: " << numberFailedTests << std::endl;
133 std::cout <<
"End Result: TEST FAILED" << std::endl;
137 if(numberFailedTests==0)
138 std::cout <<
"End Result: TEST PASSED" << std::endl;
144 template <
typename T>
147 int numberFailedTests = 0;
157 typedef typename STS::magnitudeType MagnitudeType;
161 MagnitudeType m_one = STM::one();
162 T zero = STS::zero();
167 const int n_gesv = 4;
168 std::vector<T> Ad(n_gesv*n_gesv,zero);
169 std::vector<T> bd(n_gesv,zero);
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;
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);
178 if (verbose) std::cout <<
"passed!" << std::endl;
180 if (verbose) std::cout <<
"FAILED" << std::endl;
184 if (verbose) std::cout <<
"GESV test ... ";
185 L.
GESV(n_gesv, 1, &Ad[0], n_gesv, IPIV, &bd[0], n_gesv, &info);
187 if (verbose) std::cout <<
"passed!" << std::endl;
189 if (verbose) std::cout <<
"FAILED" << std::endl;
193 #if ! (defined(__INTEL_COMPILER) && defined(_WIN32) )
198 if (verbose) std::cout <<
"ILAENV test ... ";
202 if (verbose) std::cout <<
"passed!" << std::endl;
204 if (verbose) std::cout <<
"FAILED!" << std::endl;
211 std::vector<T> Ad2_sub(n_gesv-1, zero), b2(n_gesv, one);
212 std::vector<MagnitudeType> Ad2(n_gesv, m_one);
214 if (verbose) std::cout <<
"PTTRF test ... ";
215 L.
PTTRF(n_gesv, &Ad2[0], &Ad2_sub[0], &info);
217 if (verbose) std::cout <<
"passed!" << std::endl;
219 if (verbose) std::cout <<
"FAILED" << std::endl;
224 std::vector<T> diag_a(n_potrf*n_potrf, zero);
225 for (
int i=0; i<n_potrf; i++)
228 sL.
add( i, one, tmp );
229 diag_a[i*n_potrf + i] = tmp*tmp;
232 if (verbose) std::cout <<
"POTRF test ... ";
233 L.
POTRF(char_U, n_potrf, &diag_a[0], n_potrf, &info);
237 if (verbose) std::cout <<
"FAILED" << std::endl;
242 for (
int i=0; i<n_potrf; i++)
245 sL.
add( i, one, tmp );
246 if ( diag_a[i*n_potrf + i] == tmp )
248 if (verbose && i==(n_potrf-1)) std::cout <<
"passed!" << std::endl;
252 if (verbose) std::cout <<
"FAILED" << std::endl;
259 if (verbose) std::cout <<
"POTRI test ... ";
260 std::vector<T> diag_a_trtri(diag_a);
262 L.
POTRI(char_U, n_potrf, &diag_a[0], n_potrf, &info);
266 if ( info != 0 || (diag_a[n_potrf+1] != tmp) )
268 if (verbose) std::cout <<
"FAILED" << std::endl;
272 if (verbose) std::cout <<
"passed!" << std::endl;
274 if (verbose) std::cout <<
"TRTRI test ... ";
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++)
281 sL.
divide( one, i+1.0, tmp );
287 else if ( diag_a_trtri[i*n_trtri + i] == tmp )
289 if (verbose && i==(n_trtri-1)) std::cout <<
"passed!" << std::endl;
293 if (verbose) std::cout <<
"FAILED" << std::endl;
299 #ifndef TEUCHOSNUMERICS_DISABLE_STEQR_TEST
301 if (verbose) std::cout <<
"STEQR test ... ";
303 const int n_steqr = 10;
304 std::vector<MagnitudeType> diagonal(n_steqr);
305 std::vector<MagnitudeType> subdiagonal(n_steqr-1);
307 for (
int i=0; i < n_steqr; ++i) {
308 diagonal[i] = n_steqr - i;
310 subdiagonal[i] = STM::eps() * (i+1);
313 std::vector<T> scalar_dummy(1,0.0);
314 std::vector<MagnitudeType> mag_dummy(4*n_steqr,0.0);
316 L.
STEQR (char_N, n_steqr, &diagonal[0], &subdiagonal[0],
317 &scalar_dummy[0], n_steqr, &mag_dummy[0], &info);
321 if (verbose) std::cout <<
"STEQR: compute symmetric tridiagonal eigenvalues: "
322 <<
"LAPACK's _STEQR failed with info = "
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;
333 if ((fabs(lambda_min-exp_lambda_min)<1e-12) && (fabs(lambda_max-exp_lambda_max)<1e-12))
335 if (verbose) std::cout <<
"passed!" << std::endl;
339 if (verbose) std::cout <<
"FAILED" << std::endl;
343 #endif // TEUCHOSNUMERICS_DISABLE_STEQR_TEST
347 return numberFailedTests;
355 typedef typename STS::magnitudeType MagnitudeType;
359 MagnitudeType m_one = STM::one();
360 T zero = STS::zero();
366 int numberFailedTests = 0;
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;
376 if (verbose) std::cout <<
"FAILED ( " << lapy <<
" != 5 )" << std::endl;
380 if (verbose) std::cout <<
"LAMCH test ... ";
382 T st_eps = L.
LAMCH( char_E );
384 std::cout <<
"[ eps = " << st_eps <<
" ] passed!" << std::endl;
388 std::vector<T> Ad2_sub(n-1, zero), b2(n, one);
389 std::vector<MagnitudeType> Ad2(n, m_one);
391 if (verbose) std::cout <<
"PTTRS test ... ";
392 L.
PTTRS(n, 1, &Ad2[0], &Ad2_sub[0], &b2[0], n, &info);
394 if (verbose) std::cout <<
"passed!" << std::endl;
396 if (verbose) std::cout <<
"FAILED" << std::endl;
400 if (verbose) std::cout <<
"POCON test ... ";
402 std::vector<T> diag_a(n*n);
403 for (
int i=0; i<
n; i++)
405 diag_a[i*n + i] = one;
407 MagnitudeType rcond, anorm = m_one;
408 std::vector<T> work(3*n);
409 std::vector<int> iwork(n);
411 L.
POCON(char_U, n, &diag_a[0], n, anorm, &rcond, &work[0], &iwork[0], &info);
412 if (info != 0 || (rcond != m_one))
414 if (verbose) std::cout <<
"FAILED" << std::endl;
418 if (verbose) std::cout <<
"passed!" << std::endl;
421 return numberFailedTests;
424 #ifdef HAVE_TEUCHOS_COMPLEX
431 typedef typename STS::magnitudeType MagnitudeType;
434 std::complex<T> one = STS::one();
435 MagnitudeType m_one = STM::one();
436 std::complex<T> zero = STS::zero();
442 int numberFailedTests = 0;
448 std::vector<std::complex<T> > Ad2_sub(n-1, zero), b2(n, one);
449 std::vector<MagnitudeType> Ad2(n, m_one);
451 if (verbose) std::cout <<
"PTTRS test ... ";
452 L.
PTTRS(char_L, n, 1, &Ad2[0], &Ad2_sub[0], &b2[0], n, &info);
454 if (verbose) std::cout <<
"passed!" << std::endl;
456 if (verbose) std::cout <<
"FAILED" << std::endl;
460 if (verbose) std::cout <<
"POCON test ... ";
462 std::vector<std::complex<T> > diag_a(n*n);
463 for (
int i=0; i<
n; i++)
465 diag_a[i*n + i] = one;
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);
472 L.
POCON(char_U, n, &diag_a[0], n, anorm, &rcond, &work[0], &rwork[0], &info);
473 if (info != 0 || (rcond != m_one))
475 if (verbose) std::cout <<
"FAILED" << std::endl;
479 if (verbose) std::cout <<
"passed!" << std::endl;
481 return numberFailedTests;
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' 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' 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)