5 #include "Epetra_MpiComm.h" 7 #include "Epetra_SerialComm.h" 11 #include "Epetra_Map.h" 13 #include "Epetra_Vector.h" 14 #include "Epetra_CrsMatrix.h" 16 #include "Teuchos_ParameterList.hpp" 29 bool TestAmesos(
char ProblemType[], Teuchos::ParameterList& AmesosList,
30 bool UseTranspose, Epetra_RowMatrix* A, Epetra_MultiVector* lhs,
31 Epetra_MultiVector* rhs)
33 Epetra_LinearProblem Problem;
46 A->Multiply(UseTranspose,*lhs,*rhs);
49 Problem.SetOperator(A);
58 if (Solver->
Solve())
return(
false);
61 double d = 0.0, d_tot = 0.0;
63 for (
int i=0 ; i<lhs->Map().NumMyElements() ; ++i)
64 for (
int j = 0 ; j < lhs->NumVectors() ; ++j)
65 d += ((*lhs)[j][i] - 1.0) * ((*lhs)[j][i] - 1.0);
67 A->Comm().SumAll(&d,&d_tot,1);
70 std::vector<double> Norm(rhs->NumVectors());
72 Epetra_MultiVector Ax(*rhs);
73 A->Multiply(UseTranspose, *lhs, Ax);
74 Ax.Update(1.0, *rhs, -1.0);
77 std::string msg = ProblemType;
79 if (!A->Comm().MyPID())
81 std::cout << std::endl;
82 std::cout << msg <<
" : Using " << A->Comm().NumProc() <<
" processes, UseTranspose = " << UseTranspose << std::endl;
83 for (
int j = 0 ; j < rhs->NumVectors() ; ++j)
84 std::cout << msg <<
" : eq " << j
85 <<
", ||A x - b||_2 = " << Norm[j] << std::endl;
86 std::cout << msg <<
" : ||x_exact - x||_2 = " << sqrt(d_tot) << std::endl;
91 std::cerr << std::endl << msg <<
" WARNING : TEST FAILED!" << std::endl;
104 int main(
int argc,
char *argv[])
107 MPI_Init(&argc, &argv);
108 Epetra_MpiComm Comm(MPI_COMM_WORLD);
110 Epetra_SerialComm Comm;
113 int NumGlobalRows = 100 * Comm.NumProc();
115 Epetra_Map Map(NumGlobalRows, 0, Comm);
117 Epetra_CrsMatrix Matrix(Copy, Map, 0);
119 int NumMyRows = Map.NumMyElements();
120 int* MyGlobalElements = Map.MyGlobalElements();
126 for (
int LRID = 0 ; LRID < NumMyRows ; ++LRID)
128 int GRID = MyGlobalElements[LRID];
136 Indices[NumEntries] = GRID - 1;
137 Values[NumEntries] = -1.0;
140 if (GRID != NumGlobalRows - 1)
142 Indices[NumEntries] = GRID + 1;
143 Values[NumEntries] = -1.0;
147 Matrix.InsertGlobalValues(GRID, NumEntries, Values, Indices);
150 Matrix.FillComplete();
152 Epetra_MultiVector LHS(Map, 3);
153 Epetra_MultiVector RHS(Map, 3);
157 std::vector<std::string> SolverType;
158 SolverType.push_back(
"Amesos_Lapack");
159 SolverType.push_back(
"Amesos_Klu");
160 SolverType.push_back(
"Amesos_Umfpack");
162 SolverType.push_back(
"Amesos_Taucs");
163 SolverType.push_back(
"Amesos_Superlu");
164 SolverType.push_back(
"Amesos_Superludist");
165 SolverType.push_back(
"Amesos_Mumps");
166 SolverType.push_back(
"Amesos_Dscpack");
167 SolverType.push_back(
"Amesos_Scalapack");
173 for (
unsigned int i = 0 ; i < SolverType.size() ; ++i)
175 std::string Solver = SolverType[i];
177 if (Factory.
Query((
char*)Solver.c_str()))
181 Teuchos::ParameterList AmesosList;
185 TestAmesos((
char*)Solver.c_str(), AmesosList,
false,
186 &Matrix, &LHS, &RHS);
187 assert (res ==
true);
191 if (Solver !=
"Amesos_Superludist") {
192 Teuchos::ParameterList AmesosList;
196 TestAmesos((
char*)Solver.c_str(), AmesosList,
true,
197 &Matrix, &LHS, &RHS);
198 assert (res ==
true);
205 std::cerr << std::endl;
206 std::cerr <<
"WARNING: SOLVER `" << Solver <<
"' NOT TESTED" << std::endl;
207 std::cerr << std::endl;
int main(int argc, char *argv[])
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)
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.
Factory for binding a third party direct solver to an Epetra_LinearProblem.
virtual int SetUseTranspose(bool UseTranspose)=0
If set true, X will be set to the solution of AT X = B (not A X = B)
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.
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.