33 #include "Trilinos_Util.h"
42 #ifdef HAVE_AMESOS_DSCPACK
45 #ifdef HAVE_AMESOS_LAPACK
48 #ifdef HAVE_AMESOS_SLUS
51 #ifdef HAVE_AMESOS_UMFPACK
54 #ifdef HAVE_AMESOS_KLU
57 #ifdef HAVE_AMESOS_TAUCS
60 #ifdef HAVE_AMESOS_PARDISO
63 #ifdef HAVE_AMESOS_PARAKLETE
66 #if defined(HAVE_AMESOS_MUMPS) && defined(HAVE_MPI)
69 #ifdef HAVE_AMESOS_SCALAPACK
72 #ifdef HAVE_AMESOS_SLUD
75 #ifdef HAVE_AMESOS_SUPERLU
78 #ifdef HAVE_AMESOS_SUPERLUDIST
81 #ifdef HAVE_AMESOS_SLUD2
127 std::string FileName = matrix_file ;
128 int FN_Size = FileName.size() ;
129 std::string LastFiveBytes = FileName.substr(
EPETRA_MAX(0,FN_Size-5), FN_Size );
130 std::string LastFourBytes = FileName.substr(
EPETRA_MAX(0,FN_Size-4), FN_Size );
131 bool NonContiguousMap =
false;
133 if ( LastFiveBytes ==
".triU" ) {
135 NonContiguousMap =
true;
136 EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( matrix_file,
false, Comm, readMap, readA, readx,
137 readb, readxexact, NonContiguousMap ) );
139 if ( LastFiveBytes ==
".triS" ) {
140 NonContiguousMap =
true;
142 EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( matrix_file,
true, Comm, readMap, readA, readx,
143 readb, readxexact) );
145 if ( LastFourBytes ==
".mtx" ) {
146 EPETRA_CHK_ERR( Trilinos_Util_ReadMatrixMarket2Epetra( matrix_file, Comm, readMap,
147 readA, readx, readb, readxexact) );
150 Trilinos_Util_ReadHb2Epetra( matrix_file, Comm, readMap, readA, readx,
162 serialA = &transposeA ;
172 if( NonContiguousMap ) {
177 int NumMyElements = map.NumMyElements();
178 int MyFirstElement = map.MinMyGID();
179 std::vector<int> MapMap_( NumGlobalElements );
181 Comm.
Broadcast( &MapMap_[0], NumGlobalElements, 0 ) ;
182 map_ =
new Epetra_Map( NumGlobalElements, NumMyElements, &MapMap_[MyFirstElement], 0, Comm);
213 if ( distribute_matrix ) {
218 A.Export(*serialA, exporter,
Add);
227 passxexact = &xexact;
234 passxexact = &serialxexact;
235 passresid = &serialresid;
236 passtmp = &serialtmp;
239 passxexact->SetSeed(131) ;
240 passxexact->Random();
241 passx->SetSeed(11231) ;
244 passb->PutScalar( 0.0 );
245 passA->
Multiply( transpose, *passxexact, *passb ) ;
249 double Anorm = passA->
NormInf() ;
256 double max_resid = 0.0;
257 for (
int j = 0 ; j < special+1 ; j++ ) {
265 }
else if ( SparseSolver ==
UMFPACK ) {
270 umfpack.SetTrans( transpose ) ;
274 }
else if ( SparseSolver == SuperLU ) {
275 SuperluserialOO superluserial ;
279 superluserial.SetTrans( transpose ) ;
280 superluserial.SetUseDGSSV( special == 0 ) ;
282 for (
int i= 0 ; i < numsolves ; i++ ) {
286 superluserial.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
287 superluserial.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
289 superluserial.Solve() ;
293 if ( i < numsolves-1 )
301 #ifdef HAVE_AMESOS_SLUD
302 }
else if ( SparseSolver == SuperLUdist ) {
307 for (
int i= 0 ; i < numsolves ; i++ ) {
311 Problem.
SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
312 Problem.
SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
318 if ( i < numsolves-1 )
326 #ifdef HAVE_AMESOS_SLUD2
327 }
else if ( SparseSolver == SuperLUdist2 ) {
329 superludist2.
SetTrans( transpose ) ;
332 for (
int i= 0 ; i < numsolves ; i++ ) {
336 Problem.
SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
337 Problem.
SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
343 if ( i < numsolves-1 )
351 #ifdef HAVE_AMESOS_DSCPACK
352 }
else if ( SparseSolver ==
DSCPACK ) {
355 ParamList.
set(
"MaxProcs", -3 );
358 for (
int i= 0 ; i < numsolves ; i++ ) {
362 Problem.
SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
363 Problem.
SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
368 if ( i < numsolves-1 )
376 #ifdef HAVE_AMESOS_UMFPACK
377 }
else if ( SparseSolver ==
UMFPACK ) {
380 ParamList.
set(
"MaxProcs", -3 );
384 for (
int i= 0 ; i < numsolves ; i++ ) {
388 Problem.
SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
389 Problem.
SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
394 if ( i < numsolves-1 )
402 #ifdef HAVE_AMESOS_SUPERLU
403 }
else if ( SparseSolver ==
SUPERLU ) {
406 ParamList.
set(
"MaxProcs", -3 );
412 for (
int i= 0 ; i < numsolves ; i++ ) {
416 Problem.
SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
417 Problem.
SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
422 if ( i < numsolves-1 )
430 #ifdef HAVE_AMESOS_SLUS
431 }
else if ( SparseSolver == SuperLU ) {
436 for (
int i= 0 ; i < numsolves ; i++ ) {
440 Problem.
SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
441 Problem.
SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
446 if ( i < numsolves-1 )
454 #ifdef HAVE_AMESOS_KLU
455 }
else if ( SparseSolver ==
KLU ) {
460 ParamList.
set(
"MaxProcs", -3 );
462 ParamList.
set(
"MaxProcs", -3 );
467 for (
int trials = 0 ; trials <= 1 ; trials++) {
469 for (
int i= 0 ; i < numsolves ; i++ ) {
473 Problem.
SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
474 Problem.
SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
481 if ( i < numsolves-1 )
491 #ifdef HAVE_AMESOS_LAPACK
492 }
else if ( SparseSolver ==
LAPACK ) {
499 for (
int i= 0 ; i < numsolves ; i++ ) {
503 Problem.
SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
504 Problem.
SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
509 if ( i < numsolves-1 )
517 #ifdef HAVE_AMESOS_TAUCS
518 }
else if ( SparseSolver ==
TAUCS ) {
521 ParamList.
set(
"MaxProcs", -3 );
527 for (
int i= 0 ; i < numsolves ; i++ ) {
531 Problem.
SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
532 Problem.
SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
537 if ( i < numsolves-1 )
545 #ifdef HAVE_AMESOS_PARDISO
546 }
else if ( SparseSolver ==
PARDISO ) {
549 ParamList.
set(
"MaxProcs", -3 );
555 for (
int i= 0 ; i < numsolves ; i++ ) {
559 Problem.
SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
560 Problem.
SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
565 if ( i < numsolves-1 )
573 #ifdef HAVE_AMESOS_PARAKLETE
574 }
else if ( SparseSolver ==
PARAKLETE ) {
577 ParamList.
set(
"MaxProcs", -3 );
583 for (
int i= 0 ; i < numsolves ; i++ ) {
587 Problem.
SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
588 Problem.
SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
593 if ( i < numsolves-1 )
601 #if defined(HAVE_AMESOS_MUMPS) && defined(HAVE_MPI)
602 }
else if ( SparseSolver ==
MUMPS ) {
605 ParamList.
set(
"MaxProcs", -3 );
611 for (
int i= 0 ; i < numsolves ; i++ ) {
615 Problem.
SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
616 Problem.
SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
621 if ( i < numsolves-1 )
629 #ifdef HAVE_AMESOS_SCALAPACK
630 }
else if ( SparseSolver ==
SCALAPACK ) {
633 ParamList.
set(
"MaxProcs", -3 );
639 for (
int i= 0 ; i < numsolves ; i++ ) {
643 Problem.
SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
644 Problem.
SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
649 if ( i < numsolves-1 )
657 #ifdef HAVE_AMESOS_SUPERLUDIST
660 ParamList.
set(
"MaxProcs", -3 );
668 for (
int i= 0 ; i < numsolves ; i++ ) {
672 Problem.
SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
673 Problem.
SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
675 if ( i < numsolves-1 )
682 }
else if ( SparseSolver == SPOOLES ) {
690 #ifdef TEST_SPOOLESSERIAL
691 }
else if ( SparseSolver == SPOOLESSERIAL ) {
696 spoolesserial.Solve() ;
700 std::cerr <<
"\n\n#################### Requested solver not available (Or not tested with multiple RHS) on this platform #####################\n" << std::endl ;
708 std::vector <double>
error(numsolves) ;
709 double max_error = 0.0;
711 passresid->Update(1.0, *passx, -1.0, *passxexact, 0.0);
713 passresid->Norm2(&error[0]);
714 for (
int i = 0 ; i< numsolves; i++ )
715 if ( error[i] > max_error ) max_error = error[i] ;
724 std::vector <double> residual(numsolves) ;
726 passtmp->PutScalar(0.0);
727 passA->
Multiply( transpose, *passx, *passtmp);
728 passresid->Update(1.0, *passtmp, -1.0, *passb, 0.0);
730 passresid->Norm2(&residual[0]);
732 for (
int i = 0 ; i< numsolves; i++ )
733 if ( residual[i] > max_resid ) max_resid = residual[i] ;
738 std::vector <double> bnorm(numsolves);
739 passb->Norm2( &bnorm[0] ) ;
742 std::vector <double> xnorm(numsolves);
743 passx->Norm2( &xnorm[0] ) ;
int SetUseTranspose(bool UseTranspose)
SetUseTranpose(true) is more efficient in Amesos_Scalapack.
int NumGlobalElements() const
Amesos_Klu: A serial, unblocked code ideal for getting started and for very sparse matrices...
void SetLHS(Epetra_MultiVector *X)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int SetUseTranspose(bool UseTranspose_in)
If set true, X will be set to the solution of AT X = B (not A X = B)
Epetra_SLU: An object-oriented wrapper for Xiaoye Li's serial sparse solver package: Superlu...
int Solve()
Solves A X = B (or AT x = B)
int Solve()
Solves A X = B (or AT x = B)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
void SetTrans(bool trans)
SuperludistOO: An object-oriented wrapper for Xiaoye Li's Superludist.
Amesos_Paraklete: A serial, unblocked code ideal for getting started and for very sparse matrices...
Amesos_Superludist: An object-oriented wrapper for Superludist.
int Solve()
Solves A X = B (or AT x = B)
Amesos_Mumps: An object-oriented wrapper for the double precision version of MUMPS.
int Solve()
Solves A X = B (or AT x = B)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int MyGlobalElements(int *MyGlobalElementList) const
int Solve()
Solves A X = B (or AT x = B)
void Set_Total_Time(double Total_Time_in)
double ElapsedTime(void) const
Amesos_Dscpack: An object-oriented wrapper for Dscpack.
int Amesos_TestMrhsSolver(Epetra_Comm &Comm, char *matrix_file, int numsolves, SparseSolverType SparseSolver, bool transpose, int special, AMESOS_MatrixType matrix_type)
void Set_Anorm(double anorm_in)
int SetUseTranspose(bool UseTranspose)
SetUseTranpose()
int CrsMatrixTranspose(Epetra_CrsMatrix *In, Epetra_CrsMatrix *Out)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
void SetTrans(bool trans)
Setting the transpose flag to true causes Solve() to compute A^t x = b.
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
virtual void Barrier() const =0
Amesos_Pardiso: Interface to the PARDISO package.
Amesos_Superlu: Amesos interface to Xioye Li's SuperLU 3.0 serial code.
int Solve()
Solves A X = B (or AT X = B)
int FillComplete(bool OptimizeDataStorage=true)
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
void Set_Error(double error_in)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int SetUseTranspose(bool UseTranspose_in)
SetUseTranpose()
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int SetUseTranspose(bool UseTranspose_in)
If set true, X will be set to the solution of AT X = B (not A X = B)
int SetUseTranspose(bool UseTranspose_in)
If set true, X will be set to the solution of AT X = B (not A X = B)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int SetParameters(Teuchos::ParameterList &ParameterList)
Set parameters from the input parameters list, returns 0 if successful.
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
virtual double NormInf() const =0
void Set_Middle_Time(double Middle_Time_in)
int Solve(bool Factor)
All computation is performed during the call to Solve()
int NumericFactorization()
Performs NumericFactorization on the matrix A.
#define EPETRA_CHK_ERR(xxx)
int NumericFactorization()
Performs NumericFactorization on the matrix A.
void Set_Xnorm(double xnorm_in)
int SetUseTranspose(bool useTheTranspose)
If set true, X will be set to the solution of AT X = B (not A X = B)
int Solve()
Solves A X = B (or AT x = B)
virtual int Broadcast(double *MyVals, int Count, int Root) const =0
int SetUseTranspose(bool UseTranspose_in)
SetUseTranpose(true) is more efficient in Amesos_Klu.
int SetUseTranspose(bool UseTranspose)
Amesos_Taucs supports only symmetric matrices, hence transpose is irrelevant, but harmless...
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
Amesos_Lapack: an interface to LAPACK.
void Set_Last_Time(double Last_Time_in)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
void Set_First_Time(double First_Time_in)
int NumericFactorization()
Performs NumericFactorization on the matrix A.
void SetRHS(Epetra_MultiVector *B)
int SetUseTranspose(bool UseTranspose)
Amesos_Superludist does not support transpose at this time.
int Solve()
Solves A X = B (or AT x = B)
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int Solve(bool Verbose=false, bool Equil=true, bool Factor=true, int perm_type=2, double pivot_thresh=-1, bool Refact=true, bool Trans=false)
All computation is performed during the call to Solve()
SpoolesOO: An object-oriented wrapper for Spooles.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
void Set_Residual(double residual_in)
int Solve(bool Factor)
All computation is performed during the call to Solve()
int Solve()
Solves A X = B (or AT x = B)
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
Superludist2_OO: An object-oriented wrapper for Superludist.
static std::ofstream log_file
int NumericFactorization()
Performs NumericFactorization on the matrix A.
static SparseSolverResult SS_Result
int Solve()
Solves A X = B (or AT X = B)
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
void SetTrans(bool trans)
Setting the transpose flag to true causes Solve() to compute A^t x = b.
Class Amesos_Umfpack: An object-oriented wrapper for UMFPACK.
int Solve()
Solves A X = B (or AT x = B)
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
void Set_Bnorm(double bnorm_in)
Amesos_Scalapack: A serial and parallel dense solver. For now, we implement only the unsymmetric ScaL...
Amesos_Taucs: An interface to the TAUCS package.
int NumericFactorization()
Performs NumericFactorization on the matrix A.