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