Amesos Package Browser (Single Doxygen Collection) Development
PerformOneSolveAndTest.cpp
Go to the documentation of this file.
```00001 //
00002 //  OUR_CHK_ERR always returns 1 on error.
00003 //
00004
00005 #define OUR_CHK_ERR(a) { { int epetra_err = a; \
00006                       if (epetra_err != 0) { std::cerr << "Amesos ERROR " << epetra_err << ", " \
00007                            << __FILE__ << ", line " << __LINE__ << std::endl; \
00008 relerror = 1.3e15; relresidual=1e15; return(1);}  }\
00009                    }
00010
00011 #include "Epetra_Comm.h"
00012 #include "Teuchos_ParameterList.hpp"
00013 #include "Amesos.h"
00014 #include "Epetra_CrsMatrix.h"
00015 #include "Epetra_Map.h"
00016 #include "Epetra_Vector.h"
00017 #include "Epetra_LinearProblem.h"
00018 #include "PerformOneSolveAndTest.h"
00019 #include "PartialFactorization.h"
00020 #include "NewMatNewMap.h"
00021 #include "Amesos_TestRowMatrix.h"
00022
00023 //
00024 //  Returns the number of failures.
00025 //  Note:  If AmesosClass is not supported, PerformOneSolveAndTest() will
00026 //  always return 0
00027 //
00028 //  Still have to decide where we are going to check the residual.
00029 //
00030 //  The following table shows the variable names that we use for
00031 //  each of the three phases:
00032 //     compute - which computes the correct value of b
00033 //     solve - which solves for x in  A' A' A x = b
00034 //     check - which computes bcheck = A' A' A x
00035 //
00036 //  For ill-conditioned matrices we restrict the test to one or two
00037 //  solves, by setting Levels to 1 or 2 on input to this routine.
00038 //  When Levels is less than 3, some of the transformations
00039 //  shown in the table as "->" and "<-" are not performed, instead
00040 //  a direct copy is made.
00041 //
00042 //  In the absence of roundoff, each item in a given column should
00043 //  be identical.
00044 //
00045 //  If Levels = 3, we compute and solve A' A' A x = b and hence all
00046 //  of the transformations are performed
00047 //
00048 //  If Levels = 2, the transformations shown in the first column of
00049 //  transformations (labelled Levels>=3) are replaced by a direct copy.
00050 //
00051 //  If Levels = 1, only the transformations shown in the third column
00052 //  are performed, the others being replaced by direct copies.
00053 //
00054 //                           Levels>=3    Levels>=2
00055 //                              A'         A'            A
00056 //  compute             xexact  ->  cAx    ->     cAAx   ->       b
00057 //  solve               x       <-  sAx    <-     sAAx   <-       b
00058 //  check               x       ->  kAx    ->     kAAx   ->  bcheck
00059 //
00060 //  Note that since Levels 2 and 3 use the same A, there is no need to
00061 //  call NumericFactorization() between the second and third call to Solve.
00062 //
00064 //  ====================
00065 //
00066 //  When this code was written, our intent was to have AddToDiag add a constant
00067 //  value to the diagonal even for non-existent diagonal elements.  For now, we have
00068 //  backed off of that.  If we decide to go back to the earlier semantics for
00069 //  AddToDiag (i.e. that it would force all diagonal elements to exist in the
00070 //  matrix, this can be tested by setting AddToAllDiagonalElements to true ;
00071 //
00072 //  See bug #1928 for further discussion.
00073 //
00074 //  ExpectedError
00075 //  =============
00076 //
00077 //  ExpectedError allows us to test for Singular and Structurally Singular matrices
00078 //    ExpectedError = StructurallySingularMatrix or NumericallySingularMatrix
00079
00080
00081 int PerformOneSolveAndTest( const char* AmesosClass,
00082           int EpetraMatrixType,
00083           const Epetra_Comm &Comm,
00084           bool transpose,
00085           bool verbose,
00086           Teuchos::ParameterList ParamList,
00087           Epetra_CrsMatrix *& InMat,
00088           int Levels,
00089           const double Rcond,
00090           double& relerror,
00091           double& relresidual,
00092           int ExpectedError )
00093 {
00094
00095   int ierr = 0;
00096
00098
00099
00100   bool TrustMe = ParamList.get( "TrustMe", false );
00101
00102   bool MyVerbose = false ; //  setting this equal to verbose produces too much output and exceeds the test harnesses 1 Megabyte limit
00103
00104   RCP<Epetra_CrsMatrix> MyMat ;
00105   RCP<Epetra_CrsMatrix> MyMatWithDiag ;
00106
00107   MyMat = rcp( new Epetra_CrsMatrix( *InMat ) );
00108
00109   Amesos_TestRowMatrix ATRW( &*MyMat ) ;
00110
00111   Epetra_RowMatrix* MyRowMat = 0;
00112
00113   assert ( EpetraMatrixType >= 0 && EpetraMatrixType <= 2 );
00114   switch ( EpetraMatrixType ) {
00115   case 0:
00116     MyRowMat = &*MyMat ;
00117     break;
00118   case 1:
00119     MyRowMat = &ATRW;
00120     break;
00121   case 2:
00122     MyMat->OptimizeStorage();
00123     MyRowMat = &*MyMat ;
00124     bool OptStorage = MyMat->StorageOptimized();
00125     assert( OptStorage) ;
00126     break;
00127   }
00128   bool OptStorage = MyMat->StorageOptimized();
00129
00130   Epetra_CrsMatrix* MatPtr = &*MyMat ;
00131
00132   const std::string AC = AmesosClass ;
00133   if ( ExpectedError == 0 ) {
00134     if ( AC != "Amesos_Pardiso" ) {
00135       OUR_CHK_ERR ( PartialFactorization( AmesosClass, Comm, transpose, MyVerbose,
00136             ParamList, MatPtr, Rcond ) );
00137     } else {
00138       if (MyVerbose) std::cout << " AC = "  << AC << " not tested in Partial Factorization" <<std::endl ;   // bug #1915
00139     }
00140   }
00141
00143       int oldtracebackmode = InMat->GetTracebackMode( ) ;
00144
00145     //
00146     //  If AddToDiag is set, create a matrix which is numerically identical, but structurally
00147     //  has no missing diagaonal entries.   In other words, every diagonal element in MyMayWithDiag
00148     //  has an entry in the matrix, though that entry will be zero if InMat has no entry for that
00149     //  particular diagonal element.
00150     //
00151     //  It turns out that this does not actually make sure that all diagonal entries exist because it
00152     //  does not deal with maps that are missing the diagonal element.
00153     //
00154     if ( AddToAllDiagonalElements ) {
00155       MyMatWithDiag = NewMatNewMap( *InMat, 2, 0, 0, 0, 0 );  //  Ensure that all diagonal entries exist ;
00156     } else {
00157       InMat->SetTracebackMode( EPETRA_MIN(oldtracebackmode,1) ) ;   // If the matrix is mimssing diagonal elements, the call to ReplaceDiagonalElements will return 1 indicating missing diagonal elements
00158       MyMatWithDiag = NewMatNewMap( *InMat, 0, 0, 0, 0, 0 );  //  Leave the matrix unchanged structurally
00159     }
00160     //
00162     //
00164     Epetra_Vector Diag( MyMatWithDiag->RowMap() );
00167
00168     ierr = MyMatWithDiag->ExtractDiagonalCopy( Diag );
00169     assert( ierr == 0 );
00170     Diag.Update( 1.0, AddConstVecToDiag, 1.0 ) ;
00171     ierr = MyMatWithDiag->ReplaceDiagonalValues( Diag );   // This may return 1 indicating that structurally non-zero elements were left untouched.
00172     assert( ierr >= 0 );
00173
00174       InMat->SetTracebackMode( oldtracebackmode ) ;
00175   } else {
00176     MyMatWithDiag = rcp( new Epetra_CrsMatrix( *InMat ) );
00177   }
00178
00179   if ( MyVerbose ) std::cout << " Partial Factorization complete " << std::endl ;
00180
00181   relerror = 0 ;
00182   relresidual = 0 ;
00183
00184   assert( Levels >= 1 && Levels <= 3 ) ;
00185
00186   int iam = Comm.MyPID() ;
00187   int errors = 0 ;
00188
00189   const Epetra_Map *RangeMap =
00190     transpose?&MyMat->OperatorDomainMap():&MyMat->OperatorRangeMap() ;
00191   const Epetra_Map *DomainMap =
00192     transpose?&MyMat->OperatorRangeMap():&MyMat->OperatorDomainMap() ;
00193
00194   Epetra_Vector xexact(*DomainMap);
00195   Epetra_Vector x(*DomainMap);
00196
00197   Epetra_Vector cAx(*DomainMap);
00198   Epetra_Vector sAx(*DomainMap);
00199   Epetra_Vector kAx(*DomainMap);
00200
00201   Epetra_Vector cAAx(*DomainMap);
00202   Epetra_Vector sAAx(*DomainMap);
00203   Epetra_Vector kAAx(*DomainMap);
00204
00205   Epetra_Vector FixedLHS(*DomainMap);
00206   Epetra_Vector FixedRHS(*RangeMap);
00207
00208   Epetra_Vector b(*RangeMap);
00209   Epetra_Vector bcheck(*RangeMap);
00210
00211   Epetra_Vector DomainDiff(*DomainMap);
00212   Epetra_Vector RangeDiff(*RangeMap);
00213
00214   Epetra_LinearProblem Problem;
00215   Teuchos::RCP<Amesos_BaseSolver> Abase ;
00216   Amesos Afactory;
00217
00218
00219
00220
00221   Abase = Teuchos::rcp(Afactory.Create( AmesosClass, Problem )) ;
00222
00223   relerror = 0 ;
00224   relresidual = 0 ;
00225
00226   if ( Abase == Teuchos::null )
00227     assert( false ) ;
00228   else {
00229
00230     //
00231     //  Phase 1:  Compute b = A' A' A xexact
00232     //
00233     Problem.SetOperator( MyRowMat );
00234     //    Epetra_CrsMatrix* ECM = dynamic_cast<Epetra_CrsMatrix*>(MyRowMat) ;
00235
00236     //
00237     //  We only set transpose if we have to - this allows valgrind to check
00238     //  that transpose is set to a default value before it is used.
00239     //
00240     if ( transpose ) OUR_CHK_ERR( Abase->SetUseTranspose( transpose ) );
00241     if (MyVerbose) ParamList.set( "DebugLevel", 1 );
00242     if (MyVerbose) ParamList.set( "OutputLevel", 1 );
00243     OUR_CHK_ERR( Abase->SetParameters( ParamList ) );
00244
00245     if ( TrustMe ) {
00246       Problem.SetLHS( &FixedLHS );
00247       Problem.SetRHS( &FixedRHS );
00248       assert( OptStorage) ;
00249     }
00250
00251     // bug #2184
00252     //  Structurally singular matrices are not detected by
00253     //  Amesos_Klu::SymbolicFactorization() but they are detected by
00254     //  Amesos_Klu::NumericFactorization()
00255     if ( ExpectedError == StructurallySingularMatrixError )
00256       ExpectedError = NumericallySingularMatrixError ;
00257     if ( ExpectedError == StructurallySingularMatrixError ) {
00258       Epetra_CrsMatrix* ECM = dynamic_cast<Epetra_CrsMatrix*>(MyRowMat) ;
00259       int oldtracebackmode = ECM->GetTracebackMode( ) ;
00260       ECM->SetTracebackMode(0);  // We expect an error, but we don't want it to print out
00261
00262       const int SymbolicFactorizationReturn = Abase->SymbolicFactorization(  ) ;
00263       ECM->SetTracebackMode(oldtracebackmode);
00264       // bug #2245 - Amesos fails to return error consistently across all
00265       // processes.  When this bug is fixed, remove "iam == 0 &&" from the next line
00266       if ( iam == 0 && SymbolicFactorizationReturn != ExpectedError ) {
00267   std::cout << " SymbolicFactorization returned " << SymbolicFactorizationReturn
00268        << " should be " << ExpectedError << std::endl ;
00269   OUR_CHK_ERR( 1 ) ;
00270       } else {
00271   return 0;   //  Returned the correct error code for this matrix
00272       }
00273     } else {
00274       const int SymbolicFactorizationReturn = Abase->SymbolicFactorization(  ) ;
00275       OUR_CHK_ERR( SymbolicFactorizationReturn ) ;
00276     }
00277     if ( ExpectedError == NumericallySingularMatrixError ) {
00278       Epetra_CrsMatrix* ECM = dynamic_cast<Epetra_CrsMatrix*>(MyRowMat) ;
00279       int oldtracebackmode = ECM->GetTracebackMode( ) ;
00280       ECM->SetTracebackMode(0); // We expect an error, but we don't want it to print out
00281
00282       const int NumericFactorizationReturn = Abase->NumericFactorization(  ) ;
00283       ECM->SetTracebackMode(oldtracebackmode);
00284       // bug #2245 - Amesos fails to return error consistently across all
00285       // processes.  When this bug is fixed, remove "iam == 0 &&" from the next line
00286       if ( iam == 0 && NumericFactorizationReturn != ExpectedError ) {
00287   std::cout << " NumericFactorization returned " << NumericFactorizationReturn
00288        << " should be " << ExpectedError << std::endl ;
00289   OUR_CHK_ERR( 1 ) ;
00290       } else {
00291   return 0;   //  Returned the correct error code for this matrix
00292       }
00293     } else {
00294       const int NumericFactorizationReturn = Abase->NumericFactorization(  ) ;
00295       OUR_CHK_ERR( NumericFactorizationReturn ) ;
00296     }
00297     int ind[1];
00298     double val[1];
00299     ind[0] = 0;
00300     xexact.Random();
00301     xexact.PutScalar(1.0);
00302
00303     //
00304     //  Compute cAx = A' xexact
00305     //
00306     double Value = 1.0 ;
00307     if ( Levels == 3 )
00308       {
00309   val[0] = Value ;
00310   if ( MyMatWithDiag->MyGRID( 0 ) ) {
00311     MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
00312   }
00313   MyMatWithDiag->Multiply( transpose, xexact, cAx ) ;
00314
00315   val[0] = - Value ;
00316   if ( MyMatWithDiag->MyGRID( 0 ) )
00317     MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
00318       }
00319     else
00320       {
00321   cAx = xexact ;
00322       }
00323
00324     //
00325     //  Compute cAAx = A' cAx
00326     //
00327     if ( Levels >= 2 )
00328       {
00329   val[0] =  Value ;
00330   if ( MyMatWithDiag->MyGRID( 0 ) )
00331     MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
00332   MyMatWithDiag->Multiply( transpose, cAx, cAAx ) ; //  x2 = A' x1
00333
00334   val[0] = - Value ;
00335   if ( MyMatWithDiag->MyGRID( 0 ) )
00336     MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
00337       }
00338     else
00339       {
00340   cAAx = cAx ;
00341       }
00342
00343     if ( MyVerbose ) std::cout << " Compute  b = A x2 = A A' A'' xexact  " << std::endl ;
00344
00345     MyMatWithDiag->Multiply( transpose, cAAx, b ) ;  //  b = A x2 = A A' A'' xexact
00346
00347
00348     //  Phase 2:  Solve A' A' A x = b
00349     //
00350     //
00351     //  Solve A sAAx = b
00352     //
00353     if ( TrustMe ) {
00354       FixedRHS = b;
00355     } else {
00356       Problem.SetLHS( &sAAx );
00357       Problem.SetRHS( &b );
00358     }
00359
00360
00361     OUR_CHK_ERR( Abase->SymbolicFactorization(  ) );
00362     OUR_CHK_ERR( Abase->SymbolicFactorization(  ) );     // This should be irrelevant, but should nonetheless be legal
00363     OUR_CHK_ERR( Abase->NumericFactorization(  ) );
00364     OUR_CHK_ERR( Abase->Solve(  ) );
00365     if ( TrustMe ) sAAx = FixedLHS ;
00366
00367     if ( Levels >= 2 )
00368       {
00369         OUR_CHK_ERR( Abase->SetUseTranspose( transpose ) );
00370   if ( TrustMe ) {
00371     FixedRHS = sAAx ;
00372   } else {
00373     Problem.SetLHS( &sAx );
00374     Problem.SetRHS( &sAAx );
00375   }
00376   val[0] =  Value ;
00377   if ( MyMat->MyGRID( 0 ) )
00378     MyMat->SumIntoMyValues( 0, 1, val, ind ) ;
00379   OUR_CHK_ERR( Abase->NumericFactorization(  ) );
00380   OUR_CHK_ERR( Abase->Solve(  ) );
00381   if ( TrustMe ) sAx = FixedLHS ;
00382
00383       }
00384     else
00385       {
00386   sAx = sAAx ;
00387       }
00388
00389     if ( Levels >= 3 )
00390       {
00391   if ( TrustMe ) {
00392     FixedRHS = sAx ;
00393   } else {
00394     Problem.SetLHS( &x );
00395     Problem.SetRHS( &sAx );
00396   }
00397   OUR_CHK_ERR( Abase->Solve(  ) );
00398   if ( TrustMe ) x = FixedLHS ;
00399       }
00400     else
00401       {
00402   x = sAx ;
00403       }
00404
00405     if ( Levels >= 2 )
00406       {
00407   val[0] =  -Value ;
00408   if ( MyMat->MyGRID( 0 ) ) {
00409     if ( MyMat->SumIntoMyValues( 0, 1, val, ind ) ) {
00410       std::cout << " TestOptions requires a non-zero entry in A(1,1) " << std::endl ;
00411     }
00412   }
00413       }
00414
00415     //
00416     //  Phase 3:  Check the residual: bcheck = A' A' A x
00417     //
00418
00419
00420     if ( Levels >= 3 )
00421       {
00422   val[0] =  Value ;
00423   if ( MyMatWithDiag->MyGRID( 0 ) )
00424     MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
00425   MyMatWithDiag->Multiply( transpose, x, kAx ) ;
00426   val[0] =  -Value ;
00427   if ( MyMatWithDiag->MyGRID( 0 ) )
00428     MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
00429       }
00430     else
00431       {
00432   kAx = x ;
00433       }
00434
00435     if ( Levels >= 2 )
00436       {
00437   val[0] =  Value ;
00438   if ( MyMatWithDiag->MyGRID( 0 ) )
00439     MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
00440   MyMatWithDiag->Multiply( transpose, kAx, kAAx ) ;
00441   val[0] =  -Value ;
00442   if ( MyMatWithDiag->MyGRID( 0 ) )
00443     MyMatWithDiag->SumIntoMyValues( 0, 1, val, ind ) ;
00444       }
00445     else
00446       {
00447   kAAx = kAx ;
00448       }
00449
00450     MyMatWithDiag->Multiply( transpose, kAAx, bcheck ) ; //  temp = A" x2
00451
00452     double norm_diff ;
00453     double norm_one ;
00454
00455     DomainDiff.Update( 1.0, sAAx, -1.0, cAAx, 0.0 ) ;
00456     DomainDiff.Norm2( &norm_diff ) ;
00457     sAAx.Norm2( &norm_one ) ;
00458     if (MyVerbose) std::cout << __FILE__ << "::" << __LINE__
00459       << " norm( sAAx - cAAx ) / norm(sAAx ) = "
00460       << norm_diff /norm_one << std::endl ;
00461
00462
00463
00464
00465     DomainDiff.Update( 1.0, sAx, -1.0, cAx, 0.0 ) ;
00466     DomainDiff.Norm2( &norm_diff ) ;
00467     sAx.Norm2( &norm_one ) ;
00468     if (MyVerbose) std::cout
00469       << __FILE__ << "::" << __LINE__
00470       << " norm( sAx - cAx ) / norm(sAx ) = "
00471           << norm_diff /norm_one << std::endl ;
00472
00473
00474     DomainDiff.Update( 1.0, x, -1.0, xexact, 0.0 ) ;
00475     DomainDiff.Norm2( &norm_diff ) ;
00476     x.Norm2( &norm_one ) ;
00477     if (MyVerbose) std::cout
00478       << __FILE__ << "::" << __LINE__
00479       << " norm( x - xexact ) / norm(x) = "
00480       << norm_diff /norm_one << std::endl ;
00481
00482     relerror = norm_diff / norm_one ;
00483
00484     DomainDiff.Update( 1.0, sAx, -1.0, kAx, 0.0 ) ;
00485     DomainDiff.Norm2( &norm_diff ) ;
00486     sAx.Norm2( &norm_one ) ;
00487     if (MyVerbose) std::cout
00488       << __FILE__ << "::" << __LINE__
00489       << " norm( sAx - kAx ) / norm(sAx ) = "
00490       << norm_diff /norm_one << std::endl ;
00491
00492
00493     DomainDiff.Update( 1.0, sAAx, -1.0, kAAx, 0.0 ) ;
00494     DomainDiff.Norm2( &norm_diff ) ;
00495     sAAx.Norm2( &norm_one ) ;
00496     if (MyVerbose) std::cout
00497       << __FILE__ << "::" << __LINE__
00498       << " norm( sAAx - kAAx ) / norm(sAAx ) = "
00499       << norm_diff /norm_one << std::endl ;
00500
00501
00502     RangeDiff.Update( 1.0, bcheck, -1.0, b, 0.0 ) ;
00503     RangeDiff.Norm2( &norm_diff ) ;
00504     bcheck.Norm2( &norm_one ) ;
00505     if (MyVerbose) std::cout
00506       << __FILE__ << "::" << __LINE__
00507       << " norm( bcheck - b ) / norm(bcheck ) = "
00508       << norm_diff /norm_one << std::endl ;
00509
00510     relresidual = norm_diff / norm_one ;
00511
00512     if (iam == 0 ) {
00513       if ( relresidual * Rcond < 1e-16 ) {
00514   if (MyVerbose) std::cout << " Test 1 Passed " << std::endl ;
00515       } else {
00516       std::cout <<  __FILE__ << "::"  << __LINE__ <<
00517     " relresidual = " << relresidual <<
00518     " TEST FAILED " <<
00519     " ParamList = " << ParamList << std::endl ;
00520   errors += 1 ;
00521       }
00522     }
00523   }
00524
00525   return errors;
00526
00527 }
00528
00529
```