Amesos Package Browser (Single Doxygen Collection) Development
Test_Epetra_RowMatrix/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_Vector.h"
00010 #include "Epetra_CrsMatrix.h"
00011 #include "Amesos.h"
00012 #include "Amesos_TestRowMatrix.h"
00013 #include "Teuchos_ParameterList.hpp"
00014 #include "Galeri_Maps.h"
00015 #include "Galeri_CrsMatrices.h"
00016 
00017 using namespace Galeri;
00018 
00019 bool quiet = false ;
00020 
00021 // ====================================================================== 
00022 // this function tests two things:
00023 // - the return code from Amesos;
00024 // - the true residual.
00025 // The return value is `true' if both tests are passed, `false' otherwise.
00026 //
00027 // \author Marzio Sala, SNL 9214.
00028 //
00029 // \date Last updated on 21-Apr-04.
00030 // ====================================================================== 
00031 
00032 bool TestAmesos(char ProblemType[], Teuchos::ParameterList& AmesosList,
00033                 bool UseTranspose, Epetra_RowMatrix* A, Epetra_MultiVector* lhs,
00034                 Epetra_MultiVector* rhs)
00035 {
00036   Epetra_LinearProblem Problem;
00037   Amesos A_Factory;
00038 
00039   Amesos_BaseSolver * Solver = A_Factory.Create(ProblemType, Problem);
00040   assert (Solver != 0);
00041 
00042   // Both sentences should work
00043   Solver->SetUseTranspose(UseTranspose);
00044 
00045   Solver->SetParameters(AmesosList);
00046 
00047   // create a rhs corresponding to lhs or 1's
00048   lhs->PutScalar(1.0);
00049   A->Multiply(UseTranspose,*lhs,*rhs);
00050   lhs->PutScalar(0.0);
00051 
00052   Problem.SetOperator(A);
00053 
00054   if (Solver->SymbolicFactorization()) return(false);
00055   if (Solver->NumericFactorization()) return(false);
00056 
00057   // set sol and rhs here
00058   Problem.SetLHS(lhs);
00059   Problem.SetRHS(rhs);
00060 
00061   if (Solver->Solve()) return(false);
00062 
00063   // compute difference between exact solution and Amesos
00064   double d = 0.0, d_tot = 0.0;
00065 
00066   for (int i=0 ; i<lhs->Map().NumMyElements() ; ++i)
00067     for (int j = 0 ; j < lhs->NumVectors() ; ++j)
00068       d += ((*lhs)[j][i] - 1.0) * ((*lhs)[j][i] - 1.0);
00069 
00070   A->Comm().SumAll(&d,&d_tot,1);
00071 
00072   // compute ||Ax - b||
00073   std::vector<double> Norm(rhs->NumVectors());
00074 
00075 
00076   Epetra_MultiVector Ax(*rhs);
00077   A->Multiply(UseTranspose, *lhs, Ax);
00078   Ax.Update(1.0, *rhs, -1.0);
00079   Ax.Norm2(&Norm[0]);
00080 
00081   std::string msg = ProblemType;
00082 
00083   if (!quiet && !A->Comm().MyPID()) 
00084   {
00085     std::cout << std::endl;
00086     std::cout << msg << " : Using " << A->Comm().NumProc() << " processes, UseTranspose = " << UseTranspose << std::endl;
00087     for (int j = 0 ; j < rhs->NumVectors() ; ++j)
00088       std::cout << msg << " : eq " << j 
00089      << ", ||A x - b||_2 = " << Norm[j] << std::endl;
00090     std::cout << msg << " : ||x_exact - x||_2 = " << sqrt(d_tot) << std::endl;
00091   }
00092 
00093   if (Norm[0] > 1e-9)  
00094   {
00095     std::cerr << std::endl << msg << " WARNING : TEST FAILED!" << std::endl;
00096     return(false);
00097   }
00098 
00099   delete Solver;
00100 
00101   return(true);
00102 }
00103 
00104 void driver(Epetra_Comm& Comm, const bool IsSymmetric, const bool UseTranspose, 
00105             std::vector<std::string>& SolverType)
00106 {
00107   std::string ProblemType;
00108   if (IsSymmetric)
00109     ProblemType = "Laplace2D";
00110   else
00111     ProblemType = "Recirc2D";
00112 
00113   Teuchos::ParameterList GaleriList;
00114   GaleriList.set("nx", 8);
00115   GaleriList.set("ny", 8 * Comm.NumProc());
00116   GaleriList.set("mx", 1);
00117   GaleriList.set("my", Comm.NumProc());
00118 
00119   Epetra_Map* Map = CreateMap("Cartesian2D", Comm, GaleriList);
00120   Epetra_CrsMatrix* Matrix = CreateCrsMatrix(ProblemType, Map, GaleriList);
00121   Epetra_MultiVector LHS(*Map, 3); LHS.PutScalar(0.0);
00122   Epetra_MultiVector RHS(*Map, 3); RHS.Random();
00123 
00124   Amesos_TestRowMatrix A(Matrix);
00125 
00126   Amesos Factory;  
00127   
00128   bool res;
00129 
00130   // If a given test fails, than the code stops, bue to the assert()
00131   // statement.
00132   for (unsigned int i = 0 ; i < SolverType.size() ; ++i) 
00133   {
00134     std::string Solver = SolverType[i];
00135 
00136     if (Factory.Query((char*)Solver.c_str())) 
00137     {
00138       // solve with matrix
00139       Teuchos::ParameterList AmesosList;
00140       AmesosList.set("Redistribute",true);
00141       res = TestAmesos((char*)Solver.c_str(), AmesosList, false, 
00142                        &A, &LHS, &RHS);
00143       assert (res == true);
00144       if (UseTranspose) {
00145         // solve transpose with matrix
00146         Teuchos::ParameterList AmesosList;
00147         res  = TestAmesos((char*)Solver.c_str(), AmesosList, true, 
00148                           &A, &LHS, &RHS);
00149         assert (res == true);
00150       }
00151     } 
00152     else
00153       if (!quiet && !Comm.MyPID()) 
00154       {
00155         std::cerr << std::endl;
00156         std::cerr << "WARNING: SOLVER `" << Solver << "' NOT TESTED" << std::endl;
00157         std::cerr << std::endl;
00158       }
00159   }
00160 
00161   delete Matrix;
00162   delete Map;
00163 }
00164 
00165 // =========== //
00166 // main driver //
00167 // =========== //
00168 //
00169 int main(int argc, char *argv[]) 
00170 {
00171 #ifdef HAVE_MPI
00172   MPI_Init(&argc, &argv);
00173   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00174 #else
00175   Epetra_SerialComm Comm;
00176 #endif
00177 
00178   if ( argc > 1 && argv[1][0] == '-' &&  argv[1][1] == 'q' ) quiet = true ;
00179 
00180   // NOTE: DSCPACK does not support RowMatrix's.
00181   
00182   if (true)
00183   {
00184     // non-symmetric matrix, test A and A^T
00185     std::vector<std::string> SolverType;
00186     SolverType.push_back("Amesos_Lapack");
00187     SolverType.push_back("Amesos_Klu");
00188     SolverType.push_back("Amesos_Umfpack");
00189     SolverType.push_back("Amesos_Superlu");
00190 //    SolverType.push_back("Amesos_Mumps");       Bug #1896 
00191     SolverType.push_back("Amesos_Scalapack");
00192     driver(Comm, false, true, SolverType);
00193   }
00194  
00195   if (true)
00196   {
00197     // non-symmetric matrix, test only A
00198     std::vector<std::string> SolverType;
00199     //  
00200     //    SolverType.push_back("Amesos_Pardiso");    bug #1994
00201     SolverType.push_back("Amesos_Superludist");
00202     driver(Comm, false, false, SolverType);
00203   }
00204 
00205   // I have some problems with Taucs with LAM
00206   if (true)
00207   {
00208     // symmetric
00209     std::vector<std::string> SolverType;
00210     SolverType.push_back("Amesos_Taucs");
00211     driver(Comm, true, false, SolverType);
00212   }
00213 
00214 #ifdef HAVE_MPI
00215   MPI_Finalize();
00216 #endif
00217   return(0); 
00218 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines