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