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; \
8relerror = 1.3e15; relresidual=1e15; return(1);} }\
11#include "Epetra_Comm.h"
12#include "Teuchos_ParameterList.hpp"
14#include "Epetra_CrsMatrix.h"
15#include "Epetra_Map.h"
16#include "Epetra_Vector.h"
17#include "Epetra_LinearProblem.h"
83 const Epetra_Comm &Comm,
86 Teuchos::ParameterList ParamList,
87 Epetra_CrsMatrix *& InMat,
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 ;
107 MyMat = rcp(
new Epetra_CrsMatrix( *InMat ) );
111 Epetra_RowMatrix* MyRowMat = 0;
113 assert ( EpetraMatrixType >= 0 && EpetraMatrixType <= 2 );
114 switch ( EpetraMatrixType ) {
122 MyMat->OptimizeStorage();
127 MyMat->StorageOptimized();
128 assert( OptStorage) ;
134 MyMat->StorageOptimized();
136 Epetra_CrsMatrix* MatPtr = &*MyMat ;
138 const std::string AC = AmesosClass ;
139 if ( ExpectedError == 0 ) {
140 if ( AC !=
"Amesos_Pardiso" ) {
142 ParamList, MatPtr, Rcond ) );
144 if (MyVerbose) std::cout <<
" AC = " << AC <<
" not tested in Partial Factorization" <<std::endl ;
148 if (ParamList.isParameter(
"AddToDiag") ) {
149 int oldtracebackmode = InMat->GetTracebackMode( ) ;
160 if ( AddToAllDiagonalElements ) {
163 InMat->SetTracebackMode( EPETRA_MIN(oldtracebackmode,1) ) ;
169 double AddToDiag = ParamList.get(
"AddToDiag", 0.0 );
170 Epetra_Vector Diag( MyMatWithDiag->RowMap() );
171 Epetra_Vector AddConstVecToDiag( MyMatWithDiag->RowMap() );
172 AddConstVecToDiag.PutScalar( AddToDiag );
177 MyMatWithDiag->ExtractDiagonalCopy( Diag );
179 Diag.Update( 1.0, AddConstVecToDiag, 1.0 ) ;
183 MyMatWithDiag->ReplaceDiagonalValues( Diag );
186 InMat->SetTracebackMode( oldtracebackmode ) ;
188 MyMatWithDiag = rcp(
new Epetra_CrsMatrix( *InMat ) );
191 if ( MyVerbose ) std::cout <<
" Partial Factorization complete " << std::endl ;
196 assert( Levels >= 1 && Levels <= 3 ) ;
198 int iam = Comm.MyPID() ;
201 const Epetra_Map *RangeMap =
202 transpose?&MyMat->OperatorDomainMap():&MyMat->OperatorRangeMap() ;
203 const Epetra_Map *DomainMap =
204 transpose?&MyMat->OperatorRangeMap():&MyMat->OperatorDomainMap() ;
206 Epetra_Vector xexact(*DomainMap);
207 Epetra_Vector x(*DomainMap);
209 Epetra_Vector cAx(*DomainMap);
210 Epetra_Vector sAx(*DomainMap);
211 Epetra_Vector kAx(*DomainMap);
213 Epetra_Vector cAAx(*DomainMap);
214 Epetra_Vector sAAx(*DomainMap);
215 Epetra_Vector kAAx(*DomainMap);
217 Epetra_Vector FixedLHS(*DomainMap);
218 Epetra_Vector FixedRHS(*RangeMap);
220 Epetra_Vector b(*RangeMap);
221 Epetra_Vector bcheck(*RangeMap);
223 Epetra_Vector DomainDiff(*DomainMap);
224 Epetra_Vector RangeDiff(*RangeMap);
226 Epetra_LinearProblem Problem;
227 Teuchos::RCP<Amesos_BaseSolver> Abase ;
233 Abase = Teuchos::rcp(Afactory.
Create( AmesosClass, Problem )) ;
238 if ( Abase == Teuchos::null )
245 Problem.SetOperator( MyRowMat );
252 if ( transpose )
OUR_CHK_ERR( Abase->SetUseTranspose( transpose ) );
253 if (MyVerbose) ParamList.set(
"DebugLevel", 1 );
254 if (MyVerbose) ParamList.set(
"OutputLevel", 1 );
258 Problem.SetLHS( &FixedLHS );
259 Problem.SetRHS( &FixedRHS );
260 assert( OptStorage) ;
270 Epetra_CrsMatrix* ECM =
dynamic_cast<Epetra_CrsMatrix*
>(MyRowMat) ;
271 int oldtracebackmode = ECM->GetTracebackMode( ) ;
272 ECM->SetTracebackMode(0);
274 const int SymbolicFactorizationReturn = Abase->SymbolicFactorization( ) ;
275 ECM->SetTracebackMode(oldtracebackmode);
278 if ( iam == 0 && SymbolicFactorizationReturn != ExpectedError ) {
279 std::cout <<
" SymbolicFactorization returned " << SymbolicFactorizationReturn
280 <<
" should be " << ExpectedError << std::endl ;
286 const int SymbolicFactorizationReturn = Abase->SymbolicFactorization( ) ;
290 Epetra_CrsMatrix* ECM =
dynamic_cast<Epetra_CrsMatrix*
>(MyRowMat) ;
291 int oldtracebackmode = ECM->GetTracebackMode( ) ;
292 ECM->SetTracebackMode(0);
294 const int NumericFactorizationReturn = Abase->NumericFactorization( ) ;
295 ECM->SetTracebackMode(oldtracebackmode);
298 if ( iam == 0 && NumericFactorizationReturn != ExpectedError ) {
299 std::cout <<
" NumericFactorization returned " << NumericFactorizationReturn
300 <<
" should be " << ExpectedError << std::endl ;
306 const int NumericFactorizationReturn = Abase->NumericFactorization( ) ;
313 xexact.PutScalar(1.0);
322 if ( MyMatWithDiag->MyGRID( 0 ) ) {
323 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
325 MyMatWithDiag->Multiply( transpose, xexact, cAx ) ;
328 if ( MyMatWithDiag->MyGRID( 0 ) )
329 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
342 if ( MyMatWithDiag->MyGRID( 0 ) )
343 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
344 MyMatWithDiag->Multiply( transpose, cAx, cAAx ) ;
347 if ( MyMatWithDiag->MyGRID( 0 ) )
348 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
355 if ( MyVerbose ) std::cout <<
" Compute b = A x2 = A A' A'' xexact " << std::endl ;
357 MyMatWithDiag->Multiply( transpose, cAAx, b ) ;
368 Problem.SetLHS( &sAAx );
369 Problem.SetRHS( &b );
377 if ( TrustMe ) sAAx = FixedLHS ;
381 OUR_CHK_ERR( Abase->SetUseTranspose( transpose ) );
385 Problem.SetLHS( &sAx );
386 Problem.SetRHS( &sAAx );
389 if ( MyMat->MyGRID( 0 ) )
390 MyMat->SumIntoMyValues( 0, 1, val, ind ) ;
393 if ( TrustMe ) sAx = FixedLHS ;
406 Problem.SetLHS( &x );
407 Problem.SetRHS( &sAx );
410 if ( TrustMe ) x = FixedLHS ;
420 if ( MyMat->MyGRID( 0 ) ) {
421 if ( MyMat->SumIntoMyValues( 0, 1, val, ind ) ) {
422 std::cout <<
" TestOptions requires a non-zero entry in A(1,1) " << std::endl ;
435 if ( MyMatWithDiag->MyGRID( 0 ) )
436 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
437 MyMatWithDiag->Multiply( transpose, x, kAx ) ;
439 if ( MyMatWithDiag->MyGRID( 0 ) )
440 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
450 if ( MyMatWithDiag->MyGRID( 0 ) )
451 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
452 MyMatWithDiag->Multiply( transpose, kAx, kAAx ) ;
454 if ( MyMatWithDiag->MyGRID( 0 ) )
455 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
462 MyMatWithDiag->Multiply( transpose, kAAx, bcheck ) ;
467 DomainDiff.Update( 1.0, sAAx, -1.0, cAAx, 0.0 ) ;
468 DomainDiff.Norm2( &norm_diff ) ;
469 sAAx.Norm2( &norm_one ) ;
470 if (MyVerbose) std::cout << __FILE__ <<
"::" << __LINE__
471 <<
" norm( sAAx - cAAx ) / norm(sAAx ) = "
472 << norm_diff /norm_one << std::endl ;
477 DomainDiff.Update( 1.0, sAx, -1.0, cAx, 0.0 ) ;
478 DomainDiff.Norm2( &norm_diff ) ;
479 sAx.Norm2( &norm_one ) ;
480 if (MyVerbose) std::cout
481 << __FILE__ <<
"::" << __LINE__
482 <<
" norm( sAx - cAx ) / norm(sAx ) = "
483 << norm_diff /norm_one << std::endl ;
486 DomainDiff.Update( 1.0, x, -1.0, xexact, 0.0 ) ;
487 DomainDiff.Norm2( &norm_diff ) ;
488 x.Norm2( &norm_one ) ;
489 if (MyVerbose) std::cout
490 << __FILE__ <<
"::" << __LINE__
491 <<
" norm( x - xexact ) / norm(x) = "
492 << norm_diff /norm_one << std::endl ;
494 relerror = norm_diff / norm_one ;
496 DomainDiff.Update( 1.0, sAx, -1.0, kAx, 0.0 ) ;
497 DomainDiff.Norm2( &norm_diff ) ;
498 sAx.Norm2( &norm_one ) ;
499 if (MyVerbose) std::cout
500 << __FILE__ <<
"::" << __LINE__
501 <<
" norm( sAx - kAx ) / norm(sAx ) = "
502 << norm_diff /norm_one << std::endl ;
505 DomainDiff.Update( 1.0, sAAx, -1.0, kAAx, 0.0 ) ;
506 DomainDiff.Norm2( &norm_diff ) ;
507 sAAx.Norm2( &norm_one ) ;
508 if (MyVerbose) std::cout
509 << __FILE__ <<
"::" << __LINE__
510 <<
" norm( sAAx - kAAx ) / norm(sAAx ) = "
511 << norm_diff /norm_one << std::endl ;
514 RangeDiff.Update( 1.0, bcheck, -1.0, b, 0.0 ) ;
515 RangeDiff.Norm2( &norm_diff ) ;
516 bcheck.Norm2( &norm_one ) ;
517 if (MyVerbose) std::cout
518 << __FILE__ <<
"::" << __LINE__
519 <<
" norm( bcheck - b ) / norm(bcheck ) = "
520 << norm_diff /norm_one << std::endl ;
522 relresidual = norm_diff / norm_one ;
525 if ( relresidual * Rcond < 1e-16 ) {
526 if (MyVerbose) std::cout <<
" Test 1 Passed " << std::endl ;
528 std::cout << __FILE__ <<
"::" << __LINE__ <<
529 " relresidual = " << relresidual <<
531 " ParamList = " << ParamList << std::endl ;
const int NumericallySingularMatrixError
const int StructurallySingularMatrixError
RCP< Epetra_CrsMatrix > NewMatNewMap(Epetra_CrsMatrix &In, int Diagonal, int ReindexRowMap, int ReindexColMap, int RangeMapType, int DomainMapType)
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.
Factory for binding a third party direct solver to an Epetra_LinearProblem.
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.