Amesos Package Browser (Single Doxygen Collection) Development
TestOptions.cpp
Go to the documentation of this file.
00001 //
00002 //  To run this under valgrind, try:
00003 //  valgrind --suppressions=../Test_Basic/Suppressions --gen-suppressions=yes --leak-check=yes --show-reachable=yes ./TestOptions.exe -v
00004 //
00005 //  To run this with valgrind under mpirun, 
00006 //  mpirun -np 2 valgrind --log-file=TestOpt.logfile --suppressions=../Test_Basic/Suppressions --gen-suppressions=yes --leak-check=yes --show-reachable=yes ./TestOptions.exe -v
00007 //
00008 //  test/scripts/daily/serial/TestMemoryLeaks[.exe] performs an automated test for memory leaks
00009 //  using valgrind and this code.  To run TestMemoryLeaks, cd to test/TestOptions in the
00010 //  build directory and type ../scripts/daily/serial/TestMemoryLeaks.exe.  The output is stored
00011 //  in amesos/logLinux.txt.
00012 //
00013 //
00014 
00015 //  TestOptions tests all options for each Amesos Class on a limited number 
00016 //  of matrices.  
00017 //
00018 
00019 //  TestOneMatrix - Test one matrix 
00020 //    - Distributed vs not distributed -
00021 //    - Transpose vs not transposed -
00022 //      TestAllClasses  - Test one matrix and one setting of distributed and transpose
00023 //        - Calls TestOtherClasses (one for each Amesos class) and TestSuperludist 
00024 //        TestOtherClasses
00025 //        TestSuperludist
00026 //        TestScalapack
00027 //
00028 //
00029 //  Todo:
00030 //    Write TestKlu, TestSuperlu, TestScalapack, TestUmfpack, TestDscpack
00031 //    Enable tests of various parameter options
00032 //    Make it test all four codes (DSCPACK, UMFPACK, SuperLU_DIST, KLU )
00033 //    Valgrind it
00034 //    Enable tests of transpose and distributed matrices  - DONE 
00035 //    Enable FACTOR_B - DONE 
00036 //
00037 
00038 #include "Trilinos_Util.h"
00039 #include "Trilinos_Util_ReadMatrixMarket2Epetra.h"
00040 #include "Amesos.h"
00041 #include "Teuchos_ParameterList.hpp"
00042 #include "Amesos_BaseSolver.h"
00043 #include "Epetra_LinearProblem.h"
00044 #include "Epetra_CrsMatrix.h"
00045 #include "Epetra_Map.h"
00046 #include "Epetra_Vector.h"
00047 #include "Epetra_Export.h"
00048 #include "Amesos_Umfpack.h"
00049 #include "CrsMatrixTranspose.h"
00050 #include "TestAllClasses.h"
00051 #include <string>
00052 #include "Teuchos_RCP.hpp"
00053 #include "NewMatNewMap.h"
00054 #ifdef EPETRA_MPI
00055 #include "Epetra_MpiComm.h"
00056 #else
00057 #include "Epetra_SerialComm.h"
00058 #endif
00059 
00060 #if 0
00061 
00062 #ifdef HAVE_VALGRIND_H
00063 #include <valgrind/valgrind.h>
00064 #define HAVE_VALGRIND
00065 #else
00066 #ifdef HAVE_VALGRIND_VALGRIND_H
00067 #include <valgrind/valgrind.h>
00068 #define HAVE_VALGRIND
00069 #endif 
00070 #endif 
00071 
00072 #endif 
00073 
00074 std::vector<std::string> AmesosClasses;
00075 
00076 int NumAmesosClasses;
00077 
00078 int CreateCrsMatrix( char *in_filename, const Epetra_Comm &Comm, 
00079          Epetra_Map *& readMap,
00080          const bool transpose, const bool distribute, 
00081          bool& symmetric, Epetra_CrsMatrix *& Matrix ) {
00082 
00083   Epetra_CrsMatrix * readA = 0; 
00084   Epetra_Vector * readx = 0; 
00085   Epetra_Vector * readb = 0;
00086   Epetra_Vector * readxexact = 0;
00087 
00088   //
00089   //  This hack allows TestOptions to be run from either the test/TestOptions/ directory or from 
00090   //  the test/ directory (as it is in nightly testing and in make "run-tests")
00091   //
00092   FILE *in_file = fopen( in_filename, "r");
00093 
00094   char *filename;
00095   if (in_file == NULL ) 
00096     filename = &in_filename[1] ; //  Strip off ithe "." from
00097          //  "../" and try again 
00098   else {
00099     filename = in_filename ;
00100     fclose( in_file );
00101   }
00102 
00103   symmetric = false ; 
00104   std::string FileName = filename ;
00105 
00106   int FN_Size = FileName.size() ; 
00107   std::string LastFiveBytes = FileName.substr( EPETRA_MAX(0,FN_Size-5), FN_Size );
00108   std::string LastFourBytes = FileName.substr( EPETRA_MAX(0,FN_Size-4), FN_Size );
00109 
00110   if ( LastFiveBytes == ".triU" ) { 
00111     // Call routine to read in unsymmetric Triplet matrix
00112     EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename, false, Comm, readMap, readA, readx, 
00113                   readb, readxexact) );
00114     symmetric = false; 
00115   } else {
00116     if ( LastFiveBytes == ".triS" ) { 
00117       // Call routine to read in symmetric Triplet matrix
00118       EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename, true, Comm, readMap, readA, readx, 
00119               readb, readxexact) );
00120       symmetric = true; 
00121     } else {
00122       if (  LastFourBytes == ".mtx" ) { 
00123   EPETRA_CHK_ERR( Trilinos_Util_ReadMatrixMarket2Epetra( filename, Comm, readMap, 
00124                      readA, readx, readb, readxexact) );   
00125   in_file = fopen( filename, "r");
00126   assert (in_file != NULL) ;  // Checked in Trilinos_Util_CountMatrixMarket() 
00127   const int BUFSIZE = 800 ; 
00128   char buffer[BUFSIZE] ; 
00129   fgets( buffer, BUFSIZE, in_file ) ;  // Pick symmetry info off of this std::string 
00130   std::string headerline1 = buffer;
00131 #ifdef TFLOP
00132   if ( headerline1.find("symmetric") < BUFSIZE ) symmetric = true;
00133 #else
00134   if ( headerline1.find("symmetric") != std::string::npos) symmetric = true; 
00135 
00136 #endif
00137   fclose(in_file);
00138 
00139       } else {
00140   // Call routine to read in HB problem
00141   Trilinos_Util_ReadHb2Epetra( filename, Comm, readMap, readA, readx, 
00142                  readb, readxexact) ;
00143   if (  LastFourBytes == ".rsa" ) symmetric = true ; 
00144       }
00145     }
00146   }
00147 
00148 
00149   if ( readb )  delete readb;
00150   if ( readx ) delete readx;
00151   if ( readxexact ) delete readxexact;
00152 
00153   Epetra_CrsMatrix *serialA ; 
00154   Epetra_CrsMatrix *transposeA;
00155 
00156   if ( transpose ) {
00157     transposeA = new Epetra_CrsMatrix( Copy, *readMap, 0 );
00158     assert( CrsMatrixTranspose( readA, transposeA ) == 0 ); 
00159     serialA = transposeA ; 
00160     delete readA;
00161     readA = 0 ; 
00162   } else {
00163     serialA = readA ; 
00164   }
00165 
00166   assert( (void *) &serialA->Graph() ) ;
00167   assert( (void *) &serialA->RowMap() ) ;
00168   assert( serialA->RowMap().SameAs(*readMap) ) ; 
00169 
00170   if ( distribute ) { 
00171     // Create uniform distributed map
00172     Epetra_Map DistMap(readMap->NumGlobalElements(), 0, Comm);
00173 
00174     // Create Exporter to distribute read-in matrix and vectors
00175     Epetra_Export exporter( *readMap, DistMap );
00176     
00177     Epetra_CrsMatrix *Amat = new Epetra_CrsMatrix( Copy, DistMap, 0 );
00178     Amat->Export(*serialA, exporter, Add);
00179     assert(Amat->FillComplete()==0);    
00180     
00181     Matrix = Amat; 
00182     //
00183     //  Make sure that deleting Amat->RowMap() will delete map 
00184     //
00185     //  Bug:  We can't manage to delete map his way anyway,
00186     //        and this fails on tranposes, so for now I just accept
00187     //        the memory loss.
00188     //    assert( &(Amat->RowMap()) == map ) ; 
00189     delete readMap; 
00190     readMap = 0 ; 
00191     delete serialA; 
00192   } else { 
00193 
00194     Matrix = serialA; 
00195   }
00196 
00197 
00198   return 0;
00199 }
00200 
00201 int TestErrors( const std::vector<bool> AmesosClassesInstalled, 
00202        char *filename, 
00203 #ifdef EPETRA_MPI
00204        Epetra_MpiComm& Comm,
00205 #else
00206        Epetra_SerialComm& Comm,
00207 #endif
00208        const bool verbose, 
00209        int &NumTests  ) {
00210 
00211   int NumErrors =0 ;
00212   double error = -13; 
00213   double residual = -13;
00214 
00215 
00216   for ( int iterTrans =0 ; iterTrans < 2; iterTrans++ ) {
00217     const bool transpose = iterTrans == 1 ; 
00218     
00219     const bool distribute = 1;
00220 
00221     const int iterRowindex = 0;
00222     const int iterColindex = 0 ;
00223     const int iterRangemap = 0 ;
00224     const int iterDomainmap = 0 ;
00225     const int iterDiagonalOpts = 0 ; 
00226     bool printit = true ; 
00227     if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ; 
00228     const int EpetraMatrixType = 0 ;
00229     bool symmetric = true;
00230     const int iterDist = 0 ; 
00231     
00232     Epetra_CrsMatrix *Amat = 0 ;
00233     Epetra_Map *readMap = 0 ;
00234     CreateCrsMatrix( filename, Comm, readMap, transpose, distribute, symmetric, Amat ) ;
00235     
00236     if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << 
00237         " Creating matrix from " <<
00238         " filename = " << filename <<
00239         " symmetric = " << symmetric <<
00240         " distribute = " << distribute <<
00241         " iterRowindex = " << iterRowindex <<
00242         " iterColindex = " << iterColindex <<
00243         " iterRangemap = " << iterRangemap <<
00244         " iterDomainmap = " << iterDomainmap <<
00245         " EpetraMatrixType = " << EpetraMatrixType <<
00246         " iterDiagonalOpts = " << iterDiagonalOpts <<
00247         " transpose = "  << transpose 
00248            << " iterDist = " << iterDist << std::endl ; 
00249     
00250     if ( iterDiagonalOpts )  Comm.SetTracebackMode(1);  // Turn off positive Epetra warnings (i.e. iniefficient code, such as memory re-allocation)
00251     
00252     if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ; 
00253     
00254     RCP<Epetra_CrsMatrix> Bmat = NewMatNewMap( *Amat, 
00255                    iterDiagonalOpts, 
00256                    iterRowindex,
00257                    iterColindex,
00258                    iterRangemap,
00259                    iterDomainmap
00260                    ) ; 
00261     if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ; 
00262     Comm.SetTracebackMode(2);
00263     
00264     if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << 
00265         " filename = " << filename <<
00266         " symmetric = " << symmetric <<
00267         " distribute = " << distribute <<
00268         " iterRowindex = " << iterRowindex <<
00269         " iterColindex = " << iterColindex <<
00270         " iterRangemap = " << iterRangemap <<
00271         " iterDomainmap = " << iterDomainmap <<
00272         " EpetraMatrixType = " << EpetraMatrixType <<
00273         " iterDiagonalOpts = " << iterDiagonalOpts <<
00274         " transpose = "  << transpose << " iterDist = " << iterDist << std::endl ; 
00275     
00276     //
00277     //  This causes a failure in Amesos_Superludist:
00278     Epetra_CrsMatrix* Cmat = &*Bmat;
00279     //  Epetra_CrsMatrix* Cmat = Amat ;
00280     
00281     
00282     const int Level = 1; 
00283     const double MaxError = 1e-3;
00284     
00285     int NumTheseTests = 0 ; 
00286     if ( verbose ) {
00287       std::cout << " About to test  " << filename 
00288      << __FILE__ << "::"  << __LINE__
00289      << " EpetraMatrixType = " <<  EpetraMatrixType 
00290      << (transpose?" transpose":"" ) 
00291      << (distribute?" distribute":"" )
00292      << std::endl ; 
00293     }
00294     int Error = TestAllClasses( AmesosClasses, EpetraMatrixType, 
00295         AmesosClassesInstalled, 
00296         Cmat, 
00297         transpose ,
00298         verbose, 
00299         symmetric, 
00300         Level,
00301         MaxError, 
00302         iterDiagonalOpts, 
00303         iterRowindex,
00304         iterColindex,
00305         iterRangemap,
00306         iterDomainmap,
00307         distribute,
00308         filename,
00309         error, 
00310         residual, 
00311         NumTheseTests ) ;
00312     NumTests += NumTheseTests ;
00313     NumErrors += Error ;
00314     //      BUG:  Memory leak 
00315     //      delete &(Amat->RowMap()) ; 
00316     if ( Amat ) delete Amat ; 
00317     if ( readMap ) delete readMap ; 
00318   }
00319   if ( verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ; 
00320 
00321   return NumErrors;
00322 } 
00323 
00324 int TestOneMatrix( const std::vector<bool> AmesosClassesInstalled, 
00325        char *filename, 
00326 #ifdef EPETRA_MPI
00327        Epetra_MpiComm& Comm,
00328 #else
00329        Epetra_SerialComm& Comm,
00330 #endif
00331        //      Epetra_Comm &Comm, 
00332        const bool verbose, 
00333        const bool PerformDiagonalTest, 
00334        double Rcond,
00335        int &NumTests  ) {
00336 
00337   int NumErrors =0 ;
00338   double error = -13; 
00339   double residual = -13;
00340   //  double errors[NumAmesosClasses];
00341   //  double residuals[NumAmesosClasses];
00342   //  for (int i = 0 ; i < NumAmesosClasses; i ++ ) errors[i] = residuals[i] = 0.0 ; 
00343 
00344   //#ifdef HAVE_AMESOS_UMFPACK
00345 #if 0
00346   Epetra_CrsMatrix *Amat ;
00347 
00348   //
00349   //  Compute the reciprocal condition number using Amesos_UMFPACK via the Amesos interface
00350   //
00351   Epetra_Map *readMap = 0 ; 
00352   CreateCrsMatrix( filename, Comm, readMap, false, false, symmetric, Amat ) ;
00353   Teuchos::ParameterList ParamList ;
00354   Epetra_LinearProblem Problem;
00355   Amesos Afactory;
00356 
00357   Amesos_BaseSolver* Abase ; 
00358   Abase = Afactory.Create( "Amesos_Umfpack", Problem ) ; 
00359   if ( Abase == 0 ) {
00360     std::cerr << " AMESOS_UMFPACK is required for this test " << std::endl ;
00361     exit(13);
00362   }  ;
00363 
00364   //
00365   //  Factor A to compute Rcond = reciprocal condition number estimate
00366   //
00367   Problem.SetOperator( Amat );
00368   EPETRA_CHK_ERR( Abase->SymbolicFactorization(  ) ); 
00369   EPETRA_CHK_ERR( Abase->NumericFactorization(  ) ); 
00370   Amesos_Umfpack* UmfpackOperator = dynamic_cast<Amesos_Umfpack *> (Abase) ; 
00371   //  double Rcond = UmfpackOperator->GetRcond();
00372 
00373   int ind[1];
00374   double val[1];
00375   ind[0] = 0;
00376   val[0] = 1 ; 
00377   double AnormInf =  Amat->NormInf() ;
00378   if (verbose) std::cout << " norm(Amat) = " << AnormInf << std::endl; 
00379   if ( Amat->MyGRID( 0 ) )
00380     Amat->SumIntoMyValues( 0, 1, val, ind ) ; 
00381   AnormInf =  Amat->NormInf() ;
00382   if (verbose) std::cout << " norm(Amat) = " << AnormInf << std::endl; 
00383 
00384 
00385   EPETRA_CHK_ERR( Abase->SymbolicFactorization(  ) ); 
00386   EPETRA_CHK_ERR( Abase->NumericFactorization(  ) ); 
00387   double Rcond1 = UmfpackOperator->GetRcond();
00388 
00389   if ( Amat->MyGRID( 0 ) )
00390     Amat->SumIntoMyValues( 0, 1, val, ind ) ; 
00391   AnormInf =  Amat->NormInf() ;
00392   if (verbose) std::cout << " norm(Amat) = " << AnormInf << std::endl; 
00393   EPETRA_CHK_ERR( Abase->SymbolicFactorization(  ) ); 
00394   EPETRA_CHK_ERR( Abase->NumericFactorization(  ) ); 
00395    double Rcond2 = UmfpackOperator->GetRcond();
00396 
00397   if (verbose) std::cout << " Rcond1 = " << Rcond1 << std::endl; 
00398   if (verbose) std::cout << " Rcond2 = " << Rcond2 << std::endl; 
00399 
00400   if ( readMap ) delete readMap ;
00401 #else
00402   double Rcond1 = Rcond ;
00403   double Rcond2 = Rcond ;
00404 #endif
00405 
00406   //
00407   //  Rowindex and Colindex control the maps and indices used to create the matrix
00408   //
00409   //  These tests are all disabled in TestAllClasses.cpp
00410   //
00411   const int RowindexMax = 3;   // bug should be three ( 1 based, 3 based, non contiguous )
00412   const int ColindexMax = 2;   // bug should be two:  ( row map, 4 based )
00413 
00414   //
00415   //  Rangemap and Domainmap control the Range and Domain maps used in the call to FillComplete
00416   //  If both are "no change", FillComplete is called with no parameters (i.e. without specifying maps)
00417   //  Otherwise, domain and range maps are specified in the call to FillComplete
00418   //
00419   //  These tests are all disabled in TestAllClasses.cpp
00420   //
00421   int RangemapMax = 4; // bug should be four:  ( no change, serial, bizarre dist, replicated )
00422   int DomainmapMax = 4; // bug should be four:  ( no change, serial, bizarre dist, replicated )  IRRELEVANT see ThisDomainMax  
00423 
00424   //
00425   //  DiagonalOpts controls whether diagonal elements are left alone,
00426   //  or removed from both the matrix and the underlying map
00427   //
00428   int DiagonalOptsMax = 2;   // should be two:  ( no change, elements missing from map )
00429   //
00430   //
00431   //
00432   int EpetraMatrixTypeMax = 3; // 0 = Epetra_CrsMatrix; 1 = Epetra_RowMatriw; 2 = StorageOptimized Epetra_CrsMatrix
00433   //
00434   //  No point in trying to run distributed memory tests on a serial run
00435   //
00436   int iterDistMax = 2;
00437   if ( Comm.NumProc() == 1 ) {
00438     iterDistMax = 1 ; 
00439     RangemapMax = 1 ; 
00440     DomainmapMax = 1 ; 
00441   }
00442 
00443   
00444 
00445   if (! PerformDiagonalTest ) DiagonalOptsMax = 1 ; 
00446 
00447   for ( int iterTrans =0 ; iterTrans < 2; iterTrans++ ) {
00448     bool transpose = iterTrans == 1 ; 
00449     
00450     for ( int iterDist =0; iterDist < iterDistMax; iterDist++ ) {  
00451       bool distribute = ( iterDist == 1 ); 
00452   
00453 #if 1
00454       for ( int iterRowindex = 0 ; iterRowindex < RowindexMax; iterRowindex++ ) {
00455   for ( int iterColindex = 0 ; iterColindex < ColindexMax; iterColindex++ ) {
00456     //
00457     //  The current version of NewMatNewMap.cpp supports only trivial 
00458     //  replicated maps, hence we do not allow any fancy indexing
00459     //
00460     int ThisRangemapMax = RangemapMax ;
00461               // Bug #1920 Amesos classes can't handle replicated domain or ranges  if ( iterRowindex > 0 || iterColindex > 0 ) 
00462     ThisRangemapMax = EPETRA_MIN( 3, ThisRangemapMax );
00463     int ThisDomainmapMax =  EPETRA_MIN( 3, ThisRangemapMax );  // Bug #1920 
00464     for ( int iterRangemap = 0 ; iterRangemap < ThisRangemapMax; iterRangemap++ ) {
00465       for ( int iterDomainmap = 0 ; iterDomainmap < ThisDomainmapMax; iterDomainmap++ ) {
00466         for ( int iterDiagonalOpts = 0 ; iterDiagonalOpts < DiagonalOptsMax; iterDiagonalOpts++ ) { 
00467 #else
00468     int iterRowindex = 0; { 
00469       int iterColindex = 0 ; { 
00470         int iterRangemap = 0 ; { 
00471           int iterDomainmap = 0 ; {
00472       for ( int iterDiagonalOpts = 1 ; iterDiagonalOpts < DiagonalOptsMax; iterDiagonalOpts++ ) { 
00473 #endif
00474         const bool printit = false; 
00475     //  diagonal opts testing only works on distributed matrices whose row and column indices match 
00476     //  On a serial matrix, eliminate a column from the map makes the matrix singular
00477     //  If the row and column indices don't match, eliminating a column from the map is, typically, irrelevant
00478 
00479     if ( ( iterColindex == 0 && distribute ) || iterDiagonalOpts == 0 ) { 
00480       if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ; 
00481       for ( int EpetraMatrixType = 0 ; EpetraMatrixType < EpetraMatrixTypeMax;  EpetraMatrixType++ ) {
00482 
00483       if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ; 
00484         //  These tests presently take over 7 hours on some platforms.  But, I don't want to eliminate any category of tests
00485         //  The following test will cull 90% of the tests and still cover every type of test and most combinations
00486         Epetra_Util EU;
00487         int RunTest[1] ;
00488         RunTest[0] =  (EU.RandomDouble() > 0.8)?1:0 ; 
00489         if ( iterRowindex == 0 &&
00490        iterColindex == 0 &&
00491        iterRangemap == 0 &&
00492        iterDomainmap == 0 ) RunTest[0] = 1; 
00493         Comm.Broadcast( RunTest, 1, 0 ) ; 
00494         if ( RunTest[0] ) { 
00495         //
00496         //  We test only one level for different indexing or different Range and Domain maps
00497         //  to avoid hassles of moving data from the domain space to the range space and back
00498         //
00499       if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ; 
00500         int MaxLevel = 3 ; 
00501         if ( iterRowindex > 0 ) MaxLevel = 1 ; 
00502         if ( iterColindex > 0 ) MaxLevel = 1 ; 
00503         if ( iterRangemap > 0 ) MaxLevel = 1 ; 
00504         if ( iterDomainmap > 0 ) MaxLevel = 1 ; 
00505 
00506         bool symmetric = true;
00507         
00508         Epetra_CrsMatrix *Amat = 0 ;
00509         Epetra_Map *readMap = 0 ;
00510         CreateCrsMatrix( filename, Comm, readMap, transpose, distribute, symmetric, Amat ) ;
00511         //      assert( symmetric == false ) ; 
00512         
00513         if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << 
00514              " Creating matrix from " <<
00515              " filename = " << filename <<
00516              " symmetric = " << symmetric <<
00517              " distribute = " << distribute <<
00518              " iterRowindex = " << iterRowindex <<
00519              " iterColindex = " << iterColindex <<
00520              " iterRangemap = " << iterRangemap <<
00521              " iterDomainmap = " << iterDomainmap <<
00522              " EpetraMatrixType = " << EpetraMatrixType <<
00523              " iterDiagonalOpts = " << iterDiagonalOpts <<
00524              " transpose = "  << transpose << " iterDist = " << iterDist << std::endl ; 
00525 
00526 
00527         if ( iterDiagonalOpts )  Comm.SetTracebackMode(1);  // Turn off positive Epetra warnings (i.e. iniefficient code, such as memory re-allocation)
00528 
00529       if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ; 
00530 
00531         RCP<Epetra_CrsMatrix> Bmat = NewMatNewMap( *Amat, 
00532                        iterDiagonalOpts, 
00533                        iterRowindex,
00534                        iterColindex,
00535                        iterRangemap,
00536                        iterDomainmap
00537                        ) ; 
00538       if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ; 
00539         Comm.SetTracebackMode(2);
00540 
00541         if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << 
00542              " filename = " << filename <<
00543              " symmetric = " << symmetric <<
00544              " distribute = " << distribute <<
00545              " iterRowindex = " << iterRowindex <<
00546              " iterColindex = " << iterColindex <<
00547              " iterRangemap = " << iterRangemap <<
00548              " iterDomainmap = " << iterDomainmap <<
00549              " EpetraMatrixType = " << EpetraMatrixType <<
00550              " iterDiagonalOpts = " << iterDiagonalOpts <<
00551              " transpose = "  << transpose << " iterDist = " << iterDist << std::endl ; 
00552 
00553         //
00554         //  This causes a failure in Amesos_Superludist:
00555         Epetra_CrsMatrix* Cmat = &*Bmat;
00556         //  Epetra_CrsMatrix* Cmat = Amat ;
00557      
00558 
00559         int Level ; 
00560         double MaxError ;
00561         if ( Rcond*Rcond1*Rcond2 > 1e-16 ) { 
00562           Level = EPETRA_MIN( 3, MaxLevel );
00563           MaxError = Rcond*Rcond1*Rcond2;
00564         } else if  ( Rcond*Rcond1 > 1e-16 ) {
00565           Level = EPETRA_MIN( 2, MaxLevel );
00566           MaxError = Rcond*Rcond1;
00567         } else {
00568           Level = EPETRA_MIN( 1, MaxLevel );
00569           MaxError = Rcond;
00570         }
00571 
00572         int NumTheseTests = 0 ; 
00573         if ( verbose ) {
00574           std::cout << " About to test  " << filename 
00575          << __FILE__ << "::"  << __LINE__
00576          << " EpetraMatrixType = " <<  EpetraMatrixType 
00577          << (transpose?" transpose":"" ) 
00578          << (distribute?" distribute":"" )
00579          << std::endl ; 
00580         }
00581         if ( iterDiagonalOpts == 0 ) 
00582           Comm.SetTracebackMode(2);
00583         else
00584           Comm.SetTracebackMode(1);  // In PerformOneSolveAndTest, MyMatWithDiag->ReplaceDiagonalValues may return 1 indicating that structurally non-zero elements were left untouched.
00585 
00586         int Error = TestAllClasses( AmesosClasses, EpetraMatrixType, 
00587             AmesosClassesInstalled, 
00588             Cmat, 
00589             transpose ,
00590             verbose, 
00591             symmetric, 
00592             Level,
00593             MaxError, 
00594             iterDiagonalOpts, 
00595             iterRowindex,
00596             iterColindex,
00597             iterRangemap,
00598             iterDomainmap,
00599             distribute,
00600             filename,
00601             error, 
00602             residual, 
00603             NumTheseTests ) ;
00604         NumTests += NumTheseTests ;
00605         NumErrors += Error ;
00606         if ( Comm.MyPID() == 0  && ( ( verbose && NumTheseTests ) || Error ) ) {
00607           std::cout << " Tested  " << filename 
00608          << __FILE__ << "::"  << __LINE__
00609          << " EpetraMatrixType = " <<  EpetraMatrixType 
00610          << (transpose?" transpose":"" ) 
00611          << (distribute?" distribute":"" ) << " error = " 
00612          << error 
00613          << " residual = " 
00614          << residual 
00615          << std::endl ; 
00616         }
00617         //      BUG:  Memory leak 
00618         //      delete &(Amat->RowMap()) ; 
00619         if ( Amat ) delete Amat ; 
00620         if ( readMap ) delete readMap ; 
00621 #if 0
00622         double relresidual = 
00623           errors[(int) AMESOS_SUPERLUDIST] = EPETRA_MAX( errors[ (int) AMESOS_SUPERLUDIST], error ) ; 
00624         residuals[(int) AMESOS_SUPERLUDIST] = EPETRA_MAX( residuals[ (int) AMESOS_SUPERLUDIST], residual ) ; 
00625         NumErrors += ( residual > maxresidual ) ; 
00626 #endif
00627         }
00628       }
00629       if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ; 
00630     }
00631       if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ; 
00632         }
00633       }
00634     }
00635   }
00636       }
00637     }
00638   }
00639 
00640   return NumErrors;
00641 } 
00642 
00643 #if 0
00644 #define TEST_P(variable) { { \
00645                       if ( true ) { std::cerr << "AMESOS_PRINT " << # variable << "= " << variable << std::endl; };  }\
00646                    }
00647 
00648 
00649 #define TEST_PRINT(variable) { { \
00650                       if ( true ) { std::cerr << "AMESOS_PRINT " << # variable  << "= " << variable <<  ", " \
00651                            << __FILE__ << ", line " << __LINE__ << std::endl; };  }\
00652                    }
00653 
00654 #endif
00655 
00656 //
00657 //  Usage:  TestOptions [-s] [-v] [-q]
00658 //
00659 //  -s = short
00660 //  -v = verbose
00661 //  -q = quiet 
00662 //
00663 
00664 int NextMain( int argc, char *argv[] ) {
00665 
00666 #ifdef EPETRA_MPI
00667   MPI_Init(&argc,&argv);
00668   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00669 #else
00670   Epetra_SerialComm Comm;
00671 #endif
00672 
00673 
00674   bool verbose = false; 
00675   bool small = false ; 
00676   bool quiet = false ; 
00677   if ( argc >= 2 && (argv[1][0] == '-') &&  (argv[1][1] == 'v') ) 
00678     verbose = true ; 
00679   if ( argc >= 3 && (argv[2][0] == '-') &&  (argv[2][1] == 'v') ) 
00680     verbose = true ; 
00681   if ( argc >= 4 && (argv[3][0] == '-') &&  (argv[3][1] == 'v') ) 
00682     verbose = true ; 
00683 
00684   if ( argc >= 2 && (argv[1][0] == '-') &&  (argv[1][1] == 's') ) 
00685     small = true ; 
00686   if ( argc >= 3 && (argv[2][0] == '-') &&  (argv[2][1] == 's') ) 
00687     small = true ; 
00688   if ( argc >= 4 && (argv[3][0] == '-') &&  (argv[3][1] == 's') ) 
00689     small = true ; 
00690 
00691   if ( argc >= 2 && (argv[1][0] == '-') &&  (argv[1][1] == 'q') ) 
00692     quiet = true ; 
00693   if ( argc >= 3 && (argv[2][0] == '-') &&  (argv[2][1] == 'q') ) 
00694     quiet = true ; 
00695   if ( argc >= 4 && (argv[3][0] == '-') &&  (argv[3][1] == 'q') ) 
00696     quiet = true ; 
00697 
00698 
00699   if ( argc >= 2 && (argv[1][0] == '-') &&  (argv[1][1] == 'h') ) {
00700     std::cerr << "Usage TestOptions [-s] [-v] [-q] " << std::endl ; 
00701     std::cerr << "-v:  verbose  " << std::endl ; 
00702     std::cerr << "-s:  small  " << std::endl ; 
00703     std::cerr << "-q:  quiet  " << std::endl ; 
00704     exit(-1);
00705   }
00706 
00707 
00708 
00709 #ifdef HAVE_AMESOS_KLU
00710   AmesosClasses.push_back( "Amesos_Klu" );
00711 #endif
00712 
00713 #if 1
00714 
00715 #ifdef HAVE_AMESOS_PARAKLETE
00716   AmesosClasses.push_back( "Amesos_Paraklete" );
00717 #endif
00718 
00719 
00720 
00721 
00722 
00723 #ifdef HAVE_AMESOS_PARDISO
00724   //  bug #1915  
00725   //  bug #1998 - Enabling Amesos_Pardiso causes Amesos_Klu to fail - strange but true
00726   //  AmesosClasses.push_back( "Amesos_Pardiso" );
00727 #endif
00728 #ifdef HAVE_AMESOS_SUPERLUDIST
00729   AmesosClasses.push_back( "Amesos_Superludist" );
00730 #endif
00731 
00732 
00733 
00734 
00735 #ifdef HAVE_AMESOS_LAPACK
00736   AmesosClasses.push_back( "Amesos_Lapack" );
00737 #endif
00738 
00739 #ifdef HAVE_AMESOS_SUPERLU
00740   AmesosClasses.push_back( "Amesos_Superlu" );
00741 #endif
00742 
00743 #ifdef HAVE_AMESOS_TAUCS
00744   AmesosClasses.push_back( "Amesos_Taucs" );
00745 #endif
00746 
00747 #ifdef HAVE_AMESOS_UMFPACK
00748   AmesosClasses.push_back( "Amesos_Umfpack" );
00749 #endif
00750 
00751 
00752 
00753 #ifdef HAVE_AMESOS_DSCPACK
00754   if ( ! quiet ) AmesosClasses.push_back( "Amesos_Dscpack" );         //  bug #1205 
00755 #endif
00756 
00757 #ifdef HAVE_AMESOS_MUMPS
00758   AmesosClasses.push_back( "Amesos_Mumps" );
00759 #endif
00760 
00761 #ifdef HAVE_AMESOS_SCALAPACK
00762   AmesosClasses.push_back( "Amesos_Scalapack" ) ;
00763 #endif
00764 
00765 
00766 
00767 #endif
00768 
00769   NumAmesosClasses = AmesosClasses.size();
00770   std::vector<bool> AmesosClassesInstalled( NumAmesosClasses );
00771 
00772   assert( NumAmesosClasses > 0 ) ; 
00773 
00774 
00775   if ( Comm.MyPID() != 0 ) verbose = false ; 
00776 #if 0
00777   //
00778   //  Wait for a character to allow time to attach the debugger
00779   //
00780   if ( Comm.MyPID() == 0 ) {
00781     char what = 'a'; 
00782     while ( what == 'a' )    // I don't know why we need this while loop  at least on bewoulf
00783       std::cin >> what ; 
00784   }
00785   Comm.Barrier();
00786 
00787   std::cout << __FILE__ << "::" << __LINE__ << " Comm.MyPID() = "  << Comm.MyPID() << std::endl ; 
00788 #endif
00789 
00790 
00791 
00792   Epetra_LinearProblem Problem;
00793   Amesos_BaseSolver* Abase ; 
00794   Amesos Afactory;
00795 
00796   Comm.SetTracebackMode(2);
00797 
00798 #ifndef HAVE_AMESOS_EPETRAEXT
00799     if ( ! quiet && Comm.MyPID() == 0 ) 
00800       std::cout << "Amesos has been built without epetraext, capabilites requiring epetraext, such as reindexing and Amesos_Paraklete non-transpose solves, will not be tested" << std::endl ; 
00801 #endif
00802 
00803   for (int i=0; i < NumAmesosClasses; i++ ) {
00804 
00805 
00806     Abase = Afactory.Create( &AmesosClasses[i][0], Problem ) ; 
00807     if ( Abase == 0 ) {
00808       if ( !quiet && Comm.MyPID() == 0  ) std::cout << AmesosClasses[i] << " not built in this configuration"  << std::endl ;
00809       AmesosClassesInstalled[i] = false;
00810     } else {
00811       if (  !quiet && Comm.MyPID() == 0  ) std::cout << " Testing " << AmesosClasses[i] << std::endl ;
00812       AmesosClassesInstalled[i] = true;
00813       Teuchos::ParameterList ParamList ;
00814       ParamList.set( "NoDestroy", true );    // Prevents Amesos_Mumps from deleting data
00815       Abase->SetParameters( ParamList );     // which causes Amesos_Mumps to crash  on this trivial instantiation 
00816     }
00817     delete Abase ; 
00818     }
00819 
00820   int result = 0 ; 
00821   int numtests = 0 ;
00822 
00823   //  ImpcolB.rua fails - the failure could be in the test code, in particular in NewMatNewMap.cpp
00824   //  result += TestOneMatrix( AmesosClassesInstalled, "../Test_Basic/ImpcolB.rua", Comm, verbose, false, 1e-9 , numtests ) ;
00825 
00826   //
00827   //  Khead.triS remains non-singular even after a diagaonal element is removed from the map 
00828   //
00829   // Khead.triS fails on DSCPACK 
00830   result += TestOneMatrix( AmesosClassesInstalled, (char *) "../Test_Basic/SuperLU.triU", Comm, verbose, false, 1e-6 , numtests ) ;
00831 
00832   //
00833   //  small is set by TestValgrind - keep testing to a minimum because execution time is so slow
00834   //  quiet is set by TestQuietAmesos - dscpack is not quiet at the moment, hence we can't test symmetric matrices
00835   //  in TestQuietAmesos
00836   //
00837   //  bug #1205 
00838   //
00839   if ( ! small ) {
00840   result += TestOneMatrix( AmesosClassesInstalled, (char *) "../Test_Basic/MissingADiagonal.mtx", Comm, verbose, false, 1e-2 , numtests ) ;
00841     result += TestOneMatrix( AmesosClassesInstalled, "../Test_Basic/Khead.triS", Comm, verbose, true, 1e-6 , numtests ) ;
00842     result += TestOneMatrix( AmesosClassesInstalled, (char *) "../Test_Basic/bcsstk04.mtx", Comm, verbose, false, 1e-4 , numtests ) ;
00843     //
00844     //  The file reader for .rua files is not quiet 
00845     //
00846     if (! quiet) { 
00847       result += TestOneMatrix( AmesosClassesInstalled, (char *) "../Test_Basic/Diagonal.mtx", Comm, verbose, false, 1e-1 , numtests ) ;
00848       if ( Comm.NumProc() == 1) {
00849   result += TestOneMatrix( AmesosClassesInstalled, (char *) "../Test_Basic/662_bus_out.rsa", Comm, verbose, false, 1e-5 , numtests ) ;
00850   result += TestOneMatrix( AmesosClassesInstalled, (char *) "../Test_Basic/SuperLU.rua", Comm, verbose, false, 1e-2 , numtests ) ;
00851   result += TestOneMatrix( AmesosClassesInstalled, "../Test_Basic/ImpcolB.rua", Comm, verbose, false, 1e-6 , numtests ) ;
00852   //  result += TestOneMatrix( AmesosClassesInstalled, "../Test_Basic/ImpcolB.rua", Comm, verbose, false, 1e-6 , numtests ) ;
00853       }
00854     }
00855   }
00856 
00857   // bug #2184 - Amesos_Klu fails to detect Structurally singular matrices 
00858   // This test, as modified in PerformOneSolveAndTest.cpp, tests the present
00859   // beahviour - i.e. that SymbolicFactorization() fails to detect the
00860   // structurally singular matrix, but that NumericFactorization() 
00861   // catches the singularity instead.  
00862   result+=TestErrors( AmesosClassesInstalled, (char *) "../Test_Basic/StructurallySingular.mtx",
00863           Comm, verbose, numtests ) ;
00864   result+=TestErrors( AmesosClassesInstalled, (char *) "../Test_Basic/NumericallySingular.mtx",
00865           Comm, verbose, numtests ) ;
00866  
00867   if ( ! quiet && Comm.MyPID() == 0 ) std::cout << result << " Tests failed "  << numtests << " Tests performed " << std::endl ; 
00868 
00869   if ( result == 0 && numtests > 0 ) {
00870     if (! quiet && Comm.MyPID() == 0)
00871       std::cout << std::endl << "TEST PASSED" << std::endl << std::endl;
00872   }
00873   else {
00874     if (Comm.MyPID() == 0)
00875       std::cout << std::endl << "TEST FAILED" << std::endl << std::endl;
00876     AMESOS_CHK_ERR( 1 ) ; 
00877   }
00878 
00879 #ifdef EPETRA_MPI
00880   MPI_Finalize();
00881 #endif
00882 
00883   return result ; 
00884 }
00885 
00886 //
00887 //  I put this in hoping that this would eliminate a bogus memory leak report 
00888 //  from valgrind. 
00889 //
00890 int main( int argc, char *argv[] ) {
00891   int retval = NextMain( argc, argv ) ; 
00892   return retval ;
00893 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines