5 #include "Epetra_MpiComm.h" 7 #include "Epetra_SerialComm.h" 9 #include "Epetra_Map.h" 10 #include "Epetra_Vector.h" 11 #include "Epetra_Time.h" 16 #include "Galeri_Maps.h" 17 #include "Galeri_CrsMatrices.h" 19 #ifdef HAVE_VALGRIND_H 20 #include <valgrind/valgrind.h> 23 #ifdef HAVE_VALGRIND_VALGRIND_H 24 #include <valgrind/valgrind.h> 35 const Epetra_Comm&
Comm)
39 Solver->GetTiming( TimingsList );
41 bool testPassed =
true;
54 TimingsList.
get(
"Total symbolic factorization time", 0.0 );
59 Teuchos::getParameter<double>( TimingsList,
"Total numeric factorization time" );
64 Teuchos::getParameter<double>( TimingsList,
"Total solve time" );
69 Teuchos::getParameter<double>( TimingsList,
"Total solve time" );
73 TimingsList.get(
"Total matrix redistribution time", 0.0 );
77 TimingsList.get(
"Total vector redistribution time", 0.0 );
80 catch( std::exception& e ) {
81 if (
Comm.MyPID() == 0)
82 std::cout << std::endl <<
"Exception caught in TestTiming() : " << e.what() << std::endl;
92 const std::string Descriptor,
93 const Epetra_RowMatrix&
A,
94 const Epetra_MultiVector& x,
95 const Epetra_MultiVector& b,
96 const Epetra_MultiVector& x_exact)
98 if (
A.Comm().MyPID() == 0)
99 std::cout << std::endl <<
"*** " << SolverType <<
" : " 100 << Descriptor <<
" ***" << std::endl;
101 std::vector<double> Norm;
102 int NumVectors = x.NumVectors();
103 Norm.resize(NumVectors);
104 Epetra_MultiVector Ax(x);
105 A.Multiply(
false,x,Ax);
106 Ax.Update(1.0,b,-1.0);
108 bool TestPassed =
false;
109 double TotalNorm = 0.0;
110 for (
int i = 0 ; i < NumVectors ; ++i) {
111 TotalNorm += Norm[i];
113 if (
A.Comm().MyPID() == 0)
114 std::cout <<
"||Ax - b|| = " << TotalNorm << std::endl;
115 if (TotalNorm < 1e-5 )
120 Ax.Update (1.0,x,-1.0,x_exact,0.0);
122 for (
int i = 0 ; i < NumVectors ; ++i) {
123 TotalNorm += Norm[i];
125 if (
A.Comm().MyPID() == 0)
126 std::cout <<
"||x - x_exact|| = " << TotalNorm << std::endl;
127 if (TotalNorm < 1e-5 )
137 Epetra_RowMatrix&
A, Epetra_MultiVector& x_A,
138 Epetra_MultiVector& b_A, Epetra_MultiVector& x_exactA,
139 Epetra_RowMatrix&
B, Epetra_MultiVector& x_B,
140 Epetra_MultiVector& b_B, Epetra_MultiVector& x_exactB,
141 Epetra_RowMatrix&
C, Epetra_MultiVector& x_C,
142 Epetra_MultiVector& b_C, Epetra_MultiVector& x_exactC)
144 bool TestPassed =
true;
145 Epetra_LinearProblem ProblemA;
146 Epetra_LinearProblem ProblemB;
147 Epetra_LinearProblem ProblemC;
157 ProblemA.SetOperator((Epetra_RowMatrix*)0);
158 ProblemA.SetLHS((Epetra_MultiVector*)0);
159 ProblemA.SetRHS((Epetra_MultiVector*)0);
163 ProblemA.SetOperator(&
A);
164 ProblemA.SetLHS(&x_A);
165 ProblemA.SetRHS(&b_A);
167 std::string ST = SolverType ;
168 if (! ( ST ==
"Amesos_Superludist" ) ) {
171 TestPassed = TestPassed &&
172 CheckError(SolverType,
"Solve() only",
A,x_A,b_A,x_exactA) &&
187 ProblemA.SetOperator((Epetra_RowMatrix*)0);
188 ProblemA.SetLHS((Epetra_MultiVector*)0);
189 ProblemA.SetRHS((Epetra_MultiVector*)0);
193 ProblemA.SetOperator(&
A);
195 ProblemA.SetLHS(&x_A);
196 ProblemA.SetRHS(&b_A);
199 TestPassed = TestPassed &&
200 CheckError(SolverType,
"NumFact() + Solve()",
A,x_A,b_A,x_exactA) &&
219 ProblemA.SetOperator((Epetra_RowMatrix*)0);
220 ProblemA.SetLHS((Epetra_MultiVector*)0);
221 ProblemA.SetRHS((Epetra_MultiVector*)0);
225 ProblemA.SetOperator(&
A);
233 ProblemA.SetLHS(&x_A);
234 ProblemA.SetRHS(&b_A);
243 TestPassed = TestPassed &&
244 CheckError(SolverType,
"SymFact() + NumFact() + Solve()",
245 A,x_A,b_A,x_exactA) &&
263 ProblemA.SetOperator((Epetra_RowMatrix*)0);
264 ProblemA.SetLHS((Epetra_MultiVector*)0);
265 ProblemA.SetRHS((Epetra_MultiVector*)0);
269 ProblemA.SetOperator(&
A);
271 ProblemA.SetOperator(&
B);
273 ProblemA.SetLHS(&x_B);
274 ProblemA.SetRHS(&b_B);
278 TestPassed = TestPassed &&
279 CheckError(SolverType,
"Set A, solve B",
B,x_B,b_B,x_exactB) &&
292 ProblemA.SetOperator(&
A);
293 ProblemA.SetLHS(&x_A);
294 ProblemA.SetRHS(&b_A);
298 ProblemA.SetOperator(&
C);
300 std::string ST = SolverType ;
301 if (! ( ST ==
"Amesos_Superludist" ) ) {
306 ProblemA.SetLHS(&x_C);
307 ProblemA.SetRHS(&b_C);
310 TestPassed = TestPassed &&
311 CheckError(SolverType,
"Set A, Solve C",
C,x_C,b_C,x_exactC) &&
324 ProblemA.SetOperator(&
A);
325 ProblemA.SetLHS(&x_A);
326 ProblemA.SetRHS(&b_A);
330 std::string ST = SolverType ;
331 if (! ( ST ==
"Amesos_Superludist" ) ) {
334 ProblemA.SetOperator(&
C);
338 ProblemA.SetLHS(&x_C);
339 ProblemA.SetRHS(&b_C);
342 TestPassed = TestPassed &&
343 CheckError(SolverType,
"Solve A + Solve C",
C,x_C,b_C,x_exactC) &&
363 int NumVectors_AB = 7;
364 int NumVectors_C = 13;
367 if ( RUNNING_ON_VALGRIND ) {
377 GaleriList.
set(
"n", Size_AB * Size_AB);
378 GaleriList.
set(
"nx", Size_AB);
379 GaleriList.
set(
"ny", Size_AB);
380 Epetra_Map* Map_A = CreateMap(
"Interlaced",
Comm, GaleriList);
381 Epetra_CrsMatrix* Matrix_A =
CreateCrsMatrix(
"Recirc2D", Map_A, GaleriList);
382 Epetra_CrsMatrix* Matrix_B =
CreateCrsMatrix(
"Laplace2D", Map_A, GaleriList);
384 GaleriList.
set(
"n", Size_C);
385 Epetra_Map* Map_C = CreateMap(
"Interlaced",
Comm, GaleriList);
386 Epetra_CrsMatrix* Matrix_C =
CreateCrsMatrix(
"Minij", Map_A, GaleriList);
392 Epetra_MultiVector x_A(
A.OperatorDomainMap(),NumVectors_AB);
393 Epetra_MultiVector x_exactA(
A.OperatorDomainMap(),NumVectors_AB);
394 Epetra_MultiVector b_A(
A.OperatorRangeMap(),NumVectors_AB);
396 A.Multiply(
false,x_exactA,b_A);
398 Epetra_MultiVector x_B(
B.OperatorDomainMap(),NumVectors_AB);
399 Epetra_MultiVector x_exactB(
B.OperatorDomainMap(),NumVectors_AB);
400 Epetra_MultiVector b_B(
B.OperatorRangeMap(),NumVectors_AB);
402 B.Multiply(
false,x_exactB,b_B);
404 Epetra_MultiVector x_C(
C.OperatorDomainMap(),NumVectors_C);
405 Epetra_MultiVector x_exactC(
C.OperatorDomainMap(),NumVectors_C);
406 Epetra_MultiVector b_C(
C.OperatorRangeMap(),NumVectors_C);
408 C.Multiply(
false,x_exactC,b_C);
410 std::vector<std::string> SolverType;
411 SolverType.push_back(
"Amesos_Klu");
412 SolverType.push_back(
"Amesos_Lapack");
413 SolverType.push_back(
"Amesos_Umfpack");
414 SolverType.push_back(
"Amesos_Superlu");
415 SolverType.push_back(
"Amesos_Superludist");
416 SolverType.push_back(
"Amesos_Mumps");
420 bool TestPassed =
true;
422 for (
unsigned int i = 0 ; i < SolverType.size() ; ++i) {
423 std::string
Solver = SolverType[i];
426 A, x_A, b_A, x_exactA,
427 B, x_B, b_B, x_exactB,
428 C, x_C, b_C, x_exactC);
429 TestPassed = TestPassed && ok;
433 if (
Comm.MyPID() == 0)
434 std::cout <<
"Solver " <<
Solver <<
" not available" << std::endl;
445 if (
Comm.MyPID() == 0)
446 std::cout << std::endl <<
"TEST PASSED" << std::endl << std::endl;
447 return(EXIT_SUCCESS);
450 if (
Comm.MyPID() == 0)
451 std::cout << std::endl <<
"TEST FAILED" << std::endl << std::endl;
458 int main(
int argc,
char *argv[]) {
461 MPI_Init(&argc, &argv);
462 Epetra_MpiComm
Comm(MPI_COMM_WORLD);
464 Epetra_SerialComm
Comm;
468 if (
Comm.MyPID() == 0 ) {
469 std::cout <<
"Enter a char to continue" ;
int SubMain(Epetra_Comm &Comm)
T & get(ParameterList &l, const std::string &name)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Amesos_TestRowMatrix: a class to test Epetra_RowMatrix based codes.
#define AMESOS_CHK_ERR(a)
bool TestTiming(const Amesos_BaseSolver *Solver, const Epetra_Comm &Comm)
Factory for binding a third party direct solver to an Epetra_LinearProblem.
int CreateCrsMatrix(const char *in_filename, const Epetra_Comm &Comm, Epetra_Map *&readMap, const bool transpose, const bool distribute, bool &symmetric, Epetra_CrsMatrix *&Matrix)
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.
bool CheckError(const std::string SolverType, const std::string Descriptor, const Epetra_RowMatrix &A, const Epetra_MultiVector &x, const Epetra_MultiVector &b, const Epetra_MultiVector &x_exact)
int main(int argc, char *argv[])
Amesos_BaseSolver: A pure virtual class for direct solution of real-valued double-precision operators...
bool Test(char *SolverType, Epetra_RowMatrix &A, Epetra_MultiVector &x_A, Epetra_MultiVector &b_A, Epetra_MultiVector &x_exactA, Epetra_RowMatrix &B, Epetra_MultiVector &x_B, Epetra_MultiVector &b_B, Epetra_MultiVector &x_exactB, Epetra_RowMatrix &C, Epetra_MultiVector &x_C, Epetra_MultiVector &b_C, Epetra_MultiVector &x_exactC)
bool Query(const char *ClassType)
Queries whether a given interface is available or not.