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>
29 using namespace Teuchos;
30 using namespace Galeri;
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,
99 std::cout << std::endl <<
"*** " << SolverType <<
" : "
100 << Descriptor <<
" ***" << std::endl;
101 std::vector<double> Norm;
102 int NumVectors = x.NumVectors();
103 Norm.resize(NumVectors);
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];
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];
126 std::cout <<
"||x - x_exact|| = " << TotalNorm << std::endl;
127 if (TotalNorm < 1e-5 )
144 bool TestPassed =
true;
150 std::string ST = SolverType ;
162 Solver = Factory.
Create(SolverType,ProblemA);
163 if ( ST ==
"Amesos_CssMKL" ) {
165 List.
set(
"Reindex",
true);
173 if (! ( ST ==
"Amesos_Superludist" ) ) {
176 TestPassed = TestPassed &&
177 CheckError(SolverType,
"Solve() only", A,x_A,b_A,x_exactA) &&
196 Solver = Factory.
Create(SolverType,ProblemA);
197 if ( ST ==
"Amesos_CssMKL" ) {
199 List.
set(
"Reindex",
true);
209 TestPassed = TestPassed &&
210 CheckError(SolverType,
"NumFact() + Solve()", A,x_A,b_A,x_exactA) &&
233 Solver = Factory.
Create(SolverType,ProblemA);
234 if ( ST ==
"Amesos_CssMKL" ) {
237 List.
set(
"Reindex",
true);
259 TestPassed = TestPassed &&
260 CheckError(SolverType,
"SymFact() + NumFact() + Solve()",
261 A,x_A,b_A,x_exactA) &&
283 Solver = Factory.
Create(SolverType,ProblemA);
284 if ( ST ==
"Amesos_CssMKL" ) {
287 List.
set(
"Reindex",
true);
300 TestPassed = TestPassed &&
301 CheckError(SolverType,
"Set A, solve B", B,x_B,b_B,x_exactB) &&
318 Solver = Factory.
Create(SolverType,ProblemA);
319 if ( ST ==
"Amesos_CssMKL" ) {
322 List.
set(
"Reindex",
true);
328 std::string ST = SolverType ;
329 if (! ( ST ==
"Amesos_Superludist" ) ) {
338 TestPassed = TestPassed &&
339 CheckError(SolverType,
"Set A, Solve C", C,x_C,b_C,x_exactC) &&
356 Solver = Factory.
Create(SolverType,ProblemA);
357 if ( ST ==
"Amesos_CssMKL" ) {
360 List.
set(
"Reindex",
true);
364 std::string ST = SolverType ;
365 if (! ( ST ==
"Amesos_Superludist" ) ) {
376 TestPassed = TestPassed &&
377 CheckError(SolverType,
"Solve A + Solve C", C,x_C,b_C,x_exactC) &&
397 int NumVectors_AB = 7;
398 int NumVectors_C = 13;
401 if ( RUNNING_ON_VALGRIND ) {
411 GaleriList.
set(
"n", Size_AB * Size_AB);
412 GaleriList.
set(
"nx", Size_AB);
413 GaleriList.
set(
"ny", Size_AB);
414 Epetra_Map* Map_A = CreateMap(
"Interlaced", Comm, GaleriList);
418 GaleriList.
set(
"n", Size_C);
419 Epetra_Map* Map_C = CreateMap(
"Interlaced", Comm, GaleriList);
444 std::vector<std::string> SolverType;
445 SolverType.push_back(
"Amesos_Klu");
446 SolverType.push_back(
"Amesos_Lapack");
447 SolverType.push_back(
"Amesos_Umfpack");
448 SolverType.push_back(
"Amesos_Superlu");
449 SolverType.push_back(
"Amesos_Superludist");
450 SolverType.push_back(
"Amesos_Mumps");
454 bool TestPassed =
true;
455 for (
unsigned int i = 0 ; i < SolverType.size() ; ++i) {
456 std::string
Solver = SolverType[i];
457 if (Factory.
Query((
char*)Solver.c_str())) {
458 bool ok =
Test((
char*)Solver.c_str(),
459 A, x_A, b_A, x_exactA,
460 B, x_B, b_B, x_exactB,
461 C, x_C, b_C, x_exactC);
462 TestPassed = TestPassed && ok;
466 if (Comm.
MyPID() == 0)
467 std::cout <<
"Solver " << Solver <<
" not available" << std::endl;
479 if (Factory.
Query(
"Amesos_CssMKL")) {
480 GaleriList.
set(
"n", Size_AB * Size_AB);
481 GaleriList.
set(
"nx", Size_AB);
482 GaleriList.
set(
"ny", Size_AB);
483 Epetra_Map* Map_A_css = CreateMap(
"Interlaced", Comm, GaleriList);
487 GaleriList.
set(
"n", Size_C * Size_C);
488 GaleriList.
set(
"nx", Size_C);
489 GaleriList.
set(
"ny", Size_C);
490 Epetra_Map* Map_C_css = CreateMap(
"Interlaced", Comm, GaleriList);
500 x_exactA_css.Random();
501 A_css.
Multiply(
false,x_exactA_css,b_A_css);
506 x_exactB_css.Random();
507 B_css.
Multiply(
false,x_exactB_css,b_B_css);
512 x_exactC_css.Random();
513 C_css.
Multiply(
false,x_exactC_css,b_C_css);
516 bool ok =
Test(
"Amesos_CssMKL",
517 *Matrix_A_css, x_A_css, b_A_css, x_exactA_css,
518 *Matrix_B_css, x_B_css, b_B_css, x_exactB_css,
519 *Matrix_C_css, x_C_css, b_C_css, x_exactC_css);
520 TestPassed = TestPassed && ok;
523 if (Comm.
MyPID() == 0)
524 std::cout << std::endl <<
"TEST PASSED" << std::endl << std::endl;
525 return(EXIT_SUCCESS);
528 if (Comm.
MyPID() == 0)
529 std::cout << std::endl <<
"TEST FAILED" << std::endl << std::endl;
536 int main(
int argc,
char *argv[]) {
539 MPI_Init(&argc, &argv);
546 if ( Comm.
MyPID() == 0 ) {
547 std::cout <<
"Enter a char to continue" ;
void SetLHS(Epetra_MultiVector *X)
void SetOperator(Epetra_RowMatrix *A)
int SubMain(Epetra_Comm &Comm)
virtual int Solve()=0
Solves A X = B (or AT x = B)
T & get(ParameterList &l, const std::string &name)
virtual int NumericFactorization()=0
Performs NumericFactorization on the matrix A.
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual int SymbolicFactorization()=0
Performs SymbolicFactorization on the matrix A.
virtual int SetParameters(Teuchos::ParameterList &ParameterList)=0
Updates internal variables.
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
Amesos_TestRowMatrix: a class to test Epetra_RowMatrix based codes.
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
virtual int MyPID() const =0
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_RowMatrix multiplied by a Epetra_MultiVector X in Y.
virtual const Epetra_Comm & Comm() const =0
#define AMESOS_CHK_ERR(a)
bool TestTiming(const Amesos_BaseSolver *Solver, const Epetra_Comm &Comm)
int main(int argc, char *argv[])
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)
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
virtual void GetTiming(Teuchos::ParameterList &TimingParameterList) const
Extracts timing information from the current solver and places it in the parameter list...
void SetRHS(Epetra_MultiVector *B)
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)
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.