38 #include "Trilinos_Util.h" 39 #include "Trilinos_Util_ReadMatrixMarket2Epetra.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" 49 #ifdef HAVE_AMESOS_UMFPACK 58 #include "Epetra_MpiComm.h" 60 #include "Epetra_SerialComm.h" 65 #ifdef HAVE_VALGRIND_H 66 #include <valgrind/valgrind.h> 69 #ifdef HAVE_VALGRIND_VALGRIND_H 70 #include <valgrind/valgrind.h> 82 Epetra_Map *& readMap,
83 const bool transpose,
const bool distribute,
84 bool& symmetric, Epetra_CrsMatrix *& Matrix ) {
86 Epetra_CrsMatrix * readA = 0;
87 Epetra_Vector * readx = 0;
88 Epetra_Vector * readb = 0;
89 Epetra_Vector * readxexact = 0;
95 FILE *in_file = fopen( in_filename,
"r");
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 );
113 if ( LastFiveBytes ==
".triU" ) {
116 readb, readxexact) );
119 if ( LastFiveBytes ==
".triS" ) {
122 readb, readxexact) );
125 if ( LastFourBytes ==
".mtx" ) {
127 readA, readx, readb, readxexact) );
129 assert (in_file !=
NULL) ;
130 const int BUFSIZE = 800 ;
131 char buffer[BUFSIZE] ;
132 fgets( buffer, BUFSIZE, in_file ) ;
133 std::string headerline1 = buffer;
135 if ( headerline1.find(
"symmetric") < BUFSIZE ) symmetric =
true;
137 if ( headerline1.find(
"symmetric") != std::string::npos) symmetric =
true;
144 Trilinos_Util_ReadHb2Epetra(
filename,
Comm, readMap, readA, readx,
146 if ( LastFourBytes ==
".rsa" ) symmetric = true ;
152 if ( readb )
delete readb;
153 if ( readx )
delete readx;
154 if ( readxexact )
delete readxexact;
156 Epetra_CrsMatrix *serialA ;
157 Epetra_CrsMatrix *transposeA;
162 transposeA =
new Epetra_CrsMatrix(
Copy, *readMap, 0 );
165 serialA = transposeA ;
172 assert( (
void *) &serialA->Graph() ) ;
173 assert( (
void *) &serialA->RowMap() ) ;
174 assert( serialA->RowMap().SameAs(*readMap) ) ;
178 Epetra_Map* mapPtr = 0;
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);
187 Epetra_Map& DistMap = *mapPtr;
190 Epetra_Export exporter( *readMap, DistMap );
192 Epetra_CrsMatrix *Amat =
new Epetra_CrsMatrix(
Copy, DistMap, 0 );
193 Amat->Export(*serialA, exporter, Add);
194 ierr = Amat->FillComplete();
218 int TestErrors(
const std::vector<bool> AmesosClassesInstalled,
219 const char *filename,
221 Epetra_MpiComm&
Comm,
223 Epetra_SerialComm&
Comm,
230 double residual = -13;
233 for (
int iterTrans =0 ; iterTrans < 2; iterTrans++ ) {
236 const bool distribute = 1;
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 ;
249 Epetra_CrsMatrix *Amat = 0 ;
250 Epetra_Map *readMap = 0 ;
253 if ( printit &&
verbose ) std::cout << __FILE__ <<
"::" << __LINE__ <<
254 " Creating matrix from " <<
256 " symmetric = " << symmetric <<
257 " distribute = " << distribute <<
258 " iterRowindex = " << iterRowindex <<
259 " iterColindex = " << iterColindex <<
260 " iterRangemap = " << iterRangemap <<
261 " iterDomainmap = " << iterDomainmap <<
262 " EpetraMatrixType = " << EpetraMatrixType <<
263 " iterDiagonalOpts = " << iterDiagonalOpts <<
265 <<
" iterDist = " << iterDist << std::endl ;
267 if ( iterDiagonalOpts )
Comm.SetTracebackMode(1);
269 if ( printit &&
verbose ) std::cout << __FILE__ <<
"::" << __LINE__ << std::endl ;
278 if ( printit &&
verbose ) std::cout << __FILE__ <<
"::" << __LINE__ << std::endl ;
279 Comm.SetTracebackMode(2);
281 if ( printit &&
verbose ) std::cout << __FILE__ <<
"::" << __LINE__ <<
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 ;
295 Epetra_CrsMatrix* Cmat = &*Bmat;
300 const double MaxError = 1e-3;
302 int NumTheseTests = 0 ;
304 std::cout <<
" About to test " <<
filename 305 << __FILE__ <<
"::" << __LINE__
306 <<
" EpetraMatrixType = " << EpetraMatrixType
308 << (distribute?
" distribute":
"" )
312 AmesosClassesInstalled,
329 NumTests += NumTheseTests ;
333 if ( Amat )
delete Amat ;
334 if ( readMap )
delete readMap ;
336 if (
verbose ) std::cout << __FILE__ <<
"::" << __LINE__ << std::endl ;
342 const char *filename,
344 Epetra_MpiComm&
Comm,
346 Epetra_SerialComm&
Comm,
350 const bool PerformDiagonalTest,
356 double residual = -13;
363 Epetra_CrsMatrix *Amat ;
368 Epetra_Map *readMap = 0 ;
371 Epetra_LinearProblem Problem;
375 Abase = Afactory.
Create(
"Amesos_Umfpack", Problem ) ;
377 std::cerr <<
" AMESOS_UMFPACK is required for this test " << std::endl ;
384 Problem.SetOperator( Amat );
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;
404 double Rcond1 = UmfpackOperator->GetRcond();
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();
414 if (
verbose) std::cout <<
" Rcond1 = " << Rcond1 << std::endl;
415 if (
verbose) std::cout <<
" Rcond2 = " << Rcond2 << std::endl;
417 if ( readMap )
delete readMap ;
419 double Rcond1 = Rcond ;
420 double Rcond2 = Rcond ;
428 const int RowindexMax = 3;
429 const int ColindexMax = 2;
444 int DiagonalOptsMax = 2;
448 int EpetraMatrixTypeMax = 3;
453 if (
Comm.NumProc() == 1 ) {
463 if (! PerformDiagonalTest ) DiagonalOptsMax = 1 ;
465 for (
int iterTrans =0 ; iterTrans < 2; iterTrans++ ) {
468 for (
int iterDist =0; iterDist < iterDistMax; iterDist++ ) {
469 bool distribute = ( iterDist == 1 );
472 for (
int iterRowindex = 0 ; iterRowindex < RowindexMax; iterRowindex++ ) {
473 for (
int iterColindex = 0 ; iterColindex < ColindexMax; iterColindex++ ) {
478 int ThisRangemapMax = RangemapMax ;
480 ThisRangemapMax = EPETRA_MIN( 3, ThisRangemapMax );
481 int ThisDomainmapMax = EPETRA_MIN( 3, ThisRangemapMax );
482 for (
int iterRangemap = 0 ; iterRangemap < ThisRangemapMax; iterRangemap++ ) {
483 for (
int iterDomainmap = 0 ; iterDomainmap < ThisDomainmapMax; iterDomainmap++ ) {
484 for (
int iterDiagonalOpts = 0 ; iterDiagonalOpts < DiagonalOptsMax; iterDiagonalOpts++ ) {
486 int iterRowindex = 0; {
487 int iterColindex = 0 ; {
488 int iterRangemap = 0 ; {
489 int iterDomainmap = 0 ; {
490 for (
int iterDiagonalOpts = 1 ; iterDiagonalOpts < DiagonalOptsMax; iterDiagonalOpts++ ) {
492 const bool printit =
false;
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++ ) {
501 if ( printit &&
verbose ) std::cout << __FILE__ <<
"::" << __LINE__ << std::endl ;
506 RunTest[0] = (EU.RandomDouble() > 0.8)?1:0 ;
507 if ( iterRowindex == 0 &&
510 iterDomainmap == 0 ) RunTest[0] = 1;
511 Comm.Broadcast( RunTest, 1, 0 ) ;
517 if ( printit &&
verbose ) std::cout << __FILE__ <<
"::" << __LINE__ << std::endl ;
519 if ( iterRowindex > 0 ) MaxLevel = 1 ;
520 if ( iterColindex > 0 ) MaxLevel = 1 ;
521 if ( iterRangemap > 0 ) MaxLevel = 1 ;
522 if ( iterDomainmap > 0 ) MaxLevel = 1 ;
524 bool symmetric =
true;
526 Epetra_CrsMatrix *Amat = 0 ;
527 Epetra_Map *readMap = 0 ;
531 if ( printit &&
verbose ) std::cout << __FILE__ <<
"::" << __LINE__ <<
532 " Creating matrix from " <<
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 ;
545 if ( iterDiagonalOpts )
Comm.SetTracebackMode(1);
547 if ( printit &&
verbose ) std::cout << __FILE__ <<
"::" << __LINE__ << std::endl ;
556 if ( printit &&
verbose ) std::cout << __FILE__ <<
"::" << __LINE__ << std::endl ;
557 Comm.SetTracebackMode(2);
559 if ( printit &&
verbose ) std::cout << __FILE__ <<
"::" << __LINE__ <<
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 ;
573 Epetra_CrsMatrix* Cmat = &*Bmat;
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;
586 Level = EPETRA_MIN( 1, MaxLevel );
590 int NumTheseTests = 0 ;
592 std::cout <<
" About to test " <<
filename 593 << __FILE__ <<
"::" << __LINE__
594 <<
" EpetraMatrixType = " << EpetraMatrixType
596 << (distribute?
" distribute":
"" )
599 if ( iterDiagonalOpts == 0 )
600 Comm.SetTracebackMode(2);
602 Comm.SetTracebackMode(1);
605 AmesosClassesInstalled,
622 NumTests += NumTheseTests ;
624 if (
Comm.MyPID() == 0 && ( (
verbose && NumTheseTests ) || Error ) ) {
626 << __FILE__ <<
"::" << __LINE__
627 <<
" EpetraMatrixType = " << EpetraMatrixType
629 << (distribute?
" distribute":
"" ) <<
" error = " 637 if ( Amat )
delete Amat ;
638 if ( readMap )
delete readMap ;
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 ) ;
647 if ( printit &&
verbose ) std::cout << __FILE__ <<
"::" << __LINE__ << std::endl ;
649 if ( printit &&
verbose ) std::cout << __FILE__ <<
"::" << __LINE__ << std::endl ;
662 #define TEST_P(variable) { { \ 663 if ( true ) { std::cerr << "AMESOS_PRINT " << # variable << "= " << variable << std::endl; }; }\ 667 #define TEST_PRINT(variable) { { \ 668 if ( true ) { std::cerr << "AMESOS_PRINT " << # variable << "= " << variable << ", " \ 669 << __FILE__ << ", line " << __LINE__ << std::endl; }; }\ 685 MPI_Init(&argc,&argv);
686 Epetra_MpiComm
Comm( MPI_COMM_WORLD );
688 Epetra_SerialComm
Comm;
695 if ( argc >= 2 && (argv[1][0] ==
'-') && (argv[1][1] ==
'v') )
697 if ( argc >= 3 && (argv[2][0] ==
'-') && (argv[2][1] ==
'v') )
699 if ( argc >= 4 && (argv[3][0] ==
'-') && (argv[3][1] ==
'v') )
702 if ( argc >= 2 && (argv[1][0] ==
'-') && (argv[1][1] ==
's') )
704 if ( argc >= 3 && (argv[2][0] ==
'-') && (argv[2][1] ==
's') )
706 if ( argc >= 4 && (argv[3][0] ==
'-') && (argv[3][1] ==
's') )
709 if ( argc >= 2 && (argv[1][0] ==
'-') && (argv[1][1] ==
'q') )
711 if ( argc >= 3 && (argv[2][0] ==
'-') && (argv[2][1] ==
'q') )
713 if ( argc >= 4 && (argv[3][0] ==
'-') && (argv[3][1] ==
'q') )
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 ;
727 #ifdef HAVE_AMESOS_KLU 733 #ifdef HAVE_AMESOS_PARAKLETE 741 #ifdef HAVE_AMESOS_PARDISO 746 #ifdef HAVE_AMESOS_SUPERLUDIST 753 #ifdef HAVE_AMESOS_LAPACK 757 #ifdef HAVE_AMESOS_SUPERLU 761 #ifdef HAVE_AMESOS_TAUCS 765 #ifdef HAVE_AMESOS_UMFPACK 771 #ifdef HAVE_AMESOS_DSCPACK 775 #ifdef HAVE_AMESOS_MUMPS 779 #ifdef HAVE_AMESOS_SCALAPACK 798 if (
Comm.MyPID() == 0 ) {
800 while ( what ==
'a' )
805 std::cout << __FILE__ <<
"::" << __LINE__ <<
" Comm.MyPID() = " <<
Comm.MyPID() << std::endl ;
810 Epetra_LinearProblem Problem;
814 Comm.SetTracebackMode(2);
816 #ifndef HAVE_AMESOS_EPETRAEXT 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 ;
826 if ( !
quiet &&
Comm.MyPID() == 0 ) std::cout <<
AmesosClasses[i] <<
" not built in this configuration" << std::endl ;
827 AmesosClassesInstalled[i] =
false;
830 AmesosClassesInstalled[i] =
true;
832 ParamList.
set(
"NoDestroy",
true );
848 result +=
TestOneMatrix( AmesosClassesInstalled, (
char *)
"../Test_Basic/SuperLU.triU",
Comm,
verbose,
false, 1e-6 , numtests ) ;
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 ) ;
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 ) ;
880 result+=
TestErrors( AmesosClassesInstalled, (
char *)
"../Test_Basic/StructurallySingular.mtx",
882 result+=
TestErrors( AmesosClassesInstalled, (
char *)
"../Test_Basic/NumericallySingular.mtx",
885 if ( !
quiet &&
Comm.MyPID() == 0 ) std::cout << result <<
" Tests failed " << numtests <<
" Tests performed " << std::endl ;
887 if ( result == 0 && numtests > 0 ) {
889 std::cout << std::endl <<
"TEST PASSED" << std::endl << std::endl;
892 if (
Comm.MyPID() == 0)
893 std::cout << std::endl <<
"TEST FAILED" << std::endl << std::endl;
908 int main(
int argc,
char *argv[] ) {
909 int retval =
NextMain( argc, argv ) ;
int main(int argc, char *argv[])
int NextMain(int argc, char *argv[])
cholmod_sparse *CHOLMOD() transpose(cholmod_sparse *A, int values, cholmod_common *Common)
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.
int CrsMatrixTranspose(Epetra_CrsMatrix *In, Epetra_CrsMatrix *Out)
int TestErrors(const std::vector< bool > AmesosClassesInstalled, const char *filename, Epetra_MpiComm &Comm, const bool verbose, int &NumTests)
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)
Factory for binding a third party direct solver to an Epetra_LinearProblem.
int CHOLMOD() error(int status, char *file, int line, char *message, cholmod_common *Common)
#define EPETRA_CHK_ERR(xxx)
int CreateCrsMatrix(const char *in_filename, const Epetra_Comm &Comm, Epetra_Map *&readMap, const bool transpose, const bool distribute, bool &symmetric, Epetra_CrsMatrix *&Matrix)
RCP< Epetra_CrsMatrix > NewMatNewMap(Epetra_CrsMatrix &In, int Diagonal, int ReindexRowMap, int ReindexColMap, int RangeMapType, int DomainMapType)
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.
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.
std::vector< string > AmesosClasses