Amesos Package Browser (Single Doxygen Collection) Development
Test_MultipleSolves/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_MultiVector.h"
00011 #include "Epetra_Time.h"
00012 #include "Epetra_CrsMatrix.h"
00013 #include "Amesos.h"
00014 #include "Teuchos_ParameterList.hpp"
00015 #include "Galeri_Maps.h"
00016 #include "Galeri_CrsMatrices.h"
00017 
00018 using namespace Galeri;
00019 
00020 bool TestAmesos(char ProblemType[],
00021     Epetra_RowMatrix& A,
00022     int NumVectors)
00023 {
00024 
00025   const Epetra_BlockMap& Map = A.OperatorDomainMap();
00026 
00027   Epetra_MultiVector x2(Map,NumVectors);
00028   Epetra_MultiVector x1(Map,NumVectors);
00029   Epetra_MultiVector x(Map,NumVectors);
00030   Epetra_MultiVector b(Map,NumVectors);
00031   Epetra_MultiVector residual(Map,NumVectors);
00032   Epetra_MultiVector temp(Map,NumVectors);
00033   
00034   Teuchos::ParameterList ParamList ;
00035   Epetra_LinearProblem Problem;
00036   Amesos_BaseSolver* Abase ; 
00037   Amesos Afactory;
00038 
00039   // Note that Abase is created with an empty Problem, none of A, x or b
00040   // have been specified at this point.  
00041   Abase = Afactory.Create(ProblemType,Problem ); 
00042   assert (Abase != 0);
00043 
00044   //  Factor A
00045   Problem.SetOperator(&A);
00046   EPETRA_CHK_ERR(Abase->SymbolicFactorization()); 
00047   EPETRA_CHK_ERR(Abase->NumericFactorization()); 
00048 
00049   b.Random();
00050 
00051   //  Solve Ax = b 
00052   //
00053   Problem.SetLHS(&x);
00054   Problem.SetRHS(&b);
00055   EPETRA_CHK_ERR(Abase->Solve()); 
00056 
00057   //  Solve Ax1 = x 
00058   //
00059   Problem.SetLHS(&x1);
00060   Problem.SetRHS(&x);
00061   EPETRA_CHK_ERR(Abase->Solve()); 
00062 
00063   //  Solve Ax2 = x1
00064   //
00065   Problem.SetLHS(&x2);
00066   Problem.SetRHS(&x1);
00067   EPETRA_CHK_ERR(Abase->Solve()); 
00068 
00069   //  Compute the residual: A^3 x2 - b
00070   //
00071   A.Multiply(false,x2, temp) ; //  temp = A x2
00072   A.Multiply(false,temp, x2) ; //  x2 = A^2 x2
00073   A.Multiply(false,x2, temp) ; //  temp = A^3 x2
00074   residual.Update( 1.0, temp, -1.0, b, 0.0 ) ;
00075 
00076   double norm_residual ;
00077   residual.Norm2( &norm_residual ) ; 
00078 
00079   if (A.Comm().MyPID() == 0) {
00080     std::cout << "norm2(A^3 x-b) = " << norm_residual << std::endl ; 
00081   }
00082 
00083   delete Abase;
00084 
00085   if (norm_residual < (1e-5))
00086     return(true);
00087   else
00088     return(false);
00089   
00090 }
00091 
00092 int main(int argc, char *argv[]) {
00093 
00094 #ifdef HAVE_MPI
00095   MPI_Init(&argc, &argv);
00096   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00097 #else
00098   Epetra_SerialComm Comm;
00099 #endif
00100 
00101   int NumVectors = 2;
00102 
00103   Teuchos::ParameterList GaleriList;
00104   GaleriList.set("nx", 8);
00105   GaleriList.set("ny", 8 * Comm.NumProc());
00106   GaleriList.set("mx", 1);
00107   GaleriList.set("my", Comm.NumProc());
00108 
00109   Epetra_Map* Map = CreateMap("Cartesian2D", Comm, GaleriList);
00110   Epetra_CrsMatrix* A = CreateCrsMatrix("Laplace2D", Map, GaleriList);
00111 
00112   Amesos Factory;  
00113   
00114   std::vector<std::string> SolverType;
00115   //  SolverType.push_back("Amesos_Lapack");
00116   SolverType.push_back("Amesos_Klu");
00117   SolverType.push_back("Amesos_Umfpack");
00118   SolverType.push_back("Amesos_Pardiso");
00119   SolverType.push_back("Amesos_Taucs");
00120   SolverType.push_back("Amesos_Superlu");
00121   SolverType.push_back("Amesos_Superludist");
00122   SolverType.push_back("Amesos_Mumps");
00123   SolverType.push_back("Amesos_Dscpack");
00124 
00125   bool TestPassed = true;
00126 
00127   for (unsigned int i = 0 ; i < SolverType.size() ; ++i) {
00128     std::string Solver = SolverType[i];
00129         std::cout  << Solver << " next  " << std::endl;
00130     if (Factory.Query((char*)Solver.c_str())) {
00131       if (Comm.MyPID() == 0)
00132         std::cout << "Testing " << Solver << std::endl;
00133       if(TestAmesos((char*)Solver.c_str(), *A, NumVectors) == false) {
00134         std::cout  << Solver << " Failed " << std::endl;
00135   TestPassed = false;
00136       } else { 
00137         std::cout  << Solver << " Passed " << std::endl;
00138       } 
00139     } else
00140       if (Comm.MyPID() == 0) {
00141   std::cerr << std::endl;
00142   std::cerr << "WARNING: SOLVER `" << Solver << "' NOT TESTED" << std::endl;
00143   std::cerr << std::endl;
00144       }
00145   }
00146    
00147   delete A;
00148   delete Map;
00149 
00150 #ifdef HAVE_MPI
00151   MPI_Finalize();
00152 #endif
00153 
00154   if (TestPassed) {
00155     if (Comm.MyPID() == 0) 
00156       std::cout << "TESTS PASSED!" << std::endl;
00157     return( EXIT_SUCCESS );
00158   } 
00159   else {
00160     if (Comm.MyPID() == 0) 
00161       std::cout << "TESTS FAILED!" << std::endl;
00162     return( EXIT_FAILURE );
00163   }
00164 
00165 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines