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 //   
00063 //  On testing AddToDiag
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 
00097   bool AddToAllDiagonalElements =  ParamList.get( "AddZeroToDiag", false ) ;
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 
00142   if (ParamList.isParameter("AddToDiag") ) { 
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     //
00161     //  Now add AddToDiag to each diagonal element.  
00162     //
00163     double AddToDiag = ParamList.get("AddToDiag", 0.0 );
00164     Epetra_Vector Diag( MyMatWithDiag->RowMap() );
00165     Epetra_Vector AddConstVecToDiag( MyMatWithDiag->RowMap() );
00166     AddConstVecToDiag.PutScalar( AddToDiag );
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines