56 #ifdef HAVE_TEUCHOS_COMPLEX 57 #define STYPE std::complex<double> 64 #ifdef HAVE_TEUCHOS_COMPLEX 65 #define SCALARMAX STYPE(10,0) 67 #define SCALARMAX STYPE(10) 70 template<
typename TYPE>
79 template<
typename TYPE>
91 std::complex<T>
GetRandom( std::complex<T>, std::complex<T> );
103 int main(
int argc,
char* argv[])
113 if (argc>1)
if (argv[1][0]==
'-' && argv[1][1]==
'v') verbose =
true;
118 int numberFailedTests = 0;
120 std::string testName =
"", testType =
"";
122 #ifdef HAVE_TEUCHOS_COMPLEX 123 testType =
"COMPLEX";
128 if (verbose) std::cout<<std::endl<<
"********** CHECKING TEUCHOS SERIAL QR SOLVER - " << testType <<
"-VALUED **********"<<std::endl<<std::endl;
138 testName =
"Generating right-hand side vector using A*x, where x is a random vector:";
139 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
141 #ifdef HAVE_TEUCHOS_COMPLEX 144 testName =
"Generating right-hand side vector using A^H*x, where x is a random vector:";
145 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
149 testName =
"Generating right-hand side vector using A^T*x, where x is a random vector:";
150 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
168 returnCode = solver1.
factor();
169 testName =
"Simple solve: factor() random A:";
170 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
171 returnCode = solver1.
formQ();
172 testName =
"Simple solve: formQ() random A:";
173 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
176 returnCode = solver1.
solve();
177 testName =
"Simple solve: solve() random A (NO_TRANS):";
179 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
180 testName =
"Simple solve: multiplyQ(TRANS) solveR (NO_TRANS) random A (NO_TRANS):";
183 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
184 for (
OTYPE i=0; i<
n; i++) {tmp(i) = bp(i);}
186 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
189 #ifdef HAVE_TEUCHOS_COMPLEX 190 testName =
"Simple solve: formQ() solve with explicit Q implicit R random A (NO_TRANS):";
193 for (
OTYPE i=0; i<
n; i++) {tmp(i) = bp(i);}
195 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
197 testName =
"Simple solve: formQ() solve with explicit Q explicit R random A (NO_TRANS):";
199 returnCode = solver1.
formR();
200 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
202 for (
OTYPE i=0; i<
n; i++) {tmp(i) = bp(i);}
204 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
207 testName =
"Simple solve: formQ() solve with explicit Q random A (NO_TRANS):";
210 for (
OTYPE i=0; i<
n; i++) {tmp(i) = bp(i);}
212 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
214 testName =
"Simple solve: formQ() solve with explicit Q explicit R random A (NO_TRANS):";
216 returnCode = solver1.
formR();
217 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
219 for (
OTYPE i=0; i<
n; i++) {tmp(i) = bp(i);}
221 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
225 #ifdef HAVE_TEUCHOS_COMPLEX 230 returnCode = solver1.
solve();
231 testName =
"Simple solve: solve() random A (CONJ_TRANS):";
234 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
235 testName =
"Simple solve: solveR (NO_TRANS) multiplyQ(TRANS) random A (CONJ_TRANS):";
238 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
239 for (
OTYPE i=0; i<
n; i++) {tmp2(i) = bp(i);}
241 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
245 testName =
"Simple solve: formQ() solve with explicit Q random A (CONJ_TRANS):";
249 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
260 returnCode = solver1.
solve();
261 testName =
"Simple solve: solve() random A (TRANS):";
264 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
266 testName =
"Simple solve: solveR multiplyQ(TRANS) (NO_TRANS) random A (TRANS):";
268 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
269 for (
OTYPE i=0; i<
n; i++) {tmp2(i) = bp(i);}
271 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
275 testName =
"Simple solve: formQ() solve with explicit Q random A (CONJ_TRANS):";
279 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
298 #ifdef HAVE_TEUCHOS_COMPLEX 315 MagnitudeType smlnum2 = smlnum/2;
333 returnCode = solver2.
factor();
334 testName =
"Solve with matrix equilibration: factor() random A:";
335 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
338 returnCode = solver2.
solve();
339 testName =
"Solve with matrix equilibration: solve() random A (NO_TRANS):";
341 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
343 #ifdef HAVE_TEUCHOS_COMPLEX 348 returnCode = solver2.
solve();
349 testName =
"Solve with matrix equilibration: solve() random A (CONJ_TRANS):";
352 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
359 returnCode = solver2.
solve();
360 testName =
"Solve with matrix equilibration: solve() random A (TRANS):";
363 numberFailedTests +=
ReturnCodeCheck(testName, returnCode, 0, verbose);
369 if(numberFailedTests > 0)
372 std::cout <<
"Number of failed tests: " << numberFailedTests << std::endl;
373 std::cout <<
"End Result: TEST FAILED" << std::endl;
377 if(numberFailedTests == 0)
378 std::cout <<
"End Result: TEST PASSED" << std::endl;
383 template<
typename TYPE>
384 int PrintTestResults(std::string testName, TYPE calculatedResult, TYPE expectedResult,
bool verbose)
387 if(calculatedResult == expectedResult)
389 if(verbose) std::cout << testName <<
" successful." << std::endl;
394 if(verbose) std::cout << testName <<
" unsuccessful." << std::endl;
400 int ReturnCodeCheck(std::string testName,
int returnCode,
int expectedResult,
bool verbose)
403 if(expectedResult == 0)
407 if(verbose) std::cout << testName <<
" test successful." << std::endl;
412 if(verbose) std::cout << testName <<
" test unsuccessful. Return code was " << returnCode <<
"." << std::endl;
420 if(verbose) std::cout << testName <<
" test successful -- failed as expected." << std::endl;
425 if(verbose) std::cout << testName <<
" test unsuccessful -- did not fail as expected. Return code was " << returnCode <<
"." << std::endl;
432 template<
typename TYPE>
439 std::complex<T>
GetRandom( std::complex<T> Low, std::complex<T> High)
441 T lowMag = Low.real();
442 T highMag = High.real();
445 return std::complex<T>( real, imag );
465 for (
int i=0; i<m; i++)
466 for (
int j=0; j<
n; j++) {
477 for (
int i=0; i<
n; i++) {
498 MagnitudeType temp = norm_diff;
502 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. ...
SerialDenseVector< OTYPE, STYPE > DVector
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
Templated serial dense matrix class.
int multiplyQ(ETransp transq, SerialDenseMatrix< OrdinalType, ScalarType > &C)
Left multiply the input matrix by the unitary matrix Q or its adjoint.
int solveR(ETransp transr, SerialDenseMatrix< OrdinalType, ScalarType > &C)
Solve input matrix on the left with the upper triangular matrix R or its adjoint. ...
int formQ()
Explicitly forms the unitary matrix Q.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
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.
TYPE GetRandom(TYPE, TYPE)
Templated class for solving dense linear problems.
This class creates and provides basic support for dense vectors of templated type as a specialization...
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
This structure defines some basic traits for a scalar field type.
A class for solving dense linear problems.
SerialDenseMatrix< OTYPE, STYPE > DMatrix
Teuchos::RCP< DVector > GetRandomVector(int n)
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
OrdinalType numRows() const
Returns the row dimension of this matrix.
int PrintTestResults(std::string, TYPE, TYPE, bool)
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Set all values in the matrix to a constant value.
Teuchos::RCP< DMatrix > GetRandomMatrix(int m, int n)
The Templated LAPACK Wrapper Class.
#define TEUCHOS_MAX(x, y)
std::string Teuchos_Version()
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.
int CompareVectors(const SerialDenseVector< OTYPE, STYPE > &Vector1, const SerialDenseVector< OTYPE, STYPE > &Vector2, ScalarTraits< STYPE >::magnitudeType Tolerance)
Smart reference counting pointer class for automatic garbage collection.
int formR()
Explicitly forms the upper triangular matrix R.
int factor()
Computes the in-place QR factorization of the matrix using the LAPACK routine _GETRF or the Eigen cla...
int ReturnCodeCheck(std::string, int, int, bool)
OrdinalType numCols() const
Returns the column dimension of this matrix.
Reference-counted pointer class and non-member templated function implementations.
int main(int argc, char *argv[])
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).
This class creates and provides basic support for dense rectangular matrix of templated type...