Amesos Package Browser (Single Doxygen Collection) Development
Test_Detailed/cxx_main.cpp
Go to the documentation of this file.
00001 #include "Amesos_ConfigDefs.h"
00002 
00003 #ifdef HAVE_MPI
00004 #include "mpi.h"
00005 #include "Epetra_MpiComm.h"
00006 #else
00007 #include "Epetra_SerialComm.h"
00008 #endif
00009 #include "Epetra_Map.h"
00010 #include "Epetra_Vector.h"
00011 #include "Epetra_Time.h"
00012 #include "Amesos.h"
00013 #include "Amesos_BaseSolver.h"
00014 #include "Amesos_TestRowMatrix.h"
00015 #include "Teuchos_ParameterList.hpp"
00016 #include "Galeri_Maps.h"
00017 #include "Galeri_CrsMatrices.h"
00018 
00019 #ifdef HAVE_VALGRIND_H
00020 #include <valgrind/valgrind.h>
00021 #define HAVE_VALGRIND
00022 #else
00023 #ifdef HAVE_VALGRIND_VALGRIND_H
00024 #include <valgrind/valgrind.h>
00025 #define HAVE_VALGRIND
00026 #endif 
00027 #endif 
00028 
00029 using namespace Teuchos;
00030 using namespace Galeri;
00031 
00032 //=============================================================================
00033 
00034 bool TestTiming(const Amesos_BaseSolver* Solver,
00035                 const Epetra_Comm& Comm)
00036 {
00037   // Get the timings here
00038   Teuchos::ParameterList TimingsList;
00039   Solver->GetTiming( TimingsList );
00040   
00041   bool testPassed = true;
00042   
00043   try {
00044     // you can find out how much time was spent in ...
00045     double sfact_time, nfact_time, solve_time;
00046     double mtx_conv_time, mtx_redist_time, vec_redist_time;
00047     
00048     // 1) The symbolic factorization
00049     //    (parameter doesn't always exist)
00050     sfact_time = TimingsList.get( "Total symbolic factorization time", 0.0 );
00051     
00052     // 2) The numeric factorization
00053     //    (always exists if NumericFactorization() is called)
00054     nfact_time = Teuchos::getParameter<double>( TimingsList, "Total numeric factorization time" );
00055     
00056     // 3) Solving the linear system
00057     //    (always exists if Solve() is called)
00058     solve_time = Teuchos::getParameter<double>( TimingsList, "Total solve time" );
00059     
00060     // 4) Converting the matrix to the accepted format for the solver
00061     //    (always exists if SymbolicFactorization() is called)
00062     mtx_conv_time = Teuchos::getParameter<double>( TimingsList, "Total solve time" );
00063     
00064     // 5) Redistributing the matrix for each solve to the accepted format for the solver
00065     mtx_redist_time = TimingsList.get( "Total matrix redistribution time", 0.0 );
00066     
00067     // 6) Redistributing the vector for each solve to the accepted format for the solver
00068     vec_redist_time = TimingsList.get( "Total vector redistribution time", 0.0 );
00069   }
00070   catch( exception& e ) {
00071     if (Comm.MyPID() == 0)
00072       std::cout << std::endl << "Exception caught in TestTiming() : " << e.what() << std::endl;
00073     testPassed = false;
00074   }
00075   
00076   return testPassed;
00077 
00078 }
00079 
00080 //=============================================================================
00081 bool CheckError(const std::string SolverType,
00082     const std::string Descriptor,
00083     const Epetra_RowMatrix& A,
00084     const Epetra_MultiVector& x,
00085     const Epetra_MultiVector& b,
00086     const Epetra_MultiVector& x_exact)
00087 {
00088   if (A.Comm().MyPID() == 0)
00089     std::cout << std::endl << "*** " << SolverType << " : " 
00090          << Descriptor << " ***" << std::endl;
00091   std::vector<double> Norm;
00092   int NumVectors = x.NumVectors();
00093   Norm.resize(NumVectors);
00094   Epetra_MultiVector Ax(x);
00095   A.Multiply(false,x,Ax);
00096   Ax.Update(1.0,b,-1.0);
00097   Ax.Norm2(&Norm[0]);
00098   bool TestPassed = false;
00099   double TotalNorm = 0.0;
00100   for (int i = 0 ; i < NumVectors ; ++i) {
00101     TotalNorm += Norm[i];
00102   }
00103   if (A.Comm().MyPID() == 0)
00104     std::cout << "||Ax - b||  = " << TotalNorm << std::endl;
00105   if (TotalNorm < 1e-5 )
00106     TestPassed = true;
00107   else
00108     TestPassed = false;
00109 
00110   Ax.Update (1.0,x,-1.0,x_exact,0.0);
00111   Ax.Norm2(&Norm[0]);
00112   for (int i = 0 ; i < NumVectors ; ++i) {
00113     TotalNorm += Norm[i];
00114   }
00115   if (A.Comm().MyPID() == 0)
00116     std::cout << "||x - x_exact||  = " << TotalNorm << std::endl;
00117   if (TotalNorm < 1e-5 )
00118     TestPassed = true;
00119   else
00120     TestPassed = false;
00121 
00122   return(TestPassed);
00123 }
00124 
00125 //=============================================================================
00126 bool Test(char* SolverType,
00127     Epetra_RowMatrix& A, Epetra_MultiVector& x_A,
00128     Epetra_MultiVector& b_A, Epetra_MultiVector& x_exactA,
00129     Epetra_RowMatrix& B, Epetra_MultiVector& x_B,
00130     Epetra_MultiVector& b_B, Epetra_MultiVector& x_exactB,
00131     Epetra_RowMatrix& C, Epetra_MultiVector& x_C,
00132     Epetra_MultiVector& b_C, Epetra_MultiVector& x_exactC)
00133 {
00134   bool TestPassed = true;
00135   Epetra_LinearProblem ProblemA;
00136   Epetra_LinearProblem ProblemB;
00137   Epetra_LinearProblem ProblemC;
00138   Amesos Factory;
00139   Amesos_BaseSolver* Solver;
00140 
00141   // Test simple usage:
00142   // - set problem (empty)
00143   // - set matrix and vector
00144   // - call Solve() directly
00145   if (true) {
00146     x_A.PutScalar(0.0);
00147     ProblemA.SetOperator((Epetra_RowMatrix*)0);
00148     ProblemA.SetLHS((Epetra_MultiVector*)0);
00149     ProblemA.SetRHS((Epetra_MultiVector*)0);
00150 
00151     Solver = Factory.Create(SolverType,ProblemA);
00152   
00153     ProblemA.SetOperator(&A);
00154     ProblemA.SetLHS(&x_A);
00155     ProblemA.SetRHS(&b_A);
00156 
00157     std::string ST = SolverType ;    
00158     if (! ( ST == "Amesos_Superludist" ) ) {    // Kludge see bug #1141 and bug #1138
00159       AMESOS_CHK_ERR(Solver->Solve());
00160 
00161       TestPassed = TestPassed && 
00162   CheckError(SolverType, "Solve() only", A,x_A,b_A,x_exactA) &&
00163   TestTiming(Solver, A.Comm());
00164     }
00165     delete Solver; 
00166   }
00167 
00168   // Test almost simple usage:
00169   // - set problem (empty)
00170   // - set matrix and vector
00171   // - call NumericFactorization()
00172   // - call Solve() directly
00173   
00174   if (true) 
00175   {
00176     x_A.PutScalar(0.0);
00177     ProblemA.SetOperator((Epetra_RowMatrix*)0);
00178     ProblemA.SetLHS((Epetra_MultiVector*)0);
00179     ProblemA.SetRHS((Epetra_MultiVector*)0);
00180 
00181     Solver = Factory.Create(SolverType,ProblemA);
00182   
00183     ProblemA.SetOperator(&A);
00184     AMESOS_CHK_ERR(Solver->NumericFactorization());
00185     ProblemA.SetLHS(&x_A);
00186     ProblemA.SetRHS(&b_A);
00187     AMESOS_CHK_ERR(Solver->Solve());
00188 
00189     TestPassed = TestPassed && 
00190       CheckError(SolverType, "NumFact() + Solve()", A,x_A,b_A,x_exactA) &&
00191       TestTiming(Solver, A.Comm());
00192 
00193     delete Solver; 
00194   }
00195 
00196   // Test normal usage:
00197   // - set problem (empty)
00198   // - set matrix
00199   // - call SymbolicFactorization() several times
00200   // - call NumericFactorization() several times
00201   // - set vectors
00202   // - call Solve() several times
00203   // Repeated calls should *not* break the object. Data are supposed
00204   // to be re-computed as needed.
00205   
00206   if (true) 
00207   {
00208     x_A.PutScalar(0.0);
00209     ProblemA.SetOperator((Epetra_RowMatrix*)0);
00210     ProblemA.SetLHS((Epetra_MultiVector*)0);
00211     ProblemA.SetRHS((Epetra_MultiVector*)0);
00212 
00213     Solver = Factory.Create(SolverType,ProblemA);
00214   
00215     ProblemA.SetOperator(&A);
00216     AMESOS_CHK_ERR(Solver->SymbolicFactorization());
00217     AMESOS_CHK_ERR(Solver->SymbolicFactorization());
00218     AMESOS_CHK_ERR(Solver->SymbolicFactorization());
00219     AMESOS_CHK_ERR(Solver->NumericFactorization());
00220     AMESOS_CHK_ERR(Solver->NumericFactorization());
00221     AMESOS_CHK_ERR(Solver->NumericFactorization());
00222     AMESOS_CHK_ERR(Solver->NumericFactorization());
00223     ProblemA.SetLHS(&x_A);
00224     ProblemA.SetRHS(&b_A);
00225 
00226     AMESOS_CHK_ERR(Solver->Solve());
00227     AMESOS_CHK_ERR(Solver->Solve());
00228     AMESOS_CHK_ERR(Solver->Solve());
00229     AMESOS_CHK_ERR(Solver->Solve());
00230     AMESOS_CHK_ERR(Solver->Solve());
00231     AMESOS_CHK_ERR(Solver->Solve());
00232 
00233     TestPassed = TestPassed && 
00234       CheckError(SolverType, "SymFact() + NumFact() + Solve()", 
00235      A,x_A,b_A,x_exactA) &&
00236       TestTiming(Solver,A.Comm());
00237 
00238     delete Solver; 
00239   }
00240 
00241   // Test normal usage:
00242   // - set problem (empty)
00243   // - set matrix (as A)
00244   // - call SymbolicFactorization()
00245   // - set pointer to B
00246   // - call NumericFactorization()
00247   // - set vectors for B
00248   // - call Solve()
00249   
00250   if (true)
00251   {
00252     x_A.PutScalar(0.0);
00253     ProblemA.SetOperator((Epetra_RowMatrix*)0);
00254     ProblemA.SetLHS((Epetra_MultiVector*)0);
00255     ProblemA.SetRHS((Epetra_MultiVector*)0);
00256 
00257     Solver = Factory.Create(SolverType,ProblemA);
00258  
00259     ProblemA.SetOperator(&A);
00260     AMESOS_CHK_ERR(Solver->SymbolicFactorization());
00261     ProblemA.SetOperator(&B);
00262     AMESOS_CHK_ERR(Solver->NumericFactorization());
00263     ProblemA.SetLHS(&x_B);
00264     ProblemA.SetRHS(&b_B);
00265 
00266     AMESOS_CHK_ERR(Solver->Solve());
00267 
00268     TestPassed = TestPassed && 
00269       CheckError(SolverType, "Set A, solve B", B,x_B,b_B,x_exactB) &&
00270       TestTiming(Solver, A.Comm());
00271 
00272     delete Solver; 
00273   }
00274 
00275   // Construct Solver with filled ProblemA.
00276   // Then, stick pointers for different vectors.
00277   // a different map as well.
00278   
00279   if (true)
00280   {
00281     x_C.PutScalar(0.0);
00282     ProblemA.SetOperator(&A);
00283     ProblemA.SetLHS(&x_A);
00284     ProblemA.SetRHS(&b_A);
00285 
00286     Solver = Factory.Create(SolverType,ProblemA);
00287 
00288     ProblemA.SetOperator(&C);
00289 
00290     std::string ST = SolverType ; 
00291     if (! ( ST == "Amesos_Superludist" ) ) { // Kludge see bug #1141 and bug #1138
00292 
00293       AMESOS_CHK_ERR(Solver->SymbolicFactorization());
00294       AMESOS_CHK_ERR(Solver->NumericFactorization());
00295       
00296       ProblemA.SetLHS(&x_C);
00297       ProblemA.SetRHS(&b_C);
00298       AMESOS_CHK_ERR(Solver->Solve());
00299       
00300       TestPassed = TestPassed && 
00301   CheckError(SolverType, "Set A, Solve C", C,x_C,b_C,x_exactC) &&
00302   TestTiming(Solver, A.Comm());
00303     }
00304     delete Solver; 
00305   }
00306 
00307   // Construct Solver with filled ProblemA, call Solve().
00308   // Then, replace the pointers for matrix and vectors,
00309   // and call Solve() again.
00310   
00311   if (true)
00312   {
00313     x_C.PutScalar(0.0);
00314     ProblemA.SetOperator(&A);
00315     ProblemA.SetLHS(&x_A);
00316     ProblemA.SetRHS(&b_A);
00317 
00318     Solver = Factory.Create(SolverType,ProblemA);
00319 
00320     std::string ST = SolverType ; 
00321     if (! ( ST == "Amesos_Superludist" ) ) {   // bug #1141 and bug #1138
00322       AMESOS_CHK_ERR(Solver->Solve());
00323       
00324       ProblemA.SetOperator(&C);
00325       AMESOS_CHK_ERR(Solver->SymbolicFactorization());
00326       AMESOS_CHK_ERR(Solver->NumericFactorization());
00327       
00328       ProblemA.SetLHS(&x_C);
00329       ProblemA.SetRHS(&b_C);
00330       AMESOS_CHK_ERR(Solver->Solve());
00331       
00332       TestPassed = TestPassed && 
00333   CheckError(SolverType, "Solve A + Solve C", C,x_C,b_C,x_exactC) &&
00334   TestTiming(Solver, A.Comm());
00335     }
00336     delete Solver; 
00337   }
00338 
00339   return(TestPassed);
00340 }
00341 
00342 //=============================================================================
00343 
00344 int SubMain( Epetra_Comm &Comm ) {
00345   // Creation of data
00346   // A and B refer to two different matrices, with the *same*
00347   // structure and *different* values. Clearly, the map is the same.
00348   //
00349   // C refers to a completely different matrix
00350  
00351   int Size_AB = 20; // A and B have size Size_AB * Size_AB
00352   int Size_C = 30; 
00353   int NumVectors_AB = 7;
00354   int NumVectors_C = 13;
00355   
00356 #ifdef HAVE_VALGRIND 
00357   if ( RUNNING_ON_VALGRIND ) {
00358    Size_AB = 6; // must be square
00359    Size_C = 6; 
00360    NumVectors_AB = 2;
00361    NumVectors_C = 3;
00362   }
00363 #endif
00364 
00365   Teuchos::ParameterList GaleriList;
00366 
00367   GaleriList.set("n", Size_AB * Size_AB);
00368   GaleriList.set("nx", Size_AB);
00369   GaleriList.set("ny", Size_AB);
00370   Epetra_Map* Map_A = CreateMap("Interlaced", Comm, GaleriList);
00371   Epetra_CrsMatrix* Matrix_A = CreateCrsMatrix("Recirc2D", Map_A, GaleriList);
00372   Epetra_CrsMatrix* Matrix_B = CreateCrsMatrix("Laplace2D", Map_A, GaleriList);
00373 
00374   GaleriList.set("n", Size_C);
00375   Epetra_Map* Map_C = CreateMap("Interlaced", Comm, GaleriList);
00376   Epetra_CrsMatrix* Matrix_C = CreateCrsMatrix("Minij", Map_A, GaleriList);
00377 
00378   Amesos_TestRowMatrix A(Matrix_A);
00379   Amesos_TestRowMatrix B(Matrix_B);
00380   Amesos_TestRowMatrix C(Matrix_C);
00381 
00382   Epetra_MultiVector x_A(A.OperatorDomainMap(),NumVectors_AB);
00383   Epetra_MultiVector x_exactA(A.OperatorDomainMap(),NumVectors_AB);
00384   Epetra_MultiVector b_A(A.OperatorRangeMap(),NumVectors_AB);
00385   x_exactA.Random();
00386   A.Multiply(false,x_exactA,b_A);
00387 
00388   Epetra_MultiVector x_B(B.OperatorDomainMap(),NumVectors_AB);
00389   Epetra_MultiVector x_exactB(B.OperatorDomainMap(),NumVectors_AB);
00390   Epetra_MultiVector b_B(B.OperatorRangeMap(),NumVectors_AB);
00391   x_exactB.Random();
00392   B.Multiply(false,x_exactB,b_B);
00393 
00394   Epetra_MultiVector x_C(C.OperatorDomainMap(),NumVectors_C);
00395   Epetra_MultiVector x_exactC(C.OperatorDomainMap(),NumVectors_C);
00396   Epetra_MultiVector b_C(C.OperatorRangeMap(),NumVectors_C);
00397   x_exactC.Random();
00398   C.Multiply(false,x_exactC,b_C);
00399 
00400   std::vector<std::string> SolverType;
00401   SolverType.push_back("Amesos_Klu");
00402   SolverType.push_back("Amesos_Lapack");
00403   SolverType.push_back("Amesos_Umfpack");
00404   SolverType.push_back("Amesos_Superlu");
00405   SolverType.push_back("Amesos_Superludist");
00406   SolverType.push_back("Amesos_Mumps");
00407   // NOTE: DSCPACK does not support Epetra_RowMatrix's
00408 
00409   Amesos Factory;
00410   bool TestPassed = true;
00411 
00412   for (unsigned int i = 0 ; i < SolverType.size() ; ++i) {
00413     std::string Solver = SolverType[i];
00414     if (Factory.Query((char*)Solver.c_str())) {
00415       bool ok = Test((char*)Solver.c_str(),
00416          A, x_A, b_A, x_exactA,
00417          B, x_B, b_B, x_exactB,
00418          C, x_C, b_C, x_exactC);
00419       TestPassed = TestPassed && ok;
00420     }
00421     else
00422     {
00423       if (Comm.MyPID() == 0)
00424         std::cout << "Solver " << Solver << " not available" << std::endl;
00425     }
00426   }
00427 
00428   delete Matrix_A;
00429   delete Matrix_B;
00430   delete Matrix_C;
00431   delete Map_A;
00432   delete Map_C;
00433 
00434   if (TestPassed) {
00435     if (Comm.MyPID() == 0)
00436       std::cout << std::endl << "TEST PASSED" << std::endl << std::endl;
00437     return(EXIT_SUCCESS);
00438   }
00439   else {
00440     if (Comm.MyPID() == 0)
00441       std::cout << std::endl << "TEST FAILED" << std::endl << std::endl;
00442     // exit without calling MPI_Finalize() to raise an error
00443     exit(EXIT_FAILURE);
00444   }
00445 
00446 }
00447 
00448 int main(int argc, char *argv[]) {
00449 
00450 #ifdef HAVE_MPI
00451   MPI_Init(&argc, &argv);
00452   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00453 #else
00454   Epetra_SerialComm Comm;
00455 #endif
00456 
00457 #if 0
00458   if ( Comm.MyPID() == 0 ) {
00459     std::cout << "Enter a char to continue" ;
00460     char any;
00461     std::cin >> any ; 
00462   }
00463   Comm.Barrier();
00464 #endif
00465 
00466   int retval = SubMain( Comm ) ; 
00467 
00468 #ifdef HAVE_MPI
00469   MPI_Finalize();
00470 #endif
00471 
00472   return retval ; 
00473 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines