58 #ifdef HAVE_TEUCHOS_COMPLEX 59 #define STYPE std::complex<double> 66 #ifdef HAVE_TEUCHOS_COMPLEX 67 #define SCALARMAX STYPE(10,0) 69 #define SCALARMAX STYPE(10) 72 template<
typename TYPE>
82 template<
typename TYPE>
94 std::complex<T>
GetRandom( std::complex<T>, std::complex<T> );
107 int main(
int argc,
char* argv[])
115 if (argc>1)
if (argv[1][0]==
'-' && argv[1][1]==
'v') verbose =
true;
120 int numberFailedTests = 0, failedComparison = 0;
122 std::string testName =
"", testType =
"";
124 #ifdef HAVE_TEUCHOS_COMPLEX 125 testType =
"COMPLEX";
130 if (verbose) std::cout<<std::endl<<
"********** CHECKING TEUCHOS SERIAL SPD DENSE SOLVER - " << testType <<
"-VALUED **********"<<std::endl<<std::endl;
138 for (
int j=0; j<
n; ++j) {
139 for (
int i=0; i<=j; ++i) {
140 (*A1full)(i,j) = (*A1)(i,j);
143 for (
int j=0; j<
n; ++j) {
144 for (
int i=j+1; i<
n; ++i) {
151 testName =
"Generating right-hand side vector using A*x, where x is a random vector:";
152 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
165 returnCode = solver1.
factor();
166 testName =
"Simple solve: factor() random A:";
167 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
170 returnCode = solver1.
solve();
171 testName =
"Simple solve: solve() random A (NO_TRANS):";
173 if(verbose && failedComparison>0) std::cout <<
"COMPARISON FAILED : ";
174 numberFailedTests += failedComparison;
175 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
178 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 186 for (
int j=0; j<
n; ++j) {
187 for (
int i=0; i<=j; ++i) {
188 (*A2full)(i,j) = (*A2)(i,j);
191 for (
int j=0; j<
n; ++j) {
192 for (
int i=j+1; i<
n; ++i) {
210 returnCode = solver2.
factor();
211 testName =
"Solve with iterative refinement: factor() random A:";
212 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
215 returnCode = solver2.
solve();
216 testName =
"Solve with iterative refinement: solve() random A (NO_TRANS):";
218 if(verbose && failedComparison>0) std::cout <<
"COMPARISON FAILED : ";
219 numberFailedTests += failedComparison;
220 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
231 for (
int j=0; j<
n; ++j) {
232 for (
int i=0; i<=j; ++i) {
233 (*A3full)(i,j) = (*A3)(i,j);
236 for (
int j=0; j<
n; ++j) {
237 for (
int i=j+1; i<
n; ++i) {
255 returnCode = solver3.
factor();
256 testName =
"Solve with matrix equilibration: factor() random A:";
257 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
260 returnCode = solver3.
solve();
261 testName =
"Solve with matrix equilibration: solve() random A (NO_TRANS):";
263 if(verbose && failedComparison>0) std::cout <<
"COMPARISON FAILED : ";
264 numberFailedTests += failedComparison;
265 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
270 if(numberFailedTests > 0)
273 std::cout <<
"Number of failed tests: " << numberFailedTests << std::endl;
274 std::cout <<
"End Result: TEST FAILED" << std::endl;
278 if(numberFailedTests == 0)
279 std::cout <<
"End Result: TEST PASSED" << std::endl;
285 template<
typename TYPE>
286 int PrintTestResults(std::string testName, TYPE calculatedResult, TYPE expectedResult,
bool verbose)
289 if(calculatedResult == expectedResult)
291 if(verbose) std::cout << testName <<
" successful." << std::endl;
296 if(verbose) std::cout << testName <<
" unsuccessful." << std::endl;
302 int ReturnCodeCheck(std::string testName,
int returnCode,
int expectedResult,
bool verbose)
305 if(expectedResult == 0)
309 if(verbose) std::cout << testName <<
" test successful." << std::endl;
314 if(verbose) std::cout << testName <<
" test unsuccessful. Return code was " << returnCode <<
"." << std::endl;
322 if(verbose) std::cout << testName <<
" test successful -- failed as expected." << std::endl;
327 if(verbose) std::cout << testName <<
" test unsuccessful -- did not fail as expected. Return code was " << returnCode <<
"." << std::endl;
334 template<
typename TYPE>
341 std::complex<T>
GetRandom( std::complex<T> Low, std::complex<T> High)
343 T lowMag = Low.real();
344 T highMag = High.real();
347 return std::complex<T>( real, imag );
381 for (
int i=0; i<
n; ++i) {
391 for (
int j=0; j<
n; ++j) {
392 for (
int i=0; i<=j; ++i) {
393 (*mat)(i,j) = (*tmp2)(i,j);
405 for (
int i=0; i<m; i++)
406 for (
int j=0; j<
n; j++)
417 for (
int i=0; i<
n; i++)
437 MagnitudeType temp = norm_diff;
441 if (temp > Tolerance)
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
void solveToRefinedSolution(bool flag)
Causes all solves to compute solution to best ability using iterative refinement. ...
SerialSymDenseMatrix< OTYPE, STYPE > SDMatrix
Teuchos::RCP< SDMatrix > GetRandomSpdMatrix(int n)
A class for constructing and using Hermitian positive definite dense matrices.
Non-member helper functions on the templated serial, dense matrix/vector classes. ...
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
Templated serial dense matrix class.
int formQ()
Explicitly forms the unitary matrix Q.
int ReturnCodeCheck(std::string, int, int, bool)
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
int main(int argc, char *argv[])
Teuchos::RCP< DVector > GetRandomVector(int n)
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
Multiply A * B and add them to this; this = beta * this + alpha*A*B.
int factor()
Computes the in-place Cholesky factorization of the matrix using the LAPACK routine DPOTRF...
This class creates and provides basic support for symmetric, positive-definite dense matrices of temp...
int CompareVectors(const SerialDenseVector< OTYPE, STYPE > &Vector1, const SerialDenseVector< OTYPE, STYPE > &Vector2, ScalarTraits< STYPE >::magnitudeType Tolerance)
Templated class for solving dense linear problems.
This class creates and provides basic support for dense vectors of templated type as a specialization...
This structure defines some basic traits for a scalar field type.
A class for solving dense linear problems.
void setUpper()
Specify that the upper triangle of the this matrix should be used.
int PrintTestResults(std::string, TYPE, TYPE, bool)
Teuchos::RCP< DMatrix > GetRandomMatrix(int m, int n)
TYPE GetRandom(TYPE, TYPE)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
int setMatrix(const RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > &A_in)
Sets the pointers for coefficient matrix.
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
int setVectors(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &X, const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &B)
Sets the pointers for left and right hand side vector(s).
std::string Teuchos_Version()
SerialDenseVector< OTYPE, STYPE > DVector
Templated serial, dense, symmetric matrix class.
Templated class for constructing and using Hermitian positive definite dense matrices.
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor...
Templated serial dense vector class.
Defines basic traits for the scalar field type.
Smart reference counting pointer class for automatic garbage collection.
SerialDenseMatrix< OTYPE, STYPE > DMatrix
int factor()
Computes the in-place QR factorization of the matrix using the LAPACK routine _GETRF or the Eigen cla...
Reference-counted pointer class and non-member templated function implementations.
This class creates and provides basic support for dense rectangular matrix of templated type...