Amesos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
14 #include "Amesos_TestRowMatrix.h"
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 
29 using namespace Teuchos;
30 using namespace Galeri;
31 
32 //=============================================================================
33 
34 bool 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 //=============================================================================
91 bool 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 //=============================================================================
136 bool Test(char* SolverType,
138  Epetra_MultiVector& b_A, Epetra_MultiVector& x_exactA,
140  Epetra_MultiVector& b_B, Epetra_MultiVector& x_exactB,
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;
150  std::string ST = SolverType ;
151 
152  // Test simple usage:
153  // - set problem (empty)
154  // - set matrix and vector
155  // - call Solve() directly
156  if (true) {
157  x_A.PutScalar(0.0);
158  ProblemA.SetOperator((Epetra_RowMatrix*)0);
159  ProblemA.SetLHS((Epetra_MultiVector*)0);
160  ProblemA.SetRHS((Epetra_MultiVector*)0);
161 
162  Solver = Factory.Create(SolverType,ProblemA);
163  if ( ST == "Amesos_CssMKL" ) {
165  List.set("Reindex", true);
166  Solver->SetParameters(List);
167  }
168 
169  ProblemA.SetOperator(&A);
170  ProblemA.SetLHS(&x_A);
171  ProblemA.SetRHS(&b_A);
172 
173  if (! ( ST == "Amesos_Superludist" ) ) { // Kludge see bug #1141 and bug #1138
174  AMESOS_CHK_ERR(Solver->Solve());
175 
176  TestPassed = TestPassed &&
177  CheckError(SolverType, "Solve() only", A,x_A,b_A,x_exactA) &&
178  TestTiming(Solver, A.Comm());
179  }
180  delete Solver;
181  }
182 
183  // Test almost simple usage:
184  // - set problem (empty)
185  // - set matrix and vector
186  // - call NumericFactorization()
187  // - call Solve() directly
188 
189  if (true)
190  {
191  x_A.PutScalar(0.0);
192  ProblemA.SetOperator((Epetra_RowMatrix*)0);
193  ProblemA.SetLHS((Epetra_MultiVector*)0);
194  ProblemA.SetRHS((Epetra_MultiVector*)0);
195 
196  Solver = Factory.Create(SolverType,ProblemA);
197  if ( ST == "Amesos_CssMKL" ) {
199  List.set("Reindex", true);
200  Solver->SetParameters(List);
201  }
202 
203  ProblemA.SetOperator(&A);
205  ProblemA.SetLHS(&x_A);
206  ProblemA.SetRHS(&b_A);
207  AMESOS_CHK_ERR(Solver->Solve());
208 
209  TestPassed = TestPassed &&
210  CheckError(SolverType, "NumFact() + Solve()", A,x_A,b_A,x_exactA) &&
211  TestTiming(Solver, A.Comm());
212 
213  delete Solver;
214  }
215 
216  // Test normal usage:
217  // - set problem (empty)
218  // - set matrix
219  // - call SymbolicFactorization() several times
220  // - call NumericFactorization() several times
221  // - set vectors
222  // - call Solve() several times
223  // Repeated calls should *not* break the object. Data are supposed
224  // to be re-computed as needed.
225 
226  if (true)
227  {
228  x_A.PutScalar(0.0);
229  ProblemA.SetOperator((Epetra_RowMatrix*)0);
230  ProblemA.SetLHS((Epetra_MultiVector*)0);
231  ProblemA.SetRHS((Epetra_MultiVector*)0);
232 
233  Solver = Factory.Create(SolverType,ProblemA);
234  if ( ST == "Amesos_CssMKL" ) {
235  // need to rindex since "interlaced"
237  List.set("Reindex", true);
238  Solver->SetParameters(List);
239  }
240 
241  ProblemA.SetOperator(&A);
249  ProblemA.SetLHS(&x_A);
250  ProblemA.SetRHS(&b_A);
251 
252  AMESOS_CHK_ERR(Solver->Solve());
253  AMESOS_CHK_ERR(Solver->Solve());
254  AMESOS_CHK_ERR(Solver->Solve());
255  AMESOS_CHK_ERR(Solver->Solve());
256  AMESOS_CHK_ERR(Solver->Solve());
257  AMESOS_CHK_ERR(Solver->Solve());
258 
259  TestPassed = TestPassed &&
260  CheckError(SolverType, "SymFact() + NumFact() + Solve()",
261  A,x_A,b_A,x_exactA) &&
262  TestTiming(Solver,A.Comm());
263 
264  delete Solver;
265  }
266 
267  // Test normal usage:
268  // - set problem (empty)
269  // - set matrix (as A)
270  // - call SymbolicFactorization()
271  // - set pointer to B
272  // - call NumericFactorization()
273  // - set vectors for B
274  // - call Solve()
275 
276  if (true)
277  {
278  x_A.PutScalar(0.0);
279  ProblemA.SetOperator((Epetra_RowMatrix*)0);
280  ProblemA.SetLHS((Epetra_MultiVector*)0);
281  ProblemA.SetRHS((Epetra_MultiVector*)0);
282 
283  Solver = Factory.Create(SolverType,ProblemA);
284  if ( ST == "Amesos_CssMKL" ) {
285  // need to rindex since "interlaced"
287  List.set("Reindex", true);
288  Solver->SetParameters(List);
289  }
290 
291  ProblemA.SetOperator(&A);
293  ProblemA.SetOperator(&B);
295  ProblemA.SetLHS(&x_B);
296  ProblemA.SetRHS(&b_B);
297 
298  AMESOS_CHK_ERR(Solver->Solve());
299 
300  TestPassed = TestPassed &&
301  CheckError(SolverType, "Set A, solve B", B,x_B,b_B,x_exactB) &&
302  TestTiming(Solver, A.Comm());
303 
304  delete Solver;
305  }
306 
307  // Construct Solver with filled ProblemA.
308  // Then, stick pointers for different vectors.
309  // a different map as well.
310 
311  if (true)
312  {
313  x_C.PutScalar(0.0);
314  ProblemA.SetOperator(&A);
315  ProblemA.SetLHS(&x_A);
316  ProblemA.SetRHS(&b_A);
317 
318  Solver = Factory.Create(SolverType,ProblemA);
319  if ( ST == "Amesos_CssMKL" ) {
320  // need to rindex since "interlaced"
322  List.set("Reindex", true);
323  Solver->SetParameters(List);
324  }
325 
326  ProblemA.SetOperator(&C);
327 
328  std::string ST = SolverType ;
329  if (! ( ST == "Amesos_Superludist" ) ) { // Kludge see bug #1141 and bug #1138
330 
333 
334  ProblemA.SetLHS(&x_C);
335  ProblemA.SetRHS(&b_C);
336  AMESOS_CHK_ERR(Solver->Solve());
337 
338  TestPassed = TestPassed &&
339  CheckError(SolverType, "Set A, Solve C", C,x_C,b_C,x_exactC) &&
340  TestTiming(Solver, A.Comm());
341  }
342  delete Solver;
343  }
344 
345  // Construct Solver with filled ProblemA, call Solve().
346  // Then, replace the pointers for matrix and vectors,
347  // and call Solve() again.
348 
349  if (true)
350  {
351  x_C.PutScalar(0.0);
352  ProblemA.SetOperator(&A);
353  ProblemA.SetLHS(&x_A);
354  ProblemA.SetRHS(&b_A);
355 
356  Solver = Factory.Create(SolverType,ProblemA);
357  if ( ST == "Amesos_CssMKL" ) {
358  // need to rindex since "interlaced"
360  List.set("Reindex", true);
361  Solver->SetParameters(List);
362  }
363 
364  std::string ST = SolverType ;
365  if (! ( ST == "Amesos_Superludist" ) ) { // bug #1141 and bug #1138
366  AMESOS_CHK_ERR(Solver->Solve());
367 
368  ProblemA.SetOperator(&C);
371 
372  ProblemA.SetLHS(&x_C);
373  ProblemA.SetRHS(&b_C);
374  AMESOS_CHK_ERR(Solver->Solve());
375 
376  TestPassed = TestPassed &&
377  CheckError(SolverType, "Solve A + Solve C", C,x_C,b_C,x_exactC) &&
378  TestTiming(Solver, A.Comm());
379  }
380  delete Solver;
381  }
382 
383  return(TestPassed);
384 }
385 
386 //=============================================================================
387 
389  // Creation of data
390  // A and B refer to two different matrices, with the *same*
391  // structure and *different* values. Clearly, the map is the same.
392  //
393  // C refers to a completely different matrix
394 
395  int Size_AB = 20; // A and B have size Size_AB * Size_AB
396  int Size_C = 30;
397  int NumVectors_AB = 7;
398  int NumVectors_C = 13;
399 
400 #ifdef HAVE_VALGRIND
401  if ( RUNNING_ON_VALGRIND ) {
402  Size_AB = 6; // must be square
403  Size_C = 6;
404  NumVectors_AB = 2;
405  NumVectors_C = 3;
406  }
407 #endif
408 
409  Teuchos::ParameterList GaleriList;
410 
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);
415  Epetra_CrsMatrix* Matrix_A = CreateCrsMatrix("Recirc2D", Map_A, GaleriList);
416  Epetra_CrsMatrix* Matrix_B = CreateCrsMatrix("Laplace2D", Map_A, GaleriList);
417 
418  GaleriList.set("n", Size_C);
419  Epetra_Map* Map_C = CreateMap("Interlaced", Comm, GaleriList);
420  Epetra_CrsMatrix* Matrix_C = CreateCrsMatrix("Minij", Map_C, GaleriList);
421 
422  Amesos_TestRowMatrix A(Matrix_A);
423  Amesos_TestRowMatrix B(Matrix_B);
424  Amesos_TestRowMatrix C(Matrix_C);
425 
426  Epetra_MultiVector x_A(A.OperatorDomainMap(),NumVectors_AB);
427  Epetra_MultiVector x_exactA(A.OperatorDomainMap(),NumVectors_AB);
428  Epetra_MultiVector b_A(A.OperatorRangeMap(),NumVectors_AB);
429  x_exactA.Random();
430  A.Multiply(false,x_exactA,b_A);
431 
432  Epetra_MultiVector x_B(B.OperatorDomainMap(),NumVectors_AB);
433  Epetra_MultiVector x_exactB(B.OperatorDomainMap(),NumVectors_AB);
434  Epetra_MultiVector b_B(B.OperatorRangeMap(),NumVectors_AB);
435  x_exactB.Random();
436  B.Multiply(false,x_exactB,b_B);
437 
438  Epetra_MultiVector x_C(C.OperatorDomainMap(),NumVectors_C);
439  Epetra_MultiVector x_exactC(C.OperatorDomainMap(),NumVectors_C);
440  Epetra_MultiVector b_C(C.OperatorRangeMap(),NumVectors_C);
441  x_exactC.Random();
442  C.Multiply(false,x_exactC,b_C);
443 
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");
451  // NOTE: DSCPACK does not support Epetra_RowMatrix's
452 
453  Amesos Factory;
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;
463  }
464  else
465  {
466  if (Comm.MyPID() == 0)
467  std::cout << "Solver " << Solver << " not available" << std::endl;
468  }
469  }
470 
471  delete Matrix_A;
472  delete Matrix_B;
473  delete Matrix_C;
474  delete Map_A;
475  delete Map_C;
476 
477  // Testing CSS MKL
478  // Same as above, but using Laplace2D instead of dense Minij for C
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);
484  Epetra_CrsMatrix* Matrix_A_css = CreateCrsMatrix("Recirc2D", Map_A_css, GaleriList);
485  Epetra_CrsMatrix* Matrix_B_css = CreateCrsMatrix("Laplace2D", Map_A_css, GaleriList);
486 
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);
491  Epetra_CrsMatrix* Matrix_C_css = CreateCrsMatrix("Laplace2D", Map_C_css, GaleriList);
492 
493  Amesos_TestRowMatrix A_css(Matrix_A_css);
494  Amesos_TestRowMatrix B_css(Matrix_B_css);
495  Amesos_TestRowMatrix C_css(Matrix_C_css);
496 
497  Epetra_MultiVector x_A_css(A_css.OperatorDomainMap(),NumVectors_AB);
498  Epetra_MultiVector x_exactA_css(A_css.OperatorDomainMap(),NumVectors_AB);
499  Epetra_MultiVector b_A_css(A_css.OperatorRangeMap(),NumVectors_AB);
500  x_exactA_css.Random();
501  A_css.Multiply(false,x_exactA_css,b_A_css);
502 
503  Epetra_MultiVector x_B_css(B_css.OperatorDomainMap(),NumVectors_AB);
504  Epetra_MultiVector x_exactB_css(B_css.OperatorDomainMap(),NumVectors_AB);
505  Epetra_MultiVector b_B_css(B_css.OperatorRangeMap(),NumVectors_AB);
506  x_exactB_css.Random();
507  B_css.Multiply(false,x_exactB_css,b_B_css);
508 
509  Epetra_MultiVector x_C_css(C_css.OperatorDomainMap(),NumVectors_C);
510  Epetra_MultiVector x_exactC_css(C_css.OperatorDomainMap(),NumVectors_C);
511  Epetra_MultiVector b_C_css(C_css.OperatorRangeMap(),NumVectors_C);
512  x_exactC_css.Random();
513  C_css.Multiply(false,x_exactC_css,b_C_css);
514 
515  // Passing CrsMatrix for reindexing
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;
521  }
522  if (TestPassed) {
523  if (Comm.MyPID() == 0)
524  std::cout << std::endl << "TEST PASSED" << std::endl << std::endl;
525  return(EXIT_SUCCESS);
526  }
527  else {
528  if (Comm.MyPID() == 0)
529  std::cout << std::endl << "TEST FAILED" << std::endl << std::endl;
530  // exit without calling MPI_Finalize() to raise an error
531  exit(EXIT_FAILURE);
532  }
533 
534 }
535 
536 int main(int argc, char *argv[]) {
537 
538 #ifdef HAVE_MPI
539  MPI_Init(&argc, &argv);
540  Epetra_MpiComm Comm(MPI_COMM_WORLD);
541 #else
542  Epetra_SerialComm Comm;
543 #endif
544 
545 #if 0
546  if ( Comm.MyPID() == 0 ) {
547  std::cout << "Enter a char to continue" ;
548  char any;
549  std::cin >> any ;
550  }
551  Comm.Barrier();
552 #endif
553 
554  int retval = SubMain( Comm ) ;
555 
556 #ifdef HAVE_MPI
557  MPI_Finalize();
558 #endif
559 
560  return retval ;
561 }
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)
int MyPID() const
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.
Definition: Amesos.h:50
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
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.
Definition: Amesos.cpp:72
void Barrier() const
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.
Definition: Amesos.cpp:205