Amesos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TestOptions.cpp
Go to the documentation of this file.
1 //
2 // To run this under valgrind, try:
3 // valgrind --suppressions=../Test_Basic/Suppressions --gen-suppressions=yes --leak-check=yes --show-reachable=yes ./TestOptions.exe -v
4 //
5 // To run this with valgrind under mpirun,
6 // mpirun -np 2 valgrind --log-file=TestOpt.logfile --suppressions=../Test_Basic/Suppressions --gen-suppressions=yes --leak-check=yes --show-reachable=yes ./TestOptions.exe -v
7 //
8 // test/scripts/daily/serial/TestMemoryLeaks[.exe] performs an automated test for memory leaks
9 // using valgrind and this code. To run TestMemoryLeaks, cd to test/TestOptions in the
10 // build directory and type ../scripts/daily/serial/TestMemoryLeaks.exe. The output is stored
11 // in amesos/logLinux.txt.
12 //
13 //
14 
15 // TestOptions tests all options for each Amesos Class on a limited number
16 // of matrices.
17 //
18 
19 // TestOneMatrix - Test one matrix
20 // - Distributed vs not distributed -
21 // - Transpose vs not transposed -
22 // TestAllClasses - Test one matrix and one setting of distributed and transpose
23 // - Calls TestOtherClasses (one for each Amesos class) and TestSuperludist
24 // TestOtherClasses
25 // TestSuperludist
26 // TestScalapack
27 //
28 //
29 // Todo:
30 // Write TestKlu, TestSuperlu, TestScalapack, TestUmfpack, TestDscpack
31 // Enable tests of various parameter options
32 // Make it test all four codes (DSCPACK, UMFPACK, SuperLU_DIST, KLU )
33 // Valgrind it
34 // Enable tests of transpose and distributed matrices - DONE
35 // Enable FACTOR_B - DONE
36 //
37 
38 #include "Trilinos_Util.h"
39 #include "Trilinos_Util_ReadMatrixMarket2Epetra.h"
40 #include "Amesos.h"
42 #include "Amesos_BaseSolver.h"
43 #include "Epetra_LinearProblem.h"
44 #include "Epetra_CrsMatrix.h"
45 #include "Epetra_Map.h"
46 #include "Epetra_Vector.h"
47 #include "Epetra_Export.h"
48 #include "Amesos_ConfigDefs.h"
49 #ifdef HAVE_AMESOS_UMFPACK
50 #include "Amesos_Umfpack.h"
51 #endif
52 #include "CrsMatrixTranspose.h"
53 #include "TestAllClasses.h"
54 #include <string>
55 #include "Teuchos_RCP.hpp"
56 #include "NewMatNewMap.h"
57 #ifdef EPETRA_MPI
58 #include "Epetra_MpiComm.h"
59 #else
60 #include "Epetra_SerialComm.h"
61 #endif
62 
63 #if 0
64 
65 #ifdef HAVE_VALGRIND_H
66 #include <valgrind/valgrind.h>
67 #define HAVE_VALGRIND
68 #else
69 #ifdef HAVE_VALGRIND_VALGRIND_H
70 #include <valgrind/valgrind.h>
71 #define HAVE_VALGRIND
72 #endif
73 #endif
74 
75 #endif
76 
77 std::vector<string> AmesosClasses;
78 
80 
81 int CreateCrsMatrix( const char *in_filename, const Epetra_Comm &Comm,
82  Epetra_Map *& readMap,
83  const bool transpose, const bool distribute,
84  bool& symmetric, Epetra_CrsMatrix *& Matrix ) {
85 
86  Epetra_CrsMatrix * readA = 0;
87  Epetra_Vector * readx = 0;
88  Epetra_Vector * readb = 0;
89  Epetra_Vector * readxexact = 0;
90 
91  //
92  // This hack allows TestOptions to be run from either the test/TestOptions/ directory or from
93  // the test/ directory (as it is in nightly testing and in make "run-tests")
94  //
95  FILE *in_file = fopen( in_filename, "r");
96 
97  const char *filename;
98  if (in_file == NULL )
99  filename = &in_filename[1] ; // Strip off ithe "." from
100  // "../" and try again
101  else {
102  filename = in_filename ;
103  fclose( in_file );
104  }
105 
106  symmetric = false ;
107  std::string FileName (filename);
108 
109  int FN_Size = FileName.size() ;
110  std::string LastFiveBytes = FileName.substr( EPETRA_MAX(0,FN_Size-5), FN_Size );
111  std::string LastFourBytes = FileName.substr( EPETRA_MAX(0,FN_Size-4), FN_Size );
112 
113  if ( LastFiveBytes == ".triU" ) {
114  // Call routine to read in unsymmetric Triplet matrix
115  EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename, false, Comm, readMap, readA, readx,
116  readb, readxexact) );
117  symmetric = false;
118  } else {
119  if ( LastFiveBytes == ".triS" ) {
120  // Call routine to read in symmetric Triplet matrix
121  EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename, true, Comm, readMap, readA, readx,
122  readb, readxexact) );
123  symmetric = true;
124  } else {
125  if ( LastFourBytes == ".mtx" ) {
126  EPETRA_CHK_ERR( Trilinos_Util_ReadMatrixMarket2Epetra( filename, Comm, readMap,
127  readA, readx, readb, readxexact) );
128  in_file = fopen( filename, "r");
129  assert (in_file != NULL) ; // Checked in Trilinos_Util_CountMatrixMarket()
130  const int BUFSIZE = 800 ;
131  char buffer[BUFSIZE] ;
132  fgets( buffer, BUFSIZE, in_file ) ; // Pick symmetry info off of this std::string
133  std::string headerline1 = buffer;
134 #ifdef TFLOP
135  if ( headerline1.find("symmetric") < BUFSIZE ) symmetric = true;
136 #else
137  if ( headerline1.find("symmetric") != std::string::npos) symmetric = true;
138 
139 #endif
140  fclose(in_file);
141 
142  } else {
143  // Call routine to read in HB problem
144  Trilinos_Util_ReadHb2Epetra( filename, Comm, readMap, readA, readx,
145  readb, readxexact) ;
146  if ( LastFourBytes == ".rsa" ) symmetric = true ;
147  }
148  }
149  }
150 
151 
152  if ( readb ) delete readb;
153  if ( readx ) delete readx;
154  if ( readxexact ) delete readxexact;
155 
156  Epetra_CrsMatrix *serialA ;
157  Epetra_CrsMatrix *transposeA;
158 
159  int ierr = 0;
160 
161  if ( transpose ) {
162  transposeA = new Epetra_CrsMatrix( Copy, *readMap, 0 );
163  ierr = CrsMatrixTranspose( readA, transposeA );
164  assert( ierr == 0 );
165  serialA = transposeA ;
166  delete readA;
167  readA = 0 ;
168  } else {
169  serialA = readA ;
170  }
171 
172  assert( (void *) &serialA->Graph() ) ;
173  assert( (void *) &serialA->RowMap() ) ;
174  assert( serialA->RowMap().SameAs(*readMap) ) ;
175 
176  if ( distribute ) {
177  // Create uniform distributed map
178  Epetra_Map* mapPtr = 0;
179 
180  if(readMap->GlobalIndicesInt())
181  mapPtr = new Epetra_Map((int) readMap->NumGlobalElements(), 0, Comm);
182  else if(readMap->GlobalIndicesLongLong())
183  mapPtr = new Epetra_Map(readMap->NumGlobalElements(), 0, Comm);
184  else
185  assert(false);
186 
187  Epetra_Map& DistMap = *mapPtr;
188 
189  // Create Exporter to distribute read-in matrix and vectors
190  Epetra_Export exporter( *readMap, DistMap );
191 
192  Epetra_CrsMatrix *Amat = new Epetra_CrsMatrix( Copy, DistMap, 0 );
193  Amat->Export(*serialA, exporter, Add);
194  ierr = Amat->FillComplete();
195  assert(ierr == 0);
196 
197  Matrix = Amat;
198  //
199  // Make sure that deleting Amat->RowMap() will delete map
200  //
201  // Bug: We can't manage to delete map his way anyway,
202  // and this fails on tranposes, so for now I just accept
203  // the memory loss.
204  // assert( &(Amat->RowMap()) == map ) ;
205  delete readMap;
206  readMap = 0 ;
207  delete serialA;
208  delete mapPtr;
209  } else {
210 
211  Matrix = serialA;
212  }
213 
214 
215  return 0;
216 }
217 
218 int TestErrors( const std::vector<bool> AmesosClassesInstalled,
219  const char *filename,
220 #ifdef EPETRA_MPI
222 #else
223  Epetra_SerialComm& Comm,
224 #endif
225  const bool verbose,
226  int &NumTests ) {
227 
228  int NumErrors =0 ;
229  double error = -13;
230  double residual = -13;
231 
232 
233  for ( int iterTrans =0 ; iterTrans < 2; iterTrans++ ) {
234  const bool transpose = iterTrans == 1 ;
235 
236  const bool distribute = 1;
237 
238  const int iterRowindex = 0;
239  const int iterColindex = 0 ;
240  const int iterRangemap = 0 ;
241  const int iterDomainmap = 0 ;
242  const int iterDiagonalOpts = 0 ;
243  bool printit = true ;
244  if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ;
245  const int EpetraMatrixType = 0 ;
246  bool symmetric = true;
247  const int iterDist = 0 ;
248 
249  Epetra_CrsMatrix *Amat = 0 ;
250  Epetra_Map *readMap = 0 ;
251  CreateCrsMatrix( filename, Comm, readMap, transpose, distribute, symmetric, Amat ) ;
252 
253  if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ <<
254  " Creating matrix from " <<
255  " filename = " << filename <<
256  " symmetric = " << symmetric <<
257  " distribute = " << distribute <<
258  " iterRowindex = " << iterRowindex <<
259  " iterColindex = " << iterColindex <<
260  " iterRangemap = " << iterRangemap <<
261  " iterDomainmap = " << iterDomainmap <<
262  " EpetraMatrixType = " << EpetraMatrixType <<
263  " iterDiagonalOpts = " << iterDiagonalOpts <<
264  " transpose = " << transpose
265  << " iterDist = " << iterDist << std::endl ;
266 
267  if ( iterDiagonalOpts ) Comm.SetTracebackMode(1); // Turn off positive Epetra warnings (i.e. iniefficient code, such as memory re-allocation)
268 
269  if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ;
270 
271  RCP<Epetra_CrsMatrix> Bmat = NewMatNewMap( *Amat,
272  iterDiagonalOpts,
273  iterRowindex,
274  iterColindex,
275  iterRangemap,
276  iterDomainmap
277  ) ;
278  if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ;
279  Comm.SetTracebackMode(2);
280 
281  if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ <<
282  " filename = " << filename <<
283  " symmetric = " << symmetric <<
284  " distribute = " << distribute <<
285  " iterRowindex = " << iterRowindex <<
286  " iterColindex = " << iterColindex <<
287  " iterRangemap = " << iterRangemap <<
288  " iterDomainmap = " << iterDomainmap <<
289  " EpetraMatrixType = " << EpetraMatrixType <<
290  " iterDiagonalOpts = " << iterDiagonalOpts <<
291  " transpose = " << transpose << " iterDist = " << iterDist << std::endl ;
292 
293  //
294  // This causes a failure in Amesos_Superludist:
295  Epetra_CrsMatrix* Cmat = &*Bmat;
296  // Epetra_CrsMatrix* Cmat = Amat ;
297 
298 
299  const int Level = 1;
300  const double MaxError = 1e-3;
301 
302  int NumTheseTests = 0 ;
303  if ( verbose ) {
304  std::cout << " About to test " << filename
305  << __FILE__ << "::" << __LINE__
306  << " EpetraMatrixType = " << EpetraMatrixType
307  << (transpose?" transpose":"" )
308  << (distribute?" distribute":"" )
309  << std::endl ;
310  }
311  int Error = TestAllClasses( AmesosClasses, EpetraMatrixType,
312  AmesosClassesInstalled,
313  Cmat,
314  transpose ,
315  verbose,
316  symmetric,
317  Level,
318  MaxError,
319  iterDiagonalOpts,
320  iterRowindex,
321  iterColindex,
322  iterRangemap,
323  iterDomainmap,
324  distribute,
325  filename,
326  error,
327  residual,
328  NumTheseTests ) ;
329  NumTests += NumTheseTests ;
330  NumErrors += Error ;
331  // BUG: Memory leak
332  // delete &(Amat->RowMap()) ;
333  if ( Amat ) delete Amat ;
334  if ( readMap ) delete readMap ;
335  }
336  if ( verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ;
337 
338  return NumErrors;
339 }
340 
341 int TestOneMatrix( const std::vector<bool> AmesosClassesInstalled,
342  const char *filename,
343 #ifdef EPETRA_MPI
345 #else
346  Epetra_SerialComm& Comm,
347 #endif
348  // Epetra_Comm &Comm,
349  const bool verbose,
350  const bool PerformDiagonalTest,
351  double Rcond,
352  int &NumTests ) {
353 
354  int NumErrors =0 ;
355  double error = -13;
356  double residual = -13;
357  // double errors[NumAmesosClasses];
358  // double residuals[NumAmesosClasses];
359  // for (int i = 0 ; i < NumAmesosClasses; i ++ ) errors[i] = residuals[i] = 0.0 ;
360 
361  //#ifdef HAVE_AMESOS_UMFPACK
362 #if 0
363  Epetra_CrsMatrix *Amat ;
364 
365  //
366  // Compute the reciprocal condition number using Amesos_UMFPACK via the Amesos interface
367  //
368  Epetra_Map *readMap = 0 ;
369  CreateCrsMatrix( filename, Comm, readMap, false, false, symmetric, Amat ) ;
370  Teuchos::ParameterList ParamList ;
371  Epetra_LinearProblem Problem;
372  Amesos Afactory;
373 
374  Amesos_BaseSolver* Abase ;
375  Abase = Afactory.Create( "Amesos_Umfpack", Problem ) ;
376  if ( Abase == 0 ) {
377  std::cerr << " AMESOS_UMFPACK is required for this test " << std::endl ;
378  exit(13);
379  } ;
380 
381  //
382  // Factor A to compute Rcond = reciprocal condition number estimate
383  //
384  Problem.SetOperator( Amat );
387  Amesos_Umfpack* UmfpackOperator = dynamic_cast<Amesos_Umfpack *> (Abase) ;
388  // double Rcond = UmfpackOperator->GetRcond();
389 
390  int ind[1];
391  double val[1];
392  ind[0] = 0;
393  val[0] = 1 ;
394  double AnormInf = Amat->NormInf() ;
395  if (verbose) std::cout << " norm(Amat) = " << AnormInf << std::endl;
396  if ( Amat->MyGRID( 0 ) )
397  Amat->SumIntoMyValues( 0, 1, val, ind ) ;
398  AnormInf = Amat->NormInf() ;
399  if (verbose) std::cout << " norm(Amat) = " << AnormInf << std::endl;
400 
401 
404  double Rcond1 = UmfpackOperator->GetRcond();
405 
406  if ( Amat->MyGRID( 0 ) )
407  Amat->SumIntoMyValues( 0, 1, val, ind ) ;
408  AnormInf = Amat->NormInf() ;
409  if (verbose) std::cout << " norm(Amat) = " << AnormInf << std::endl;
412  double Rcond2 = UmfpackOperator->GetRcond();
413 
414  if (verbose) std::cout << " Rcond1 = " << Rcond1 << std::endl;
415  if (verbose) std::cout << " Rcond2 = " << Rcond2 << std::endl;
416 
417  if ( readMap ) delete readMap ;
418 #else
419  double Rcond1 = Rcond ;
420  double Rcond2 = Rcond ;
421 #endif
422 
423  //
424  // Rowindex and Colindex control the maps and indices used to create the matrix
425  //
426  // These tests are all disabled in TestAllClasses.cpp
427  //
428  const int RowindexMax = 3; // bug should be three ( 1 based, 3 based, non contiguous )
429  const int ColindexMax = 2; // bug should be two: ( row map, 4 based )
430 
431  //
432  // Rangemap and Domainmap control the Range and Domain maps used in the call to FillComplete
433  // If both are "no change", FillComplete is called with no parameters (i.e. without specifying maps)
434  // Otherwise, domain and range maps are specified in the call to FillComplete
435  //
436  // These tests are all disabled in TestAllClasses.cpp
437  //
438  int RangemapMax = 4; // bug should be four: ( no change, serial, bizarre dist, replicated )
439 
440  //
441  // DiagonalOpts controls whether diagonal elements are left alone,
442  // or removed from both the matrix and the underlying map
443  //
444  int DiagonalOptsMax = 2; // should be two: ( no change, elements missing from map )
445  //
446  //
447  //
448  int EpetraMatrixTypeMax = 3; // 0 = Epetra_CrsMatrix; 1 = Epetra_RowMatriw; 2 = StorageOptimized Epetra_CrsMatrix
449  //
450  // No point in trying to run distributed memory tests on a serial run
451  //
452  int iterDistMax = 2;
453  if ( Comm.NumProc() == 1 ) {
454  iterDistMax = 1 ;
455  RangemapMax = 1 ;
456  }
457  else{
458  RangemapMax = 4;
459  }
460 
461 
462 
463  if (! PerformDiagonalTest ) DiagonalOptsMax = 1 ;
464 
465  for ( int iterTrans =0 ; iterTrans < 2; iterTrans++ ) {
466  bool transpose = iterTrans == 1 ;
467 
468  for ( int iterDist =0; iterDist < iterDistMax; iterDist++ ) {
469  bool distribute = ( iterDist == 1 );
470 
471 #if 1
472  for ( int iterRowindex = 0 ; iterRowindex < RowindexMax; iterRowindex++ ) {
473  for ( int iterColindex = 0 ; iterColindex < ColindexMax; iterColindex++ ) {
474  //
475  // The current version of NewMatNewMap.cpp supports only trivial
476  // replicated maps, hence we do not allow any fancy indexing
477  //
478  int ThisRangemapMax = RangemapMax ;
479  // Bug #1920 Amesos classes can't handle replicated domain or ranges if ( iterRowindex > 0 || iterColindex > 0 )
480  ThisRangemapMax = EPETRA_MIN( 3, ThisRangemapMax );
481  int ThisDomainmapMax = EPETRA_MIN( 3, ThisRangemapMax ); // Bug #1920
482  for ( int iterRangemap = 0 ; iterRangemap < ThisRangemapMax; iterRangemap++ ) {
483  for ( int iterDomainmap = 0 ; iterDomainmap < ThisDomainmapMax; iterDomainmap++ ) {
484  for ( int iterDiagonalOpts = 0 ; iterDiagonalOpts < DiagonalOptsMax; iterDiagonalOpts++ ) {
485 #else
486  int iterRowindex = 0; {
487  int iterColindex = 0 ; {
488  int iterRangemap = 0 ; {
489  int iterDomainmap = 0 ; {
490  for ( int iterDiagonalOpts = 1 ; iterDiagonalOpts < DiagonalOptsMax; iterDiagonalOpts++ ) {
491 #endif
492  const bool printit = false;
493  // diagonal opts testing only works on distributed matrices whose row and column indices match
494  // On a serial matrix, eliminate a column from the map makes the matrix singular
495  // If the row and column indices don't match, eliminating a column from the map is, typically, irrelevant
496 
497  if ( ( iterColindex == 0 && distribute ) || iterDiagonalOpts == 0 ) {
498  if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ;
499  for ( int EpetraMatrixType = 0 ; EpetraMatrixType < EpetraMatrixTypeMax; EpetraMatrixType++ ) {
500 
501  if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ;
502  // These tests presently take over 7 hours on some platforms. But, I don't want to eliminate any category of tests
503  // The following test will cull 90% of the tests and still cover every type of test and most combinations
504  Epetra_Util EU;
505  int RunTest[1] ;
506  RunTest[0] = (EU.RandomDouble() > 0.8)?1:0 ;
507  if ( iterRowindex == 0 &&
508  iterColindex == 0 &&
509  iterRangemap == 0 &&
510  iterDomainmap == 0 ) RunTest[0] = 1;
511  Comm.Broadcast( RunTest, 1, 0 ) ;
512  if ( RunTest[0] ) {
513  //
514  // We test only one level for different indexing or different Range and Domain maps
515  // to avoid hassles of moving data from the domain space to the range space and back
516  //
517  if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ;
518  int MaxLevel = 3 ;
519  if ( iterRowindex > 0 ) MaxLevel = 1 ;
520  if ( iterColindex > 0 ) MaxLevel = 1 ;
521  if ( iterRangemap > 0 ) MaxLevel = 1 ;
522  if ( iterDomainmap > 0 ) MaxLevel = 1 ;
523 
524  bool symmetric = true;
525 
526  Epetra_CrsMatrix *Amat = 0 ;
527  Epetra_Map *readMap = 0 ;
528  CreateCrsMatrix( filename, Comm, readMap, transpose, distribute, symmetric, Amat ) ;
529  // assert( symmetric == false ) ;
530 
531  if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ <<
532  " Creating matrix from " <<
533  " filename = " << filename <<
534  " symmetric = " << symmetric <<
535  " distribute = " << distribute <<
536  " iterRowindex = " << iterRowindex <<
537  " iterColindex = " << iterColindex <<
538  " iterRangemap = " << iterRangemap <<
539  " iterDomainmap = " << iterDomainmap <<
540  " EpetraMatrixType = " << EpetraMatrixType <<
541  " iterDiagonalOpts = " << iterDiagonalOpts <<
542  " transpose = " << transpose << " iterDist = " << iterDist << std::endl ;
543 
544 
545  if ( iterDiagonalOpts ) Comm.SetTracebackMode(1); // Turn off positive Epetra warnings (i.e. iniefficient code, such as memory re-allocation)
546 
547  if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ;
548 
549  RCP<Epetra_CrsMatrix> Bmat = NewMatNewMap( *Amat,
550  iterDiagonalOpts,
551  iterRowindex,
552  iterColindex,
553  iterRangemap,
554  iterDomainmap
555  ) ;
556  if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ;
557  Comm.SetTracebackMode(2);
558 
559  if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ <<
560  " filename = " << filename <<
561  " symmetric = " << symmetric <<
562  " distribute = " << distribute <<
563  " iterRowindex = " << iterRowindex <<
564  " iterColindex = " << iterColindex <<
565  " iterRangemap = " << iterRangemap <<
566  " iterDomainmap = " << iterDomainmap <<
567  " EpetraMatrixType = " << EpetraMatrixType <<
568  " iterDiagonalOpts = " << iterDiagonalOpts <<
569  " transpose = " << transpose << " iterDist = " << iterDist << std::endl ;
570 
571  //
572  // This causes a failure in Amesos_Superludist:
573  Epetra_CrsMatrix* Cmat = &*Bmat;
574  // Epetra_CrsMatrix* Cmat = Amat ;
575 
576 
577  int Level ;
578  double MaxError ;
579  if ( Rcond*Rcond1*Rcond2 > 1e-16 ) {
580  Level = EPETRA_MIN( 3, MaxLevel );
581  MaxError = Rcond*Rcond1*Rcond2;
582  } else if ( Rcond*Rcond1 > 1e-16 ) {
583  Level = EPETRA_MIN( 2, MaxLevel );
584  MaxError = Rcond*Rcond1;
585  } else {
586  Level = EPETRA_MIN( 1, MaxLevel );
587  MaxError = Rcond;
588  }
589 
590  int NumTheseTests = 0 ;
591  if ( verbose ) {
592  std::cout << " About to test " << filename
593  << __FILE__ << "::" << __LINE__
594  << " EpetraMatrixType = " << EpetraMatrixType
595  << (transpose?" transpose":"" )
596  << (distribute?" distribute":"" )
597  << std::endl ;
598  }
599  if ( iterDiagonalOpts == 0 )
600  Comm.SetTracebackMode(2);
601  else
602  Comm.SetTracebackMode(1); // In PerformOneSolveAndTest, MyMatWithDiag->ReplaceDiagonalValues may return 1 indicating that structurally non-zero elements were left untouched.
603 
604  int Error = TestAllClasses( AmesosClasses, EpetraMatrixType,
605  AmesosClassesInstalled,
606  Cmat,
607  transpose ,
608  verbose,
609  symmetric,
610  Level,
611  MaxError,
612  iterDiagonalOpts,
613  iterRowindex,
614  iterColindex,
615  iterRangemap,
616  iterDomainmap,
617  distribute,
618  filename,
619  error,
620  residual,
621  NumTheseTests ) ;
622  NumTests += NumTheseTests ;
623  NumErrors += Error ;
624  if ( Comm.MyPID() == 0 && ( ( verbose && NumTheseTests ) || Error ) ) {
625  std::cout << " Tested " << filename
626  << __FILE__ << "::" << __LINE__
627  << " EpetraMatrixType = " << EpetraMatrixType
628  << (transpose?" transpose":"" )
629  << (distribute?" distribute":"" ) << " error = "
630  << error
631  << " residual = "
632  << residual
633  << std::endl ;
634  }
635  // BUG: Memory leak
636  // delete &(Amat->RowMap()) ;
637  if ( Amat ) delete Amat ;
638  if ( readMap ) delete readMap ;
639 #if 0
640  double relresidual =
641  errors[(int) AMESOS_SUPERLUDIST] = EPETRA_MAX( errors[ (int) AMESOS_SUPERLUDIST], error ) ;
642  residuals[(int) AMESOS_SUPERLUDIST] = EPETRA_MAX( residuals[ (int) AMESOS_SUPERLUDIST], residual ) ;
643  NumErrors += ( residual > maxresidual ) ;
644 #endif
645  }
646  }
647  if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ;
648  }
649  if ( printit && verbose ) std::cout << __FILE__ << "::" << __LINE__ << std::endl ;
650  }
651  }
652  }
653  }
654  }
655  }
656  }
657 
658  return NumErrors;
659 }
660 
661 #if 0
662 #define TEST_P(variable) { { \
663  if ( true ) { std::cerr << "AMESOS_PRINT " << # variable << "= " << variable << std::endl; }; }\
664  }
665 
666 
667 #define TEST_PRINT(variable) { { \
668  if ( true ) { std::cerr << "AMESOS_PRINT " << # variable << "= " << variable << ", " \
669  << __FILE__ << ", line " << __LINE__ << std::endl; }; }\
670  }
671 
672 #endif
673 
674 //
675 // Usage: TestOptions [-s] [-v] [-q]
676 //
677 // -s = short
678 // -v = verbose
679 // -q = quiet
680 //
681 
682 int NextMain( int argc, char *argv[] ) {
683 
684 #ifdef EPETRA_MPI
685  MPI_Init(&argc,&argv);
686  Epetra_MpiComm Comm( MPI_COMM_WORLD );
687 #else
688  Epetra_SerialComm Comm;
689 #endif
690 
691 
692  bool verbose = false;
693  bool small = false ;
694  bool quiet = false ;
695  if ( argc >= 2 && (argv[1][0] == '-') && (argv[1][1] == 'v') )
696  verbose = true ;
697  if ( argc >= 3 && (argv[2][0] == '-') && (argv[2][1] == 'v') )
698  verbose = true ;
699  if ( argc >= 4 && (argv[3][0] == '-') && (argv[3][1] == 'v') )
700  verbose = true ;
701 
702  if ( argc >= 2 && (argv[1][0] == '-') && (argv[1][1] == 's') )
703  small = true ;
704  if ( argc >= 3 && (argv[2][0] == '-') && (argv[2][1] == 's') )
705  small = true ;
706  if ( argc >= 4 && (argv[3][0] == '-') && (argv[3][1] == 's') )
707  small = true ;
708 
709  if ( argc >= 2 && (argv[1][0] == '-') && (argv[1][1] == 'q') )
710  quiet = true ;
711  if ( argc >= 3 && (argv[2][0] == '-') && (argv[2][1] == 'q') )
712  quiet = true ;
713  if ( argc >= 4 && (argv[3][0] == '-') && (argv[3][1] == 'q') )
714  quiet = true ;
715 
716 
717  if ( argc >= 2 && (argv[1][0] == '-') && (argv[1][1] == 'h') ) {
718  std::cerr << "Usage TestOptions [-s] [-v] [-q] " << std::endl ;
719  std::cerr << "-v: verbose " << std::endl ;
720  std::cerr << "-s: small " << std::endl ;
721  std::cerr << "-q: quiet " << std::endl ;
722  exit(-1);
723  }
724 
725 
726 
727 #ifdef HAVE_AMESOS_KLU
728  AmesosClasses.push_back( "Amesos_Klu" );
729 #endif
730 
731 #if 1
732 
733 #ifdef HAVE_AMESOS_PARAKLETE
734  AmesosClasses.push_back( "Amesos_Paraklete" );
735 #endif
736 
737 
738 
739 
740 
741 #ifdef HAVE_AMESOS_PARDISO
742  // bug #1915
743  // bug #1998 - Enabling Amesos_Pardiso causes Amesos_Klu to fail - strange but true
744  // AmesosClasses.push_back( "Amesos_Pardiso" );
745 #endif
746 #ifdef HAVE_AMESOS_SUPERLUDIST
747  AmesosClasses.push_back( "Amesos_Superludist" );
748 #endif
749 
750 
751 
752 
753 #ifdef HAVE_AMESOS_LAPACK
754  AmesosClasses.push_back( "Amesos_Lapack" );
755 #endif
756 
757 #ifdef HAVE_AMESOS_SUPERLU
758  AmesosClasses.push_back( "Amesos_Superlu" );
759 #endif
760 
761 #ifdef HAVE_AMESOS_TAUCS
762  AmesosClasses.push_back( "Amesos_Taucs" );
763 #endif
764 
765 #ifdef HAVE_AMESOS_UMFPACK
766  AmesosClasses.push_back( "Amesos_Umfpack" );
767 #endif
768 
769 
770 
771 #ifdef HAVE_AMESOS_DSCPACK
772  if ( ! quiet ) AmesosClasses.push_back( "Amesos_Dscpack" ); // bug #1205
773 #endif
774 
775 #ifdef HAVE_AMESOS_MUMPS
776  AmesosClasses.push_back( "Amesos_Mumps" );
777 #endif
778 
779 #ifdef HAVE_AMESOS_SCALAPACK
780  AmesosClasses.push_back( "Amesos_Scalapack" ) ;
781 #endif
782 
783 
784 
785 #endif
786 
788  std::vector<bool> AmesosClassesInstalled( NumAmesosClasses );
789 
790  assert( NumAmesosClasses > 0 ) ;
791 
792 
793  if ( Comm.MyPID() != 0 ) verbose = false ;
794 #if 0
795  //
796  // Wait for a character to allow time to attach the debugger
797  //
798  if ( Comm.MyPID() == 0 ) {
799  char what = 'a';
800  while ( what == 'a' ) // I don't know why we need this while loop at least on bewoulf
801  std::cin >> what ;
802  }
803  Comm.Barrier();
804 
805  std::cout << __FILE__ << "::" << __LINE__ << " Comm.MyPID() = " << Comm.MyPID() << std::endl ;
806 #endif
807 
808 
809 
810  Epetra_LinearProblem Problem;
811  Amesos_BaseSolver* Abase ;
812  Amesos Afactory;
813 
814  Comm.SetTracebackMode(2);
815 
816 #ifndef HAVE_AMESOS_EPETRAEXT
817  if ( ! quiet && Comm.MyPID() == 0 )
818  std::cout << "Amesos has been built without epetraext, capabilites requiring epetraext, such as reindexing and Amesos_Paraklete non-transpose solves, will not be tested" << std::endl ;
819 #endif
820 
821  for (int i=0; i < NumAmesosClasses; i++ ) {
822 
823 
824  Abase = Afactory.Create( &AmesosClasses[i][0], Problem ) ;
825  if ( Abase == 0 ) {
826  if ( !quiet && Comm.MyPID() == 0 ) std::cout << AmesosClasses[i] << " not built in this configuration" << std::endl ;
827  AmesosClassesInstalled[i] = false;
828  } else {
829  if ( !quiet && Comm.MyPID() == 0 ) std::cout << " Testing " << AmesosClasses[i] << std::endl ;
830  AmesosClassesInstalled[i] = true;
831  Teuchos::ParameterList ParamList ;
832  ParamList.set( "NoDestroy", true ); // Prevents Amesos_Mumps from deleting data
833  Abase->SetParameters( ParamList ); // which causes Amesos_Mumps to crash on this trivial instantiation
834  }
835  delete Abase ;
836  }
837 
838  int result = 0 ;
839  int numtests = 0 ;
840 
841  // ImpcolB.rua fails - the failure could be in the test code, in particular in NewMatNewMap.cpp
842  // result += TestOneMatrix( AmesosClassesInstalled, "../Test_Basic/ImpcolB.rua", Comm, verbose, false, 1e-9 , numtests ) ;
843 
844  //
845  // Khead.triS remains non-singular even after a diagaonal element is removed from the map
846  //
847  // Khead.triS fails on DSCPACK
848  result += TestOneMatrix( AmesosClassesInstalled, (char *) "../Test_Basic/SuperLU.triU", Comm, verbose, false, 1e-6 , numtests ) ;
849 
850  //
851  // small is set by TestValgrind - keep testing to a minimum because execution time is so slow
852  // quiet is set by TestQuietAmesos - dscpack is not quiet at the moment, hence we can't test symmetric matrices
853  // in TestQuietAmesos
854  //
855  // bug #1205
856  //
857  if ( ! small ) {
858  result += TestOneMatrix( AmesosClassesInstalled, (char *) "../Test_Basic/MissingADiagonal.mtx", Comm, verbose, false, 1e-2 , numtests ) ;
859  result += TestOneMatrix( AmesosClassesInstalled, "../Test_Basic/Khead.triS", Comm, verbose, true, 1e-6 , numtests ) ;
860  result += TestOneMatrix( AmesosClassesInstalled, (char *) "../Test_Basic/bcsstk04.mtx", Comm, verbose, false, 1e-4 , numtests ) ;
861  //
862  // The file reader for .rua files is not quiet
863  //
864  if (! quiet) {
865  result += TestOneMatrix( AmesosClassesInstalled, (char *) "../Test_Basic/Diagonal.mtx", Comm, verbose, false, 1e-1 , numtests ) ;
866  if ( Comm.NumProc() == 1) {
867  result += TestOneMatrix( AmesosClassesInstalled, (char *) "../Test_Basic/662_bus_out.rsa", Comm, verbose, false, 1e-5 , numtests ) ;
868  result += TestOneMatrix( AmesosClassesInstalled, (char *) "../Test_Basic/SuperLU.rua", Comm, verbose, false, 1e-2 , numtests ) ;
869  result += TestOneMatrix( AmesosClassesInstalled, "../Test_Basic/ImpcolB.rua", Comm, verbose, false, 1e-6 , numtests ) ;
870  // result += TestOneMatrix( AmesosClassesInstalled, "../Test_Basic/ImpcolB.rua", Comm, verbose, false, 1e-6 , numtests ) ;
871  }
872  }
873  }
874 
875  // bug #2184 - Amesos_Klu fails to detect Structurally singular matrices
876  // This test, as modified in PerformOneSolveAndTest.cpp, tests the present
877  // beahviour - i.e. that SymbolicFactorization() fails to detect the
878  // structurally singular matrix, but that NumericFactorization()
879  // catches the singularity instead.
880  result+=TestErrors( AmesosClassesInstalled, (char *) "../Test_Basic/StructurallySingular.mtx",
881  Comm, verbose, numtests ) ;
882  result+=TestErrors( AmesosClassesInstalled, (char *) "../Test_Basic/NumericallySingular.mtx",
883  Comm, verbose, numtests ) ;
884 
885  if ( ! quiet && Comm.MyPID() == 0 ) std::cout << result << " Tests failed " << numtests << " Tests performed " << std::endl ;
886 
887  if ( result == 0 && numtests > 0 ) {
888  if (! quiet && Comm.MyPID() == 0)
889  std::cout << std::endl << "TEST PASSED" << std::endl << std::endl;
890  }
891  else {
892  if (Comm.MyPID() == 0)
893  std::cout << std::endl << "TEST FAILED" << std::endl << std::endl;
894  AMESOS_CHK_ERR( 1 ) ;
895  }
896 
897 #ifdef EPETRA_MPI
898  MPI_Finalize();
899 #endif
900 
901  return result ;
902 }
903 
904 //
905 // I put this in hoping that this would eliminate a bogus memory leak report
906 // from valgrind.
907 //
908 int main( int argc, char *argv[] ) {
909  int retval = NextMain( argc, argv ) ;
910  return retval ;
911 }
int NumGlobalElements() const
bool MyGRID(int GRID_in) const
void SetOperator(Epetra_RowMatrix *A)
int NextMain(int argc, char *argv[])
bool GlobalIndicesLongLong() const
bool SameAs(const Epetra_BlockMap &Map) const
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
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.
double NormInf() const
static bool verbose
Definition: Amesos.cpp:67
int CrsMatrixTranspose(Epetra_CrsMatrix *In, Epetra_CrsMatrix *Out)
#define EPETRA_MIN(x, y)
bool GlobalIndicesInt() const
int MyPID() const
int FillComplete(bool OptimizeDataStorage=true)
int TestErrors(const std::vector< bool > AmesosClassesInstalled, const char *filename, Epetra_MpiComm &Comm, const bool verbose, int &NumTests)
const Epetra_Map & RowMap() const
std::string filename
int TestAllClasses(const std::vector< std::string > AmesosClasses, int EpetraMatrixType, const std::vector< bool > AmesosClassesInstalled, Epetra_CrsMatrix *&Amat, const bool transpose, const bool verbose, const bool symmetric, const int Levels, const double Rcond, int Diagonal, int ReindexRowMap, int ReindexColMap, int RangeMapType, int DomainMapType, bool distribute, const char *filename, double &maxrelerror, double &maxrelresidual, int &NumTests)
#define AMESOS_CHK_ERR(a)
int main(int argc, char *argv[])
Factory for binding a third party direct solver to an Epetra_LinearProblem.
Definition: Amesos.h:44
#define BUFSIZE
#define EPETRA_CHK_ERR(xxx)
int NumProc() const
std::string error
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 NumAmesosClasses
Definition: TestOptions.cpp:79
RCP< Epetra_CrsMatrix > NewMatNewMap(Epetra_CrsMatrix &In, int Diagonal, int ReindexRowMap, int ReindexColMap, int RangeMapType, int DomainMapType)
int SumIntoMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.
Definition: Amesos.cpp:69
void Barrier() const
int TestOneMatrix(const std::vector< bool > AmesosClassesInstalled, const char *filename, Epetra_MpiComm &Comm, const bool verbose, const bool PerformDiagonalTest, double Rcond, int &NumTests)
double RandomDouble()
const Epetra_CrsGraph & Graph() const
Amesos_BaseSolver: A pure virtual class for direct solution of real-valued double-precision operators...
Class Amesos_Umfpack: An object-oriented wrapper for UMFPACK.
#define EPETRA_MAX(x, y)
std::vector< string > AmesosClasses
Definition: TestOptions.cpp:77