57#ifdef HAVE_TEUCHOS_COMPLEX
58#define STYPE std::complex<double>
65#ifdef HAVE_TEUCHOS_COMPLEX
66#define SCALARMAX STYPE(10,0)
68#define SCALARMAX STYPE(10)
71template<
typename TYPE>
81template<
typename TYPE>
93std::complex<T>
GetRandom( std::complex<T>, std::complex<T> );
111int main(
int argc,
char* argv[])
117 int n=10, kl=2, ku=3;
121 if (argc>1)
if (argv[1][0]==
'-' && argv[1][1]==
'v') verbose =
true;
126 int numberFailedTests = 0;
128 std::string testName =
"", testType =
"";
130#ifdef HAVE_TEUCHOS_COMPLEX
131 testType =
"COMPLEX";
136 if (verbose) std::cout<<std::endl<<
"********** CHECKING TEUCHOS SERIAL BAND DENSE SOLVER - " << testType <<
"-VALUED **********"<<std::endl<<std::endl;
150 AB1 = Teuchos::generalToBanded( A1, kl, ku,
true );
153 C1 = Teuchos::bandedToGeneral( AB1 );
154 testName =
"Converting matrix formats: generalToBanded and bandedToGeneral random A:";
156 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
160 testName =
"Generating right-hand side vector using A*x, where x is a random vector:";
161 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
163 testName =
"Generating right-hand side vector using A^T*x, where x is a random vector:";
164 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
166#ifdef HAVE_TEUCHOS_COMPLEX
169 testName =
"Generating right-hand side vector using A^H*x, where x is a random vector:";
170 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
181 returnCode = solver1.
factor();
182 testName =
"Simple solve: factor() random A:";
183 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
186 returnCode = solver1.
solve();
187 testName =
"Simple solve: solve() random A (NO_TRANS):";
189 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
195 returnCode = solver1.
solve();
196 testName =
"Simple solve: solve() random A (TRANS):";
198 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
200#ifdef HAVE_TEUCHOS_COMPLEX
205 returnCode = solver1.
solve();
206 testName =
"Simple solve: solve() random A (CONJ_TRANS):";
208 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
212#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
223#ifdef HAVE_TEUCHOS_COMPLEX
234 AB2 = Teuchos::generalToBanded( A2, kl, ku,
true );
237 C2 = Teuchos::bandedToGeneral( AB2 );
238 testName =
"Converting matrix formats: generalToBanded and bandedToGeneral random A:";
240 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
247 returnCode = solver2.
factor();
248 testName =
"Solve with iterative refinement: factor() random A:";
249 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
252 returnCode = solver2.
solve();
253 testName =
"Solve with iterative refinement: solve() random A (NO_TRANS):";
255 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
261 returnCode = solver2.
solve();
262 testName =
"Solve with iterative refinement: solve() random A (TRANS):";
264 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
266#ifdef HAVE_TEUCHOS_COMPLEX
271 returnCode = solver2.
solve();
272 testName =
"Solve with iterative refinement: solve() random A (CONJ_TRANS):";
274 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
288#ifdef HAVE_TEUCHOS_COMPLEX
299 AB3 = Teuchos::generalToBanded( A3, kl, ku,
true );
302 C3 = Teuchos::bandedToGeneral( AB3 );
303 testName =
"Converting matrix formats: generalToBanded and bandedToGeneral random A:";
305 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
315 returnCode = solver3.
factor();
316 testName =
"Solve with matrix equilibration: factor() random A:";
317 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
320 returnCode = solver3.
solve();
321 testName =
"Solve with matrix equilibration: solve() random A (NO_TRANS):";
323 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
329 returnCode = solver3.
solve();
330 testName =
"Solve with matrix equilibration: solve() random A (TRANS):";
332 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
334#ifdef HAVE_TEUCHOS_COMPLEX
339 returnCode = solver3.
solve();
340 testName =
"Solve with matrix equilibration: solve() random A (CONJ_TRANS):";
342 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
350 returnCode = solver3.
solve();
351 testName =
"Solve with matrix equilibration: solve() without factor() random A (NO_TRANS):";
353 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
358 if(numberFailedTests > 0)
361 std::cout <<
"Number of failed tests: " << numberFailedTests << std::endl;
362 std::cout <<
"End Result: TEST FAILED" << std::endl;
366 if(numberFailedTests == 0)
367 std::cout <<
"End Result: TEST PASSED!" << std::endl;
372template<
typename TYPE>
373int PrintTestResults(std::string testName, TYPE calculatedResult, TYPE expectedResult,
bool verbose)
376 if(calculatedResult == expectedResult)
378 if(verbose) std::cout << testName <<
" successful." << std::endl;
383 if(verbose) std::cout << testName <<
" unsuccessful." << std::endl;
389int ReturnCodeCheck(std::string testName,
int returnCode,
int expectedResult,
bool verbose)
392 if(expectedResult == 0)
396 if(verbose) std::cout << testName <<
" test successful." << std::endl;
401 if(verbose) std::cout << testName <<
" test unsuccessful. Return code was " << returnCode <<
"." << std::endl;
409 if(verbose) std::cout << testName <<
" test successful -- failed as expected." << std::endl;
414 if(verbose) std::cout << testName <<
" test unsuccessful -- did not fail as expected. Return code was " << returnCode <<
"." << std::endl;
421template<
typename TYPE>
428std::complex<T>
GetRandom( std::complex<T> Low, std::complex<T> High)
430 T lowMag = Low.real();
431 T highMag = High.real();
434 return std::complex<T>( real, imag );
454 for (
int i=0; i<m; i++)
455 for (
int j=0; j<n; j++)
456 if (j>= i-kl && j<=i+ku)
467 for (
int i=0; i<n; i++)
487 MagnitudeType temp = norm_diff;
491 if (temp > Tolerance)
509 MagnitudeType norm_diff = diff.
normInf();
510 MagnitudeType norm_m2 = Matrix2.
normInf();
511 MagnitudeType temp = norm_diff;
515 if (temp > Tolerance)
Reference-counted pointer class and non-member templated function implementations.
Defines basic traits for the scalar field type.
Templated serial dense matrix class.
Non-member helper functions on the templated serial, dense matrix/vector classes.
Templated serial dense matrix class.
Templated serial dense vector class.
Smart reference counting pointer class for automatic garbage collection.
This class creates and provides basic support for banded dense matrices of templated type.
A class for representing and solving banded dense linear systems.
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).
int setMatrix(const RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > &AB)
Sets the pointer for coefficient matrix.
void factorWithEquilibration(bool flag)
int factor()
Computes the in-place LU factorization of the matrix using the LAPACK routine _GBTRF.
void solveWithTransposeFlag(Teuchos::ETransp trans)
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
void solveToRefinedSolution(bool flag)
Set whether or not to use iterative refinement to improve solutions to linear systems.
This class creates and provides basic support for dense rectangular matrix of templated type.
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
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.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
This class creates and provides basic support for dense vectors of templated type as a specialization...
Teuchos::RCP< DMatrix > GetRandomBandedMatrix(int m, int n, int kl, int ku)
TYPE GetRandom(TYPE, TYPE)
SerialDenseVector< OTYPE, STYPE > DVector
int PrintTestResults(std::string, TYPE, TYPE, bool)
SerialDenseMatrix< OTYPE, STYPE > DMatrix
int ReturnCodeCheck(std::string, int, int, bool)
int CompareVectors(const SerialDenseVector< OTYPE, STYPE > &Vector1, const SerialDenseVector< OTYPE, STYPE > &Vector2, ScalarTraits< STYPE >::magnitudeType Tolerance)
Teuchos::RCP< DVector > GetRandomVector(int n)
int CompareMatrices(const SerialDenseMatrix< OTYPE, STYPE > &Matrix1, const SerialDenseMatrix< OTYPE, STYPE > &Matrix2, ScalarTraits< STYPE >::magnitudeType Tolerance)
SerialBandDenseMatrix< OTYPE, STYPE > BDMatrix
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
std::string Teuchos_Version()
This structure defines some basic traits for a scalar field type.
static T one()
Returns representation of one for this scalar type.