Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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"
13#include "Teuchos_ParameterList.hpp"
14#include "Galeri_Maps.h"
15#include "Galeri_CrsMatrices.h"
16
17using namespace Galeri;
18
19bool 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
32bool 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
104void 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#ifndef NDEBUG
128 bool res;
129#endif
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#ifndef NDEBUG
142 res =
143#endif
144 TestAmesos((char*)Solver.c_str(), AmesosList, false,
145 &A, &LHS, &RHS);
146 assert (res == true);
147 if (UseTranspose) {
148 // solve transpose with matrix
149 Teuchos::ParameterList AmesosList;
150#ifndef NDEBUG
151 res =
152#endif
153 TestAmesos((char*)Solver.c_str(), AmesosList, true,
154 &A, &LHS, &RHS);
155 assert (res == true);
156 }
157 }
158 else
159 if (!quiet && !Comm.MyPID())
160 {
161 std::cerr << std::endl;
162 std::cerr << "WARNING: SOLVER `" << Solver << "' NOT TESTED" << std::endl;
163 std::cerr << std::endl;
164 }
165 }
166
167 delete Matrix;
168 delete Map;
169}
170
171// =========== //
172// main driver //
173// =========== //
174//
175int main(int argc, char *argv[])
176{
177#ifdef HAVE_MPI
178 MPI_Init(&argc, &argv);
179 Epetra_MpiComm Comm(MPI_COMM_WORLD);
180#else
181 Epetra_SerialComm Comm;
182#endif
183
184 if ( argc > 1 && argv[1][0] == '-' && argv[1][1] == 'q' ) quiet = true ;
185
186 // NOTE: DSCPACK does not support RowMatrix's.
187
188 if (true)
189 {
190 // non-symmetric matrix, test A and A^T
191 std::vector<std::string> SolverType;
192 SolverType.push_back("Amesos_Lapack");
193 SolverType.push_back("Amesos_Klu");
194 SolverType.push_back("Amesos_Umfpack");
195 SolverType.push_back("Amesos_Superlu");
196// SolverType.push_back("Amesos_Mumps"); Bug #1896
197 SolverType.push_back("Amesos_Scalapack");
198 driver(Comm, false, true, SolverType);
199 }
200
201 if (true)
202 {
203 // non-symmetric matrix, test only A
204 std::vector<std::string> SolverType;
205 //
206 // SolverType.push_back("Amesos_Pardiso"); bug #1994
207 SolverType.push_back("Amesos_Superludist");
208 driver(Comm, false, false, SolverType);
209 }
210
211 // I have some problems with Taucs with LAM
212 if (true)
213 {
214 // symmetric
215 std::vector<std::string> SolverType;
216 SolverType.push_back("Amesos_Taucs");
217 driver(Comm, true, false, SolverType);
218 }
219
220#ifdef HAVE_MPI
221 MPI_Finalize();
222#endif
223 return(0);
224}
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
int main(int argc, char *argv[])
void driver(Epetra_Comm &Comm, const bool IsSymmetric, const bool UseTranspose, std::vector< std::string > &SolverType)
bool TestAmesos(char ProblemType[], Teuchos::ParameterList &AmesosList, bool UseTranspose, Epetra_RowMatrix *A, Epetra_MultiVector *lhs, Epetra_MultiVector *rhs)
Amesos_BaseSolver: A pure virtual class for direct solution of real-valued double-precision operators...
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 SetUseTranspose(bool UseTranspose)=0
If set true, X will be set to the solution of AT X = B (not A X = B)
virtual int SetParameters(Teuchos::ParameterList &ParameterList)=0
Updates internal variables.
Amesos_TestRowMatrix: a class to test Epetra_RowMatrix based codes.
Factory for binding a third party direct solver to an Epetra_LinearProblem.
Definition: Amesos.h:44
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.
Definition: Amesos.cpp:69
bool Query(const char *ClassType)
Queries whether a given interface is available or not.
Definition: Amesos.cpp:193