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() );
00165     AddConstVecToDiag.PutScalar( AddToDiag );
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines