30#include "Teuchos_ParameterList.hpp"
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
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 ) ;
270 umfpack.SetTrans( transpose ) ;
274 }
else if ( SparseSolver == SuperLU ) {
275 SuperluserialOO superluserial ;
276 superluserial.SetUserMatrix( (Epetra_RowMatrix *) passA) ;
279 superluserial.SetTrans( transpose ) ;
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 ) {
329 superludist2.
SetTrans( transpose ) ;
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 ) {
353 Teuchos::ParameterList ParamList ;
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 ) {
378 Teuchos::ParameterList ParamList ;
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 ) {
404 Teuchos::ParameterList ParamList ;
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 ) {
456 Teuchos::ParameterList ParamList ;
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 ) {
493 Teuchos::ParameterList ParamList ;
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 ) {
519 Teuchos::ParameterList ParamList ;
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 ) {
547 Teuchos::ParameterList ParamList ;
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 ) {
575 Teuchos::ParameterList ParamList ;
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 ) {
603 Teuchos::ParameterList ParamList ;
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 ) {
631 Teuchos::ParameterList ParamList ;
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
659 Teuchos::ParameterList ParamList ;
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 Amesos_TestMrhsSolver(Epetra_Comm &Comm, char *matrix_file, int numsolves, SparseSolverType SparseSolver, bool transpose, int special, AMESOS_MatrixType matrix_type)
int CrsMatrixTranspose(Epetra_CrsMatrix *In, Epetra_CrsMatrix *Out)
#define EPETRA_CHK_ERR(xxx)
Amesos_Dscpack: An object-oriented wrapper for Dscpack.
int Solve()
Solves A X = B (or AT x = B)
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
Amesos_Klu: A serial, unblocked code ideal for getting started and for very sparse matrices,...
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
int SetUseTranspose(bool UseTranspose_in)
SetUseTranpose(true) is more efficient in Amesos_Klu.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int Solve()
Solves A X = B (or AT x = B)
Amesos_Lapack: an interface to LAPACK.
int SetUseTranspose(bool UseTranspose_in)
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 NumericFactorization()
Performs NumericFactorization on the matrix A.
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
Amesos_Mumps: An object-oriented wrapper for the double precision version of MUMPS.
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)
int Solve()
Solves A X = B (or AT x = B)
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
Amesos_Paraklete: A serial, unblocked code ideal for getting started and for very sparse matrices,...
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int SetUseTranspose(bool UseTranspose_in)
SetUseTranpose()
int Solve()
Solves A X = B (or AT 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.
int SetParameters(Teuchos::ParameterList &ParameterList)
Set parameters from the input parameters list, returns 0 if successful.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int Solve()
Solves A X = B (or AT X = B)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int SetUseTranspose(bool UseTranspose)
SetUseTranpose()
Amesos_Scalapack: A serial and parallel dense solver. For now, we implement only the unsymmetric ScaL...
int SetUseTranspose(bool UseTranspose)
SetUseTranpose(true) is more efficient in Amesos_Scalapack.
int Solve()
Solves A X = B (or AT X = B)
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
Amesos_Superlu: Amesos interface to Xioye Li's SuperLU 3.0 serial code.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
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)
Amesos_Superludist: An object-oriented wrapper for Superludist.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int Solve()
Solves A X = B (or AT x = B)
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int SetUseTranspose(bool UseTranspose)
Amesos_Superludist does not support transpose at this time.
Amesos_Taucs: An interface to the TAUCS package.
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int SetUseTranspose(bool UseTranspose)
Amesos_Taucs supports only symmetric matrices, hence transpose is irrelevant, but harmless.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
int Solve()
Solves A X = B (or AT x = B)
Class Amesos_Umfpack: An object-oriented wrapper for UMFPACK.
int SetUseTranspose(bool UseTranspose_in)
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 SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
Epetra_SLU: An object-oriented wrapper for Xiaoye Li's serial sparse solver package: Superlu.
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()
static std::ofstream log_file
static SparseSolverResult SS_Result
void Set_Xnorm(double xnorm_in)
void Set_Last_Time(double Last_Time_in)
void Set_Error(double error_in)
void Set_Anorm(double anorm_in)
void Set_Bnorm(double bnorm_in)
void Set_Total_Time(double Total_Time_in)
void Set_First_Time(double First_Time_in)
void Set_Residual(double residual_in)
void Set_Middle_Time(double Middle_Time_in)
SpoolesOO: An object-oriented wrapper for Spooles.
void SetTrans(bool trans)
Superludist2_OO: An object-oriented wrapper for Superludist.
int Solve(bool Factor)
All computation is performed during the call to Solve()
void SetTrans(bool trans)
Setting the transpose flag to true causes Solve() to compute A^t x = b.
SuperludistOO: An object-oriented wrapper for Xiaoye Li's Superludist.
void SetTrans(bool trans)
Setting the transpose flag to true causes Solve() to compute A^t x = b.
int Solve(bool Factor)
All computation is performed during the call to Solve()