Amesos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Test_Epetra_RowMatrix/cxx_main.cpp
Go to the documentation of this file.
1 #include "Amesos_ConfigDefs.h"
2 
3 #ifdef HAVE_MPI
4 #include "mpi.h"
5 #include "Epetra_MpiComm.h"
6 #else
7 #include "Epetra_SerialComm.h"
8 #endif
9 #include "Epetra_Vector.h"
10 #include "Epetra_CrsMatrix.h"
11 #include "Amesos.h"
12 #include "Amesos_TestRowMatrix.h"
14 #include "Galeri_Maps.h"
15 #include "Galeri_CrsMatrices.h"
16 
17 using namespace Galeri;
18 
19 bool quiet = false ;
20 
21 // ======================================================================
22 // this function tests two things:
23 // - the return code from Amesos;
24 // - the true residual.
25 // The return value is `true' if both tests are passed, `false' otherwise.
26 //
27 // \author Marzio Sala, SNL 9214.
28 //
29 // \date Last updated on 21-Apr-04.
30 // ======================================================================
31 
32 bool TestAmesos(char ProblemType[], Teuchos::ParameterList& AmesosList,
33  bool UseTranspose, Epetra_RowMatrix* A, Epetra_MultiVector* lhs,
34  Epetra_MultiVector* rhs)
35 {
36  Epetra_LinearProblem Problem;
37  Amesos A_Factory;
38 
39  Amesos_BaseSolver * Solver = A_Factory.Create(ProblemType, Problem);
40  assert (Solver != 0);
41 
42  // Both sentences should work
43  Solver->SetUseTranspose(UseTranspose);
44 
45  Solver->SetParameters(AmesosList);
46 
47  // create a rhs corresponding to lhs or 1's
48  lhs->PutScalar(1.0);
49  A->Multiply(UseTranspose,*lhs,*rhs);
50  lhs->PutScalar(0.0);
51 
52  Problem.SetOperator(A);
53 
54  if (Solver->SymbolicFactorization()) return(false);
55  if (Solver->NumericFactorization()) return(false);
56 
57  // set sol and rhs here
58  Problem.SetLHS(lhs);
59  Problem.SetRHS(rhs);
60 
61  if (Solver->Solve()) return(false);
62 
63  // compute difference between exact solution and Amesos
64  double d = 0.0, d_tot = 0.0;
65 
66  for (int i=0 ; i<lhs->Map().NumMyElements() ; ++i)
67  for (int j = 0 ; j < lhs->NumVectors() ; ++j)
68  d += ((*lhs)[j][i] - 1.0) * ((*lhs)[j][i] - 1.0);
69 
70  A->Comm().SumAll(&d,&d_tot,1);
71 
72  // compute ||Ax - b||
73  std::vector<double> Norm(rhs->NumVectors());
74 
75 
76  Epetra_MultiVector Ax(*rhs);
77  A->Multiply(UseTranspose, *lhs, Ax);
78  Ax.Update(1.0, *rhs, -1.0);
79  Ax.Norm2(&Norm[0]);
80 
81  std::string msg = ProblemType;
82 
83  if (!quiet && !A->Comm().MyPID())
84  {
85  std::cout << std::endl;
86  std::cout << msg << " : Using " << A->Comm().NumProc() << " processes, UseTranspose = " << UseTranspose << std::endl;
87  for (int j = 0 ; j < rhs->NumVectors() ; ++j)
88  std::cout << msg << " : eq " << j
89  << ", ||A x - b||_2 = " << Norm[j] << std::endl;
90  std::cout << msg << " : ||x_exact - x||_2 = " << sqrt(d_tot) << std::endl;
91  }
92 
93  if (Norm[0] > 1e-9)
94  {
95  std::cerr << std::endl << msg << " WARNING : TEST FAILED!" << std::endl;
96  return(false);
97  }
98 
99  delete Solver;
100 
101  return(true);
102 }
103 
104 void driver(Epetra_Comm& Comm, const bool IsSymmetric, const bool UseTranspose,
105  std::vector<std::string>& SolverType)
106 {
107  std::string ProblemType;
108  if (IsSymmetric)
109  ProblemType = "Laplace2D";
110  else
111  ProblemType = "Recirc2D";
112 
113  Teuchos::ParameterList GaleriList;
114  GaleriList.set("nx", 8);
115  GaleriList.set("ny", 8 * Comm.NumProc());
116  GaleriList.set("mx", 1);
117  GaleriList.set("my", Comm.NumProc());
118 
119  Epetra_Map* Map = CreateMap("Cartesian2D", Comm, GaleriList);
120  Epetra_CrsMatrix* Matrix = CreateCrsMatrix(ProblemType, Map, GaleriList);
121  Epetra_MultiVector LHS(*Map, 3); LHS.PutScalar(0.0);
122  Epetra_MultiVector RHS(*Map, 3); RHS.Random();
123 
124  Amesos_TestRowMatrix A(Matrix);
125 
126  Amesos Factory;
127 
128  bool res;
129 
130  // If a given test fails, than the code stops, bue to the assert()
131  // statement.
132  for (unsigned int i = 0 ; i < SolverType.size() ; ++i)
133  {
134  std::string Solver = SolverType[i];
135 
136  if (Factory.Query((char*)Solver.c_str()))
137  {
138  // solve with matrix
139  Teuchos::ParameterList AmesosList;
140  AmesosList.set("Redistribute",true);
141  res = TestAmesos((char*)Solver.c_str(), AmesosList, false,
142  &A, &LHS, &RHS);
143  assert (res == true);
144  if (UseTranspose) {
145  // solve transpose with matrix
146  Teuchos::ParameterList AmesosList;
147  res = TestAmesos((char*)Solver.c_str(), AmesosList, true,
148  &A, &LHS, &RHS);
149  assert (res == true);
150  }
151  }
152  else
153  if (!quiet && !Comm.MyPID())
154  {
155  std::cerr << std::endl;
156  std::cerr << "WARNING: SOLVER `" << Solver << "' NOT TESTED" << std::endl;
157  std::cerr << std::endl;
158  }
159  }
160 
161  delete Matrix;
162  delete Map;
163 }
164 
165 // =========== //
166 // main driver //
167 // =========== //
168 //
169 int main(int argc, char *argv[])
170 {
171 #ifdef HAVE_MPI
172  MPI_Init(&argc, &argv);
173  Epetra_MpiComm Comm(MPI_COMM_WORLD);
174 #else
175  Epetra_SerialComm Comm;
176 #endif
177 
178  if ( argc > 1 && argv[1][0] == '-' && argv[1][1] == 'q' ) quiet = true ;
179 
180  // NOTE: DSCPACK does not support RowMatrix's.
181 
182  if (true)
183  {
184  // non-symmetric matrix, test A and A^T
185  std::vector<std::string> SolverType;
186  SolverType.push_back("Amesos_Lapack");
187  SolverType.push_back("Amesos_Klu");
188  SolverType.push_back("Amesos_Umfpack");
189  SolverType.push_back("Amesos_Superlu");
190 // SolverType.push_back("Amesos_Mumps"); Bug #1896
191  SolverType.push_back("Amesos_Scalapack");
192  driver(Comm, false, true, SolverType);
193  }
194 
195  if (true)
196  {
197  // non-symmetric matrix, test only A
198  std::vector<std::string> SolverType;
199  //
200  // SolverType.push_back("Amesos_Pardiso"); bug #1994
201  SolverType.push_back("Amesos_Superludist");
202  driver(Comm, false, false, SolverType);
203  }
204 
205  // I have some problems with Taucs with LAM
206  if (true)
207  {
208  // symmetric
209  std::vector<std::string> SolverType;
210  SolverType.push_back("Amesos_Taucs");
211  driver(Comm, true, false, SolverType);
212  }
213 
214 #ifdef HAVE_MPI
215  MPI_Finalize();
216 #endif
217  return(0);
218 }
void driver(Epetra_Comm &Comm, const bool IsSymmetric, const bool UseTranspose, std::vector< std::string > &SolverType)
void SetLHS(Epetra_MultiVector *X)
void SetOperator(Epetra_RowMatrix *A)
bool TestAmesos(char ProblemType[], Teuchos::ParameterList &AmesosList, bool UseTranspose, Epetra_RowMatrix *A, Epetra_MultiVector *lhs, Epetra_MultiVector *rhs)
virtual int Solve()=0
Solves A X = B (or AT x = B)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
virtual int NumericFactorization()=0
Performs NumericFactorization on the matrix A.
virtual int SymbolicFactorization()=0
Performs SymbolicFactorization on the matrix A.
virtual int SetParameters(Teuchos::ParameterList &ParameterList)=0
Updates internal variables.
Amesos_TestRowMatrix: a class to test Epetra_RowMatrix based codes.
virtual int MyPID() const =0
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
int NumMyElements() const
virtual const Epetra_Comm & Comm() const =0
virtual const Epetra_BlockMap & Map() const =0
int main(int argc, char *argv[])
Factory for binding a third party direct solver to an Epetra_LinearProblem.
Definition: Amesos.h:44
virtual int SetUseTranspose(bool UseTranspose)=0
If set true, X will be set to the solution of AT X = B (not A X = B)
int CreateCrsMatrix(const char *in_filename, const Epetra_Comm &Comm, Epetra_Map *&readMap, const bool transpose, const bool distribute, bool &symmetric, Epetra_CrsMatrix *&Matrix)
Definition: TestOptions.cpp:81
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
void SetRHS(Epetra_MultiVector *B)
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.
Definition: Amesos.cpp:69
virtual int NumProc() const =0
Amesos_BaseSolver: A pure virtual class for direct solution of real-valued double-precision operators...
bool Query(const char *ClassType)
Queries whether a given interface is available or not.
Definition: Amesos.cpp:193