5 #define OUR_CHK_ERR(a) { { int epetra_err = a; \ 6 if (epetra_err != 0) { std::cerr << "Amesos ERROR " << epetra_err << ", " \ 7 << __FILE__ << ", line " << __LINE__ << std::endl; \ 8 relerror = 1.3e15; relresidual=1e15; return(1);} }\ 12 #include "Teuchos_ParameterList.hpp" 86 Teuchos::ParameterList ParamList,
97 bool AddToAllDiagonalElements = ParamList.get(
"AddZeroToDiag",
false ) ;
100 bool TrustMe = ParamList.get(
"TrustMe",
false );
102 bool MyVerbose = false ;
104 RCP<Epetra_CrsMatrix> MyMat ;
105 RCP<Epetra_CrsMatrix> MyMatWithDiag ;
113 assert ( EpetraMatrixType >= 0 && EpetraMatrixType <= 2 );
114 switch ( EpetraMatrixType ) {
122 MyMat->OptimizeStorage();
124 bool OptStorage = MyMat->StorageOptimized();
125 assert( OptStorage) ;
128 bool OptStorage = MyMat->StorageOptimized();
132 const std::string AC = AmesosClass ;
133 if ( ExpectedError == 0 ) {
134 if ( AC !=
"Amesos_Pardiso" ) {
136 ParamList, MatPtr, Rcond ) );
138 if (MyVerbose) std::cout <<
" AC = " << AC <<
" not tested in Partial Factorization" <<std::endl ;
142 if (ParamList.isParameter(
"AddToDiag") ) {
154 if ( AddToAllDiagonalElements ) {
163 double AddToDiag = ParamList.get(
"AddToDiag", 0.0 );
166 AddConstVecToDiag.
PutScalar( AddToDiag );
168 ierr = MyMatWithDiag->ExtractDiagonalCopy( Diag );
170 Diag.Update( 1.0, AddConstVecToDiag, 1.0 ) ;
171 ierr = MyMatWithDiag->ReplaceDiagonalValues( Diag );
179 if ( MyVerbose ) std::cout <<
" Partial Factorization complete " << std::endl ;
184 assert( Levels >= 1 && Levels <= 3 ) ;
186 int iam = Comm.
MyPID() ;
190 transpose?&MyMat->OperatorDomainMap():&MyMat->OperatorRangeMap() ;
192 transpose?&MyMat->OperatorRangeMap():&MyMat->OperatorDomainMap() ;
215 Teuchos::RCP<Amesos_BaseSolver> Abase ;
221 Abase = Teuchos::rcp(Afactory.
Create( AmesosClass, Problem )) ;
226 if ( Abase == Teuchos::null )
240 if ( transpose )
OUR_CHK_ERR( Abase->SetUseTranspose( transpose ) );
241 if (MyVerbose) ParamList.set(
"DebugLevel", 1 );
242 if (MyVerbose) ParamList.set(
"OutputLevel", 1 );
246 Problem.
SetLHS( &FixedLHS );
247 Problem.
SetRHS( &FixedRHS );
248 assert( OptStorage) ;
262 const int SymbolicFactorizationReturn = Abase->SymbolicFactorization( ) ;
266 if ( iam == 0 && SymbolicFactorizationReturn != ExpectedError ) {
267 std::cout <<
" SymbolicFactorization returned " << SymbolicFactorizationReturn
268 <<
" should be " << ExpectedError << std::endl ;
274 const int SymbolicFactorizationReturn = Abase->SymbolicFactorization( ) ;
282 const int NumericFactorizationReturn = Abase->NumericFactorization( ) ;
286 if ( iam == 0 && NumericFactorizationReturn != ExpectedError ) {
287 std::cout <<
" NumericFactorization returned " << NumericFactorizationReturn
288 <<
" should be " << ExpectedError << std::endl ;
294 const int NumericFactorizationReturn = Abase->NumericFactorization( ) ;
310 if ( MyMatWithDiag->MyGRID( 0 ) ) {
311 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
313 MyMatWithDiag->Multiply( transpose, xexact, cAx ) ;
316 if ( MyMatWithDiag->MyGRID( 0 ) )
317 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
330 if ( MyMatWithDiag->MyGRID( 0 ) )
332 MyMatWithDiag->Multiply( transpose, cAx, cAAx ) ;
335 if ( MyMatWithDiag->MyGRID( 0 ) )
336 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
343 if ( MyVerbose ) std::cout <<
" Compute b = A x2 = A A' A'' xexact " << std::endl ;
345 MyMatWithDiag->
Multiply( transpose, cAAx, b ) ;
365 if ( TrustMe ) sAAx = FixedLHS ;
369 OUR_CHK_ERR( Abase->SetUseTranspose( transpose ) );
377 if ( MyMat->MyGRID( 0 ) )
378 MyMat->SumIntoMyValues( 0, 1, val, ind ) ;
381 if ( TrustMe ) sAx = FixedLHS ;
398 if ( TrustMe ) x = FixedLHS ;
408 if ( MyMat->MyGRID( 0 ) ) {
409 if ( MyMat->SumIntoMyValues( 0, 1, val, ind ) ) {
410 std::cout <<
" TestOptions requires a non-zero entry in A(1,1) " << std::endl ;
423 if ( MyMatWithDiag->MyGRID( 0 ) )
425 MyMatWithDiag->Multiply( transpose, x, kAx ) ;
427 if ( MyMatWithDiag->MyGRID( 0 ) )
428 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
438 if ( MyMatWithDiag->MyGRID( 0 ) )
440 MyMatWithDiag->Multiply( transpose, kAx, kAAx ) ;
442 if ( MyMatWithDiag->MyGRID( 0 ) )
443 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
450 MyMatWithDiag->
Multiply( transpose, kAAx, bcheck ) ;
455 DomainDiff.
Update( 1.0, sAAx, -1.0, cAAx, 0.0 ) ;
456 DomainDiff.
Norm2( &norm_diff ) ;
457 sAAx.
Norm2( &norm_one ) ;
458 if (MyVerbose) std::cout << __FILE__ <<
"::" << __LINE__
459 <<
" norm( sAAx - cAAx ) / norm(sAAx ) = " 460 << norm_diff /norm_one << std::endl ;
465 DomainDiff.
Update( 1.0, sAx, -1.0, cAx, 0.0 ) ;
466 DomainDiff.
Norm2( &norm_diff ) ;
467 sAx.
Norm2( &norm_one ) ;
468 if (MyVerbose) std::cout
469 << __FILE__ <<
"::" << __LINE__
470 <<
" norm( sAx - cAx ) / norm(sAx ) = " 471 << norm_diff /norm_one << std::endl ;
474 DomainDiff.
Update( 1.0, x, -1.0, xexact, 0.0 ) ;
475 DomainDiff.
Norm2( &norm_diff ) ;
476 x.
Norm2( &norm_one ) ;
477 if (MyVerbose) std::cout
478 << __FILE__ <<
"::" << __LINE__
479 <<
" norm( x - xexact ) / norm(x) = " 480 << norm_diff /norm_one << std::endl ;
482 relerror = norm_diff / norm_one ;
484 DomainDiff.
Update( 1.0, sAx, -1.0, kAx, 0.0 ) ;
485 DomainDiff.
Norm2( &norm_diff ) ;
486 sAx.
Norm2( &norm_one ) ;
487 if (MyVerbose) std::cout
488 << __FILE__ <<
"::" << __LINE__
489 <<
" norm( sAx - kAx ) / norm(sAx ) = " 490 << norm_diff /norm_one << std::endl ;
493 DomainDiff.
Update( 1.0, sAAx, -1.0, kAAx, 0.0 ) ;
494 DomainDiff.
Norm2( &norm_diff ) ;
495 sAAx.
Norm2( &norm_one ) ;
496 if (MyVerbose) std::cout
497 << __FILE__ <<
"::" << __LINE__
498 <<
" norm( sAAx - kAAx ) / norm(sAAx ) = " 499 << norm_diff /norm_one << std::endl ;
502 RangeDiff.
Update( 1.0, bcheck, -1.0, b, 0.0 ) ;
503 RangeDiff.
Norm2( &norm_diff ) ;
504 bcheck.
Norm2( &norm_one ) ;
505 if (MyVerbose) std::cout
506 << __FILE__ <<
"::" << __LINE__
507 <<
" norm( bcheck - b ) / norm(bcheck ) = " 508 << norm_diff /norm_one << std::endl ;
510 relresidual = norm_diff / norm_one ;
513 if ( relresidual * Rcond < 1e-16 ) {
514 if (MyVerbose) std::cout <<
" Test 1 Passed " << std::endl ;
516 std::cout << __FILE__ <<
"::" << __LINE__ <<
517 " relresidual = " << relresidual <<
519 " ParamList = " << ParamList << std::endl ;
int Norm2(double *Result) const
void SetLHS(Epetra_MultiVector *X)
void SetOperator(Epetra_RowMatrix *A)
int PutScalar(double ScalarConstant)
static void SetTracebackMode(int TracebackModeValue)
int PartialFactorization(const char *AmesosClass, const Epetra_Comm &Comm, bool transpose, bool verbose, Teuchos::ParameterList ParamList, Epetra_CrsMatrix *&Amat, double Rcond)
Amesos_TestRowMatrix: a class to test Epetra_RowMatrix based codes.
virtual int MyPID() const=0
int SumIntoMyValues(int NumEntries, const double *Values, const int *Indices)
Factory for binding a third party direct solver to an Epetra_LinearProblem.
RCP< Epetra_CrsMatrix > NewMatNewMap(Epetra_CrsMatrix &In, int Diagonal, int ReindexRowMap, int ReindexColMap, int RangeMapType, int DomainMapType)
void SetRHS(Epetra_MultiVector *B)
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
static int GetTracebackMode()
const int NumericallySingularMatrixError
const int StructurallySingularMatrixError