55 #ifdef HAVE_TEUCHOS_COMPLEX 56 #define STYPE std::complex<double> 63 #ifdef HAVE_TEUCHOS_COMPLEX 64 #define SCALARMAX STYPE(10,0) 66 #define SCALARMAX STYPE(10) 69 template<
typename TYPE>
78 template<
typename TYPE>
90 std::complex<T>
GetRandom( std::complex<T>, std::complex<T> );
102 int main(
int argc,
char* argv[])
111 if (argc>1)
if (argv[1][0]==
'-' && argv[1][1]==
'v') verbose =
true;
116 int numberFailedTests = 0;
118 std::string testName =
"", testType =
"";
120 #ifdef HAVE_TEUCHOS_COMPLEX 121 testType =
"COMPLEX";
126 if (verbose) std::cout<<std::endl<<
"********** CHECKING TEUCHOS SERIAL DENSE SOLVER - " << testType <<
"-VALUED **********"<<std::endl<<std::endl;
135 testName =
"Generating right-hand side vector using A*x, where x is a random vector:";
136 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
139 testName =
"Generating right-hand side vector using A^T*x, where x is a random vector:";
140 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
142 #ifdef HAVE_TEUCHOS_COMPLEX 145 testName =
"Generating right-hand side vector using A^H*x, where x is a random vector:";
146 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
160 returnCode = solver1.
factor();
161 testName =
"Simple solve: factor() random A:";
162 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
165 returnCode = solver1.
solve();
166 testName =
"Simple solve: solve() random A (NO_TRANS):";
168 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
174 returnCode = solver1.
solve();
175 testName =
"Simple solve: solve() random A (TRANS):";
177 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
179 #ifdef HAVE_TEUCHOS_COMPLEX 184 returnCode = solver1.
solve();
185 testName =
"Simple solve: solve() random A (CONJ_TRANS):";
187 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
191 returnCode = solver1.
invert();
192 testName =
"Simple solve: invert() random A:";
193 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
197 testName =
"Computing solution using inverted random A (NO_TRANS):";
199 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
202 testName =
"Computing solution using inverted random A (TRANS):";
204 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
206 #ifdef HAVE_TEUCHOS_COMPLEX 208 testName =
"Computing solution using inverted random A (CONJ_TRANS):";
210 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
214 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN 225 #ifdef HAVE_TEUCHOS_COMPLEX 238 returnCode = solver2.
factor();
239 testName =
"Solve with iterative refinement: factor() random A:";
240 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
243 returnCode = solver2.
solve();
244 testName =
"Solve with iterative refinement: solve() random A (NO_TRANS):";
246 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
252 returnCode = solver2.
solve();
253 testName =
"Solve with iterative refinement: solve() random A (TRANS):";
255 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
257 #ifdef HAVE_TEUCHOS_COMPLEX 262 returnCode = solver2.
solve();
263 testName =
"Solve with iterative refinement: solve() random A (CONJ_TRANS):";
265 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
279 #ifdef HAVE_TEUCHOS_COMPLEX 292 returnCode = solver3.
factor();
293 testName =
"Solve with matrix equilibration: factor() random A:";
294 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
297 returnCode = solver3.
solve();
298 testName =
"Solve with matrix equilibration: solve() random A (NO_TRANS):";
300 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
306 returnCode = solver3.
solve();
307 testName =
"Solve with matrix equilibration: solve() random A (TRANS):";
309 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
311 #ifdef HAVE_TEUCHOS_COMPLEX 316 returnCode = solver3.
solve();
317 testName =
"Solve with matrix equilibration: solve() random A (CONJ_TRANS):";
319 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
325 if(numberFailedTests > 0)
328 std::cout <<
"Number of failed tests: " << numberFailedTests << std::endl;
329 std::cout <<
"End Result: TEST FAILED" << std::endl;
333 if(numberFailedTests == 0)
334 std::cout <<
"End Result: TEST PASSED" << std::endl;
339 template<
typename TYPE>
340 int PrintTestResults(std::string testName, TYPE calculatedResult, TYPE expectedResult,
bool verbose)
343 if(calculatedResult == expectedResult)
345 if(verbose) std::cout << testName <<
" successful." << std::endl;
350 if(verbose) std::cout << testName <<
" unsuccessful." << std::endl;
356 int ReturnCodeCheck(std::string testName,
int returnCode,
int expectedResult,
bool verbose)
359 if(expectedResult == 0)
363 if(verbose) std::cout << testName <<
" test successful." << std::endl;
368 if(verbose) std::cout << testName <<
" test unsuccessful. Return code was " << returnCode <<
"." << std::endl;
376 if(verbose) std::cout << testName <<
" test successful -- failed as expected." << std::endl;
381 if(verbose) std::cout << testName <<
" test unsuccessful -- did not fail as expected. Return code was " << returnCode <<
"." << std::endl;
388 template<
typename TYPE>
395 std::complex<T>
GetRandom( std::complex<T> Low, std::complex<T> High)
397 T lowMag = Low.real();
398 T highMag = High.real();
401 return std::complex<T>( real, imag );
421 for (
int i=0; i<m; i++)
422 for (
int j=0; j<
n; j++)
433 for (
int i=0; i<
n; i++)
453 MagnitudeType temp = norm_diff;
457 if (temp > Tolerance)
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
Non-member helper functions on the templated serial, dense matrix/vector classes. ...
Templated serial dense matrix class.
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
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.
Teuchos::RCP< DVector > GetRandomVector(int n)
int CompareVectors(const SerialDenseVector< OTYPE, STYPE > &Vector1, const SerialDenseVector< OTYPE, STYPE > &Vector2, ScalarTraits< STYPE >::magnitudeType Tolerance)
int main(int argc, char *argv[])
TYPE GetRandom(TYPE, TYPE)
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.
void solveToRefinedSolution(bool flag)
Causes all solves to compute solution to best ability using iterative refinement. ...
Templated class for solving dense linear problems.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor...
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...
int invert()
Inverts the this matrix.
std::string Teuchos_Version()
int PrintTestResults(std::string, TYPE, TYPE, bool)
int ReturnCodeCheck(std::string, int, int, bool)
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).
Teuchos::RCP< DMatrix > GetRandomMatrix(int m, int n)
SerialDenseMatrix< OTYPE, STYPE > DMatrix
Templated serial dense vector class.
Defines basic traits for the scalar field type.
SerialDenseVector< OTYPE, STYPE > DVector
Smart reference counting pointer class for automatic garbage collection.
int factor()
Computes the in-place LU factorization of the matrix using the LAPACK routine _GETRF.
Reference-counted pointer class and non-member templated function implementations.
A class for solving dense linear problems.
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
This class creates and provides basic support for dense rectangular matrix of templated type...