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