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);} }\ 11 #include "Epetra_Comm.h" 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,
87 Epetra_CrsMatrix *& InMat,
97 bool AddToAllDiagonalElements = ParamList.
get(
"AddZeroToDiag",
false ) ;
100 bool TrustMe = ParamList.
get(
"TrustMe",
false );
102 bool MyVerbose = false ;
107 MyMat =
rcp(
new Epetra_CrsMatrix( *InMat ) );
111 Epetra_RowMatrix* MyRowMat = 0;
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();
130 Epetra_CrsMatrix* MatPtr = &*MyMat ;
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 ;
143 int oldtracebackmode = InMat->GetTracebackMode( ) ;
154 if ( AddToAllDiagonalElements ) {
157 InMat->SetTracebackMode( EPETRA_MIN(oldtracebackmode,1) ) ;
163 double AddToDiag = ParamList.
get(
"AddToDiag", 0.0 );
164 Epetra_Vector Diag( MyMatWithDiag->RowMap() );
165 Epetra_Vector AddConstVecToDiag( MyMatWithDiag->RowMap() );
166 AddConstVecToDiag.PutScalar( AddToDiag );
168 ierr = MyMatWithDiag->ExtractDiagonalCopy( Diag );
170 Diag.Update( 1.0, AddConstVecToDiag, 1.0 ) ;
171 ierr = MyMatWithDiag->ReplaceDiagonalValues( Diag );
174 InMat->SetTracebackMode( oldtracebackmode ) ;
176 MyMatWithDiag =
rcp(
new Epetra_CrsMatrix( *InMat ) );
179 if ( MyVerbose ) std::cout <<
" Partial Factorization complete " << std::endl ;
184 assert( Levels >= 1 && Levels <= 3 ) ;
186 int iam =
Comm.MyPID() ;
189 const Epetra_Map *RangeMap =
190 transpose?&MyMat->OperatorDomainMap():&MyMat->OperatorRangeMap() ;
191 const Epetra_Map *DomainMap =
192 transpose?&MyMat->OperatorRangeMap():&MyMat->OperatorDomainMap() ;
194 Epetra_Vector xexact(*DomainMap);
195 Epetra_Vector x(*DomainMap);
197 Epetra_Vector cAx(*DomainMap);
198 Epetra_Vector sAx(*DomainMap);
199 Epetra_Vector kAx(*DomainMap);
201 Epetra_Vector cAAx(*DomainMap);
202 Epetra_Vector sAAx(*DomainMap);
203 Epetra_Vector kAAx(*DomainMap);
205 Epetra_Vector FixedLHS(*DomainMap);
206 Epetra_Vector FixedRHS(*RangeMap);
208 Epetra_Vector b(*RangeMap);
209 Epetra_Vector bcheck(*RangeMap);
211 Epetra_Vector DomainDiff(*DomainMap);
212 Epetra_Vector RangeDiff(*RangeMap);
214 Epetra_LinearProblem Problem;
226 if ( Abase == Teuchos::null )
233 Problem.SetOperator( MyRowMat );
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) ;
258 Epetra_CrsMatrix* ECM =
dynamic_cast<Epetra_CrsMatrix*
>(MyRowMat) ;
259 int oldtracebackmode = ECM->GetTracebackMode( ) ;
260 ECM->SetTracebackMode(0);
263 ECM->SetTracebackMode(oldtracebackmode);
266 if ( iam == 0 && SymbolicFactorizationReturn != ExpectedError ) {
267 std::cout <<
" SymbolicFactorization returned " << SymbolicFactorizationReturn
268 <<
" should be " << ExpectedError << std::endl ;
278 Epetra_CrsMatrix* ECM =
dynamic_cast<Epetra_CrsMatrix*
>(MyRowMat) ;
279 int oldtracebackmode = ECM->GetTracebackMode( ) ;
280 ECM->SetTracebackMode(0);
283 ECM->SetTracebackMode(oldtracebackmode);
286 if ( iam == 0 && NumericFactorizationReturn != ExpectedError ) {
287 std::cout <<
" NumericFactorization returned " << NumericFactorizationReturn
288 <<
" should be " << ExpectedError << std::endl ;
301 xexact.PutScalar(1.0);
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 ) )
331 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
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 ) ;
356 Problem.SetLHS( &sAAx );
357 Problem.SetRHS( &b );
365 if ( TrustMe ) sAAx = FixedLHS ;
373 Problem.SetLHS( &sAx );
374 Problem.SetRHS( &sAAx );
377 if ( MyMat->MyGRID( 0 ) )
378 MyMat->SumIntoMyValues( 0, 1, val, ind ) ;
381 if ( TrustMe ) sAx = FixedLHS ;
394 Problem.SetLHS( &x );
395 Problem.SetRHS( &sAx );
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 ) )
424 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
425 MyMatWithDiag->Multiply( transpose, x, kAx ) ;
427 if ( MyMatWithDiag->MyGRID( 0 ) )
428 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
438 if ( MyMatWithDiag->MyGRID( 0 ) )
439 MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
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 ;
virtual int Solve()=0
Solves A X = B (or AT x = B)
T & get(ParameterList &l, const std::string &name)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
virtual int NumericFactorization()=0
Performs NumericFactorization on the matrix A.
virtual int SymbolicFactorization()=0
Performs SymbolicFactorization on the matrix A.
virtual int SetParameters(Teuchos::ParameterList &ParameterList)=0
Updates internal variables.
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Factory for binding a third party direct solver to an Epetra_LinearProblem.
virtual int SetUseTranspose(bool UseTranspose)=0
If set true, X will be set to the solution of AT X = B (not A X = B)
RCP< Epetra_CrsMatrix > NewMatNewMap(Epetra_CrsMatrix &In, int Diagonal, int ReindexRowMap, int ReindexColMap, int RangeMapType, int DomainMapType)
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.
const int NumericallySingularMatrixError
const int StructurallySingularMatrixError
bool isParameter(const std::string &name) const