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