33 #include "Trilinos_Util.h" 34 #include "Epetra_LocalMap.h" 35 #include "Epetra_Map.h" 36 #include "Epetra_Vector.h" 37 #include "Epetra_MultiVector.h" 38 #include "Epetra_Export.h" 39 #include "Epetra_CrsMatrix.h" 40 #include "Epetra_LinearProblem.h" 41 #include "Epetra_Time.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 121 Epetra_Map * readMap;
122 Epetra_CrsMatrix * readA;
123 Epetra_Vector * readx;
124 Epetra_Vector * readb;
125 Epetra_Vector * readxexact;
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,
157 Epetra_CrsMatrix transposeA(
Copy, *readMap, 0);
158 Epetra_CrsMatrix *serialA ;
162 serialA = &transposeA ;
169 Epetra_Map map(readMap->NumGlobalElements(), 0, Comm);
172 if( NonContiguousMap ) {
176 int NumGlobalElements = readMap->NumGlobalElements();
177 int NumMyElements = map.NumMyElements();
178 int MyFirstElement = map.MinMyGID();
179 std::vector<int> MapMap_( NumGlobalElements );
180 readMap->MyGlobalElements( &MapMap_[0] ) ;
181 Comm.Broadcast( &MapMap_[0], NumGlobalElements, 0 ) ;
182 map_ =
new Epetra_Map( NumGlobalElements, NumMyElements, &MapMap_[MyFirstElement], 0, Comm);
184 map_ =
new Epetra_Map( map ) ;
189 Epetra_Export exporter(*readMap, *map_);
190 Epetra_CrsMatrix
A(
Copy, *map_, 0);
192 Epetra_RowMatrix * passA = 0;
193 Epetra_MultiVector * passx = 0;
194 Epetra_MultiVector * passb = 0;
195 Epetra_MultiVector * passxexact = 0;
196 Epetra_MultiVector * passresid = 0;
197 Epetra_MultiVector * passtmp = 0;
199 Epetra_MultiVector x(*map_,numsolves);
200 Epetra_MultiVector b(*map_,numsolves);
201 Epetra_MultiVector xexact(*map_,numsolves);
202 Epetra_MultiVector resid(*map_,numsolves);
203 Epetra_MultiVector tmp(*map_,numsolves);
206 Epetra_MultiVector serialx(*readMap,numsolves);
207 Epetra_MultiVector serialb(*readMap,numsolves);
208 Epetra_MultiVector serialxexact(*readMap,numsolves);
209 Epetra_MultiVector serialresid(*readMap,numsolves);
210 Epetra_MultiVector serialtmp(*readMap,numsolves);
213 if ( distribute_matrix ) {
218 A.Export(*serialA, exporter, Add);
221 assert(
A.FillComplete()==0);
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 ) ;
247 Epetra_MultiVector CopyB( *passb ) ;
249 double Anorm = passA->NormInf() ;
252 Epetra_LinearProblem Problem( (Epetra_RowMatrix *) passA,
253 (Epetra_MultiVector *) passx,
254 (Epetra_MultiVector *) passb );
256 double max_resid = 0.0;
257 for (
int j = 0 ; j < special+1 ; j++ ) {
259 Epetra_Time TotalTime( Comm ) ;
265 }
else if ( SparseSolver ==
UMFPACK ) {
266 UmfpackOO umfpack( (Epetra_RowMatrix *) passA,
267 (Epetra_MultiVector *) passx,
268 (Epetra_MultiVector *) passb ) ;
274 }
else if ( SparseSolver == SuperLU ) {
275 SuperluserialOO superluserial ;
276 superluserial.SetUserMatrix( (Epetra_RowMatrix *) passA) ;
280 superluserial.SetUseDGSSV( special == 0 ) ;
282 for (
int i= 0 ; i < numsolves ; i++ ) {
284 Epetra_Vector *passb_i = (*passb)(i) ;
285 Epetra_Vector *passx_i = (*passx)(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++ ) {
309 Epetra_Vector *passb_i = (*passb)(i) ;
310 Epetra_Vector *passx_i = (*passx)(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 ) {
332 for (
int i= 0 ; i < numsolves ; i++ ) {
334 Epetra_Vector *passb_i = (*passb)(i) ;
335 Epetra_Vector *passx_i = (*passx)(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++ ) {
360 Epetra_Vector *passb_i = (*passb)(i) ;
361 Epetra_Vector *passx_i = (*passx)(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++ ) {
386 Epetra_Vector *passb_i = (*passb)(i) ;
387 Epetra_Vector *passx_i = (*passx)(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++ ) {
414 Epetra_Vector *passb_i = (*passb)(i) ;
415 Epetra_Vector *passx_i = (*passx)(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++ ) {
438 Epetra_Vector *passb_i = (*passb)(i) ;
439 Epetra_Vector *passx_i = (*passx)(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++ ) {
471 Epetra_Vector *passb_i = (*passb)(i) ;
472 Epetra_Vector *passx_i = (*passx)(i) ;
473 Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
474 Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
479 TotalTime.ElapsedTime() );
481 if ( i < numsolves-1 )
483 TotalTime.ElapsedTime() );
486 TotalTime.ElapsedTime() );
491 #ifdef HAVE_AMESOS_LAPACK 492 }
else if ( SparseSolver ==
LAPACK ) {
499 for (
int i= 0 ; i < numsolves ; i++ ) {
501 Epetra_Vector *passb_i = (*passb)(i) ;
502 Epetra_Vector *passx_i = (*passx)(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++ ) {
529 Epetra_Vector *passb_i = (*passb)(i) ;
530 Epetra_Vector *passx_i = (*passx)(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++ ) {
557 Epetra_Vector *passb_i = (*passb)(i) ;
558 Epetra_Vector *passx_i = (*passx)(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++ ) {
585 Epetra_Vector *passb_i = (*passb)(i) ;
586 Epetra_Vector *passx_i = (*passx)(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++ ) {
613 Epetra_Vector *passb_i = (*passb)(i) ;
614 Epetra_Vector *passx_i = (*passx)(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++ ) {
641 Epetra_Vector *passb_i = (*passb)(i) ;
642 Epetra_Vector *passx_i = (*passx)(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++ ) {
670 Epetra_Vector *passb_i = (*passb)(i) ;
671 Epetra_Vector *passx_i = (*passx)(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 ) {
683 SpoolesOO spooles( (Epetra_RowMatrix *) passA,
684 (Epetra_MultiVector *) passx,
685 (Epetra_MultiVector *) passb ) ;
690 #ifdef TEST_SPOOLESSERIAL 691 }
else if ( SparseSolver == SPOOLESSERIAL ) {
692 SpoolesserialOO spoolesserial( (Epetra_RowMatrix *) passA,
693 (Epetra_MultiVector *) passx,
694 (Epetra_MultiVector *) passb ) ;
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.
Amesos_Klu: A serial, unblocked code ideal for getting started and for very sparse matrices...
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 Solve()
Solves A X = B (or AT x = B)
void Set_Total_Time(double Total_Time_in)
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)
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.
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
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 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.
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.
int CHOLMOD() error(int status, char *file, int line, char *message, cholmod_common *Common)
#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)
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...
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.
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.