Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Test_Detailed/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_Map.h"
10#include "Epetra_Vector.h"
11#include "Epetra_Time.h"
12#include "Amesos.h"
13#include "Amesos_BaseSolver.h"
15#include "Teuchos_ParameterList.hpp"
16#include "Galeri_Maps.h"
17#include "Galeri_CrsMatrices.h"
18
19#ifdef HAVE_VALGRIND_H
20#include <valgrind/valgrind.h>
21#define HAVE_VALGRIND
22#else
23#ifdef HAVE_VALGRIND_VALGRIND_H
24#include <valgrind/valgrind.h>
25#define HAVE_VALGRIND
26#endif
27#endif
28
29using namespace Teuchos;
30using namespace Galeri;
31
32//=============================================================================
33
34bool TestTiming(const Amesos_BaseSolver* Solver,
35 const Epetra_Comm& Comm)
36{
37 // Get the timings here
38 Teuchos::ParameterList TimingsList;
39 Solver->GetTiming( TimingsList );
40
41 bool testPassed = true;
42
43 try {
44
45 // Test passes if no exception is caught
46
47 // you can find out how much time was spent in ...
48 //double sfact_time, nfact_time, solve_time;
49 //double mtx_conv_time, mtx_redist_time, vec_redist_time;
50
51 // 1) The symbolic factorization
52 // (parameter doesn't always exist)
53 //sfact_time = TimingsList.get( "Total symbolic factorization time", 0.0 );
54 TimingsList.get( "Total symbolic factorization time", 0.0 );
55
56 // 2) The numeric factorization
57 // (always exists if NumericFactorization() is called)
58 //nfact_time = Teuchos::getParameter<double>( TimingsList, "Total numeric factorization time" );
59 Teuchos::getParameter<double>( TimingsList, "Total numeric factorization time" );
60
61 // 3) Solving the linear system
62 // (always exists if Solve() is called)
63 //solve_time = Teuchos::getParameter<double>( TimingsList, "Total solve time" );
64 Teuchos::getParameter<double>( TimingsList, "Total solve time" );
65
66 // 4) Converting the matrix to the accepted format for the solver
67 // (always exists if SymbolicFactorization() is called)
68 //mtx_conv_time = Teuchos::getParameter<double>( TimingsList, "Total solve time" );
69 Teuchos::getParameter<double>( TimingsList, "Total solve time" );
70
71 // 5) Redistributing the matrix for each solve to the accepted format for the solver
72 //mtx_redist_time = TimingsList.get( "Total matrix redistribution time", 0.0 );
73 TimingsList.get( "Total matrix redistribution time", 0.0 );
74
75 // 6) Redistributing the vector for each solve to the accepted format for the solver
76 //vec_redist_time = TimingsList.get( "Total vector redistribution time", 0.0 );
77 TimingsList.get( "Total vector redistribution time", 0.0 );
78
79 }
80 catch( std::exception& e ) {
81 if (Comm.MyPID() == 0)
82 std::cout << std::endl << "Exception caught in TestTiming() : " << e.what() << std::endl;
83 testPassed = false;
84 }
85
86 return testPassed;
87
88}
89
90//=============================================================================
91bool CheckError(const std::string SolverType,
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)
97{
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);
107 Ax.Norm2(&Norm[0]);
108 bool TestPassed = false;
109 double TotalNorm = 0.0;
110 for (int i = 0 ; i < NumVectors ; ++i) {
111 TotalNorm += Norm[i];
112 }
113 if (A.Comm().MyPID() == 0)
114 std::cout << "||Ax - b|| = " << TotalNorm << std::endl;
115 if (TotalNorm < 1e-5 )
116 TestPassed = true;
117 else
118 TestPassed = false;
119
120 Ax.Update (1.0,x,-1.0,x_exact,0.0);
121 Ax.Norm2(&Norm[0]);
122 for (int i = 0 ; i < NumVectors ; ++i) {
123 TotalNorm += Norm[i];
124 }
125 if (A.Comm().MyPID() == 0)
126 std::cout << "||x - x_exact|| = " << TotalNorm << std::endl;
127 if (TotalNorm < 1e-5 )
128 TestPassed = true;
129 else
130 TestPassed = false;
131
132 return(TestPassed);
133}
134
135//=============================================================================
136bool Test(char* SolverType,
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)
143{
144 bool TestPassed = true;
145 Epetra_LinearProblem ProblemA;
146 Epetra_LinearProblem ProblemB;
147 Epetra_LinearProblem ProblemC;
148 Amesos Factory;
149 Amesos_BaseSolver* Solver;
150
151 // Test simple usage:
152 // - set problem (empty)
153 // - set matrix and vector
154 // - call Solve() directly
155 if (true) {
156 x_A.PutScalar(0.0);
157 ProblemA.SetOperator((Epetra_RowMatrix*)0);
158 ProblemA.SetLHS((Epetra_MultiVector*)0);
159 ProblemA.SetRHS((Epetra_MultiVector*)0);
160
161 Solver = Factory.Create(SolverType,ProblemA);
162
163 ProblemA.SetOperator(&A);
164 ProblemA.SetLHS(&x_A);
165 ProblemA.SetRHS(&b_A);
166
167 std::string ST = SolverType ;
168 if (! ( ST == "Amesos_Superludist" ) ) { // Kludge see bug #1141 and bug #1138
169 AMESOS_CHK_ERR(Solver->Solve());
170
171 TestPassed = TestPassed &&
172 CheckError(SolverType, "Solve() only", A,x_A,b_A,x_exactA) &&
173 TestTiming(Solver, A.Comm());
174 }
175 delete Solver;
176 }
177
178 // Test almost simple usage:
179 // - set problem (empty)
180 // - set matrix and vector
181 // - call NumericFactorization()
182 // - call Solve() directly
183
184 if (true)
185 {
186 x_A.PutScalar(0.0);
187 ProblemA.SetOperator((Epetra_RowMatrix*)0);
188 ProblemA.SetLHS((Epetra_MultiVector*)0);
189 ProblemA.SetRHS((Epetra_MultiVector*)0);
190
191 Solver = Factory.Create(SolverType,ProblemA);
192
193 ProblemA.SetOperator(&A);
195 ProblemA.SetLHS(&x_A);
196 ProblemA.SetRHS(&b_A);
197 AMESOS_CHK_ERR(Solver->Solve());
198
199 TestPassed = TestPassed &&
200 CheckError(SolverType, "NumFact() + Solve()", A,x_A,b_A,x_exactA) &&
201 TestTiming(Solver, A.Comm());
202
203 delete Solver;
204 }
205
206 // Test normal usage:
207 // - set problem (empty)
208 // - set matrix
209 // - call SymbolicFactorization() several times
210 // - call NumericFactorization() several times
211 // - set vectors
212 // - call Solve() several times
213 // Repeated calls should *not* break the object. Data are supposed
214 // to be re-computed as needed.
215
216 if (true)
217 {
218 x_A.PutScalar(0.0);
219 ProblemA.SetOperator((Epetra_RowMatrix*)0);
220 ProblemA.SetLHS((Epetra_MultiVector*)0);
221 ProblemA.SetRHS((Epetra_MultiVector*)0);
222
223 Solver = Factory.Create(SolverType,ProblemA);
224
225 ProblemA.SetOperator(&A);
233 ProblemA.SetLHS(&x_A);
234 ProblemA.SetRHS(&b_A);
235
236 AMESOS_CHK_ERR(Solver->Solve());
237 AMESOS_CHK_ERR(Solver->Solve());
238 AMESOS_CHK_ERR(Solver->Solve());
239 AMESOS_CHK_ERR(Solver->Solve());
240 AMESOS_CHK_ERR(Solver->Solve());
241 AMESOS_CHK_ERR(Solver->Solve());
242
243 TestPassed = TestPassed &&
244 CheckError(SolverType, "SymFact() + NumFact() + Solve()",
245 A,x_A,b_A,x_exactA) &&
246 TestTiming(Solver,A.Comm());
247
248 delete Solver;
249 }
250
251 // Test normal usage:
252 // - set problem (empty)
253 // - set matrix (as A)
254 // - call SymbolicFactorization()
255 // - set pointer to B
256 // - call NumericFactorization()
257 // - set vectors for B
258 // - call Solve()
259
260 if (true)
261 {
262 x_A.PutScalar(0.0);
263 ProblemA.SetOperator((Epetra_RowMatrix*)0);
264 ProblemA.SetLHS((Epetra_MultiVector*)0);
265 ProblemA.SetRHS((Epetra_MultiVector*)0);
266
267 Solver = Factory.Create(SolverType,ProblemA);
268
269 ProblemA.SetOperator(&A);
271 ProblemA.SetOperator(&B);
273 ProblemA.SetLHS(&x_B);
274 ProblemA.SetRHS(&b_B);
275
276 AMESOS_CHK_ERR(Solver->Solve());
277
278 TestPassed = TestPassed &&
279 CheckError(SolverType, "Set A, solve B", B,x_B,b_B,x_exactB) &&
280 TestTiming(Solver, A.Comm());
281
282 delete Solver;
283 }
284
285 // Construct Solver with filled ProblemA.
286 // Then, stick pointers for different vectors.
287 // a different map as well.
288
289 if (true)
290 {
291 x_C.PutScalar(0.0);
292 ProblemA.SetOperator(&A);
293 ProblemA.SetLHS(&x_A);
294 ProblemA.SetRHS(&b_A);
295
296 Solver = Factory.Create(SolverType,ProblemA);
297
298 ProblemA.SetOperator(&C);
299
300 std::string ST = SolverType ;
301 if (! ( ST == "Amesos_Superludist" ) ) { // Kludge see bug #1141 and bug #1138
302
305
306 ProblemA.SetLHS(&x_C);
307 ProblemA.SetRHS(&b_C);
308 AMESOS_CHK_ERR(Solver->Solve());
309
310 TestPassed = TestPassed &&
311 CheckError(SolverType, "Set A, Solve C", C,x_C,b_C,x_exactC) &&
312 TestTiming(Solver, A.Comm());
313 }
314 delete Solver;
315 }
316
317 // Construct Solver with filled ProblemA, call Solve().
318 // Then, replace the pointers for matrix and vectors,
319 // and call Solve() again.
320
321 if (true)
322 {
323 x_C.PutScalar(0.0);
324 ProblemA.SetOperator(&A);
325 ProblemA.SetLHS(&x_A);
326 ProblemA.SetRHS(&b_A);
327
328 Solver = Factory.Create(SolverType,ProblemA);
329
330 std::string ST = SolverType ;
331 if (! ( ST == "Amesos_Superludist" ) ) { // bug #1141 and bug #1138
332 AMESOS_CHK_ERR(Solver->Solve());
333
334 ProblemA.SetOperator(&C);
337
338 ProblemA.SetLHS(&x_C);
339 ProblemA.SetRHS(&b_C);
340 AMESOS_CHK_ERR(Solver->Solve());
341
342 TestPassed = TestPassed &&
343 CheckError(SolverType, "Solve A + Solve C", C,x_C,b_C,x_exactC) &&
344 TestTiming(Solver, A.Comm());
345 }
346 delete Solver;
347 }
348
349 return(TestPassed);
350}
351
352//=============================================================================
353
354int SubMain( Epetra_Comm &Comm ) {
355 // Creation of data
356 // A and B refer to two different matrices, with the *same*
357 // structure and *different* values. Clearly, the map is the same.
358 //
359 // C refers to a completely different matrix
360
361 int Size_AB = 20; // A and B have size Size_AB * Size_AB
362 int Size_C = 30;
363 int NumVectors_AB = 7;
364 int NumVectors_C = 13;
365
366#ifdef HAVE_VALGRIND
367 if ( RUNNING_ON_VALGRIND ) {
368 Size_AB = 6; // must be square
369 Size_C = 6;
370 NumVectors_AB = 2;
371 NumVectors_C = 3;
372 }
373#endif
374
375 Teuchos::ParameterList GaleriList;
376
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);
383
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);
387
388 Amesos_TestRowMatrix A(Matrix_A);
389 Amesos_TestRowMatrix B(Matrix_B);
390 Amesos_TestRowMatrix C(Matrix_C);
391
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);
395 x_exactA.Random();
396 A.Multiply(false,x_exactA,b_A);
397
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);
401 x_exactB.Random();
402 B.Multiply(false,x_exactB,b_B);
403
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);
407 x_exactC.Random();
408 C.Multiply(false,x_exactC,b_C);
409
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");
417 // NOTE: DSCPACK does not support Epetra_RowMatrix's
418
419 Amesos Factory;
420 bool TestPassed = true;
421
422 for (unsigned int i = 0 ; i < SolverType.size() ; ++i) {
423 std::string Solver = SolverType[i];
424 if (Factory.Query((char*)Solver.c_str())) {
425 bool ok = Test((char*)Solver.c_str(),
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;
430 }
431 else
432 {
433 if (Comm.MyPID() == 0)
434 std::cout << "Solver " << Solver << " not available" << std::endl;
435 }
436 }
437
438 delete Matrix_A;
439 delete Matrix_B;
440 delete Matrix_C;
441 delete Map_A;
442 delete Map_C;
443
444 if (TestPassed) {
445 if (Comm.MyPID() == 0)
446 std::cout << std::endl << "TEST PASSED" << std::endl << std::endl;
447 return(EXIT_SUCCESS);
448 }
449 else {
450 if (Comm.MyPID() == 0)
451 std::cout << std::endl << "TEST FAILED" << std::endl << std::endl;
452 // exit without calling MPI_Finalize() to raise an error
453 exit(EXIT_FAILURE);
454 }
455
456}
457
458int main(int argc, char *argv[]) {
459
460#ifdef HAVE_MPI
461 MPI_Init(&argc, &argv);
462 Epetra_MpiComm Comm(MPI_COMM_WORLD);
463#else
464 Epetra_SerialComm Comm;
465#endif
466
467#if 0
468 if ( Comm.MyPID() == 0 ) {
469 std::cout << "Enter a char to continue" ;
470 char any;
471 std::cin >> any ;
472 }
473 Comm.Barrier();
474#endif
475
476 int retval = SubMain( Comm ) ;
477
478#ifdef HAVE_MPI
479 MPI_Finalize();
480#endif
481
482 return retval ;
483}
#define AMESOS_CHK_ERR(a)
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[])
int SubMain(Epetra_Comm &Comm)
bool TestTiming(const Amesos_BaseSolver *Solver, const Epetra_Comm &Comm)
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)
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)
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 void GetTiming(Teuchos::ParameterList &TimingParameterList) const
Extracts timing information from the current solver and places it in the parameter list....
virtual const Epetra_Comm & Comm() const =0
Returns a pointer to the Epetra_Comm communicator associated with this operator.
virtual int NumericFactorization()=0
Performs NumericFactorization on the matrix A.
virtual int SymbolicFactorization()=0
Performs SymbolicFactorization on the matrix A.
Amesos_TestRowMatrix: a class to test Epetra_RowMatrix based codes.
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.
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
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