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