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_Export.h" 38 #include "Epetra_CrsMatrix.h" 39 #include "Epetra_LinearProblem.h" 40 #include "Epetra_Time.h" 41 #ifdef HAVE_AMESOS_SLUD 44 #ifdef HAVE_AMESOS_SLUS 47 #ifdef HAVE_AMESOS_SLUD2 50 #ifdef HAVE_AMESOS_DSCPACK 53 #ifdef HAVE_AMESOS_UMFPACK 56 #ifdef HAVE_AMESOS_KLU 59 #ifdef HAVE_AMESOS_LAPACK 62 #ifdef HAVE_AMESOS_TAUCS 65 #ifdef HAVE_AMESOS_PARDISO 68 #ifdef HAVE_AMESOS_PARAKLETE 71 #if defined(HAVE_AMESOS_MUMPS) && defined(HAVE_MPI) 74 #ifdef HAVE_AMESOS_SCALAPACK 77 #ifdef HAVE_AMESOS_SUPERLUDIST 80 #ifdef HAVE_AMESOS_SUPERLU 117 int iam = Comm.MyPID() ;
125 Epetra_Map * readMap;
126 Epetra_CrsMatrix * readA;
127 Epetra_Vector * readx;
128 Epetra_Vector * readb;
129 Epetra_Vector * readxexact;
131 std::string FileName = matrix_file ;
132 int FN_Size = FileName.size() ;
133 std::string LastFiveBytes = FileName.substr( EPETRA_MAX(0,FN_Size-5), FN_Size );
134 std::string LastFourBytes = FileName.substr( EPETRA_MAX(0,FN_Size-4), FN_Size );
135 bool NonContiguousMap =
false;
137 if ( LastFiveBytes ==
".triU" ) {
138 NonContiguousMap =
true;
140 EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( matrix_file,
false, Comm, readMap, readA, readx,
141 readb, readxexact, NonContiguousMap ) );
143 if ( LastFiveBytes ==
".triS" ) {
144 NonContiguousMap =
true;
146 EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( matrix_file,
true, Comm,
147 readMap, readA, readx,
148 readb, readxexact, NonContiguousMap ) );
150 if ( LastFourBytes ==
".mtx" ) {
151 EPETRA_CHK_ERR( Trilinos_Util_ReadMatrixMarket2Epetra( matrix_file, Comm, readMap,
152 readA, readx, readb, readxexact) );
155 Trilinos_Util_ReadHb2Epetra( matrix_file, Comm, readMap, readA, readx,
161 Epetra_CrsMatrix transposeA(Copy, *readMap, 0);
162 Epetra_CrsMatrix *serialA ;
166 serialA = &transposeA ;
172 Epetra_Map map(readMap->NumGlobalElements(), 0, Comm);
175 if( NonContiguousMap ) {
179 int NumGlobalElements = readMap->NumGlobalElements();
180 int NumMyElements = map.NumMyElements();
181 int MyFirstElement = map.MinMyGID();
182 std::vector<int> MapMap_( NumGlobalElements );
183 readMap->MyGlobalElements( &MapMap_[0] ) ;
184 Comm.Broadcast( &MapMap_[0], NumGlobalElements, 0 ) ;
185 map_ =
new Epetra_Map( NumGlobalElements, NumMyElements, &MapMap_[MyFirstElement], 0, Comm);
187 map_ =
new Epetra_Map( map ) ;
192 Epetra_Export exporter(*readMap, *map_);
193 Epetra_CrsMatrix A(Copy, *map_, 0);
195 Epetra_RowMatrix * passA = 0;
196 Epetra_MultiVector * passx = 0;
197 Epetra_MultiVector * passb = 0;
198 Epetra_MultiVector * passxexact = 0;
199 Epetra_MultiVector * passresid = 0;
200 Epetra_MultiVector * passtmp = 0;
202 Epetra_MultiVector x(*map_,numsolves);
203 Epetra_MultiVector b(*map_,numsolves);
204 Epetra_MultiVector xexact(*map_,numsolves);
205 Epetra_MultiVector resid(*map_,numsolves);
206 Epetra_MultiVector tmp(*map_,numsolves);
208 Epetra_MultiVector serialx(*readMap,numsolves);
209 Epetra_MultiVector serialb(*readMap,numsolves);
210 Epetra_MultiVector serialxexact(*readMap,numsolves);
211 Epetra_MultiVector serialresid(*readMap,numsolves);
212 Epetra_MultiVector serialtmp(*readMap,numsolves);
215 if ( distribute_matrix ) {
220 A.Export(*serialA, exporter, Add);
223 assert(A.FillComplete()==0);
229 passxexact = &xexact;
236 passxexact = &serialxexact;
237 passresid = &serialresid;
238 passtmp = &serialtmp;
241 passxexact->SetSeed(131) ;
242 passxexact->Random();
243 passx->SetSeed(11231) ;
246 passb->PutScalar( 0.0 );
247 passA->Multiply( transpose, *passxexact, *passb ) ;
249 Epetra_MultiVector CopyB( *passb ) ;
251 double Anorm = passA->NormInf() ;
254 Epetra_LinearProblem Problem( (Epetra_RowMatrix *) passA,
255 (Epetra_MultiVector *) passx,
256 (Epetra_MultiVector *) passb );
258 double max_resid = 0.0;
259 for (
int j = 0 ; j < special+1 ; j++ ) {
261 Epetra_Time TotalTime( Comm ) ;
267 }
else if ( SparseSolver ==
UMFPACK ) {
268 UmfpackOO umfpack( (Epetra_RowMatrix *) passA,
269 (Epetra_MultiVector *) passx,
270 (Epetra_MultiVector *) passb ) ;
272 umfpack.SetTrans( transpose ) ;
276 }
else if ( SparseSolver == SuperLU ) {
277 SuperluserialOO superluserial( (Epetra_RowMatrix *) passA,
278 (Epetra_MultiVector *) passx,
279 (Epetra_MultiVector *) passb ) ;
282 superluserial.SetTrans( transpose ) ;
283 superluserial.SetUseDGSSV( special == 0 ) ;
284 superluserial.Solve() ;
286 #ifdef HAVE_AMESOS_SLUD 287 }
else if ( SparseSolver == SuperLUdist ) {
292 #ifdef HAVE_AMESOS_SLUD2 293 }
else if ( SparseSolver == SuperLUdist2 ) {
295 superludist2.
SetTrans( transpose ) ;
299 }
else if ( SparseSolver == SPOOLES ) {
300 SpoolesOO spooles( (Epetra_RowMatrix *) passA,
301 (Epetra_MultiVector *) passx,
302 (Epetra_MultiVector *) passb ) ;
307 #ifdef HAVE_AMESOS_DSCPACK 308 }
else if ( SparseSolver ==
DSCPACK ) {
309 Teuchos::ParameterList ParamList ;
311 ParamList.set(
"MaxProcs", -3 );
316 #ifdef HAVE_AMESOS_UMFPACK 317 }
else if ( SparseSolver ==
UMFPACK ) {
318 Teuchos::ParameterList ParamList ;
320 ParamList.set(
"MaxProcs", -3 );
326 #ifdef HAVE_AMESOS_KLU 327 }
else if ( SparseSolver ==
KLU ) {
328 Teuchos::ParameterList ParamList ;
330 ParamList.set(
"MaxProcs", -3 );
338 #ifdef HAVE_AMESOS_PARAKLETE 339 }
else if ( SparseSolver ==
PARAKLETE ) {
340 Teuchos::ParameterList ParamList ;
342 ParamList.set(
"MaxProcs", -3 );
350 #ifdef HAVE_AMESOS_SLUS 351 }
else if ( SparseSolver == SuperLU ) {
360 #ifdef HAVE_AMESOS_LAPACK 361 }
else if ( SparseSolver ==
LAPACK ) {
362 Teuchos::ParameterList ParamList ;
363 ParamList.set(
"MaxProcs", -3 );
371 #ifdef HAVE_AMESOS_TAUCS 372 }
else if ( SparseSolver ==
TAUCS ) {
373 Teuchos::ParameterList ParamList ;
375 ParamList.set(
"MaxProcs", -3 );
383 #ifdef HAVE_AMESOS_PARDISO 384 }
else if ( SparseSolver ==
PARDISO ) {
385 Teuchos::ParameterList ParamList ;
387 ParamList.set(
"MaxProcs", -3 );
395 #ifdef HAVE_AMESOS_PARKLETE 396 }
else if ( SparseSolver == PARKLETE ) {
397 Teuchos::ParameterList ParamList ;
398 Amesos_Parklete parklete( Problem ) ;
399 ParamList.set(
"MaxProcs", -3 );
407 #if defined(HAVE_AMESOS_MUMPS) && defined(HAVE_MPI) 408 }
else if ( SparseSolver ==
MUMPS ) {
409 Teuchos::ParameterList ParamList ;
411 ParamList.set(
"MaxProcs", -3 );
419 #ifdef HAVE_AMESOS_SCALAPACK 420 }
else if ( SparseSolver ==
SCALAPACK ) {
421 Teuchos::ParameterList ParamList ;
423 ParamList.set(
"MaxProcs", -3 );
431 #ifdef HAVE_AMESOS_SUPERLUDIST 433 Teuchos::ParameterList ParamList ;
435 ParamList.set(
"MaxProcs", -3 );
444 #ifdef HAVE_AMESOS_SUPERLU 445 }
else if ( SparseSolver ==
SUPERLU ) {
446 Teuchos::ParameterList ParamList ;
448 ParamList.set(
"MaxProcs", -3 );
456 #ifdef TEST_SPOOLESSERIAL 457 }
else if ( SparseSolver == SPOOLESSERIAL ) {
458 SpoolesserialOO spoolesserial( (Epetra_RowMatrix *) passA,
459 (Epetra_MultiVector *) passx,
460 (Epetra_MultiVector *) passb ) ;
462 spoolesserial.Solve() ;
466 std::cerr <<
"\n\n#################### Requested solver not available (Or not tested with blocked RHS) on this platform #####################\n" << std::endl ;
477 std::vector <double> error(numsolves) ;
478 double max_error = 0.0;
480 passresid->Update(1.0, *passx, -1.0, *passxexact, 0.0);
482 passresid->Norm2(&error[0]);
483 for (
int i = 0 ; i< numsolves; i++ )
484 if ( error[i] > max_error ) max_error = error[i] ;
493 std::vector <double> residual(numsolves) ;
495 passtmp->PutScalar(0.0);
496 passA->Multiply( transpose, *passx, *passtmp);
497 passresid->Update(1.0, *passtmp, -1.0, *passb, 0.0);
499 passresid->Norm2(&residual[0]);
501 for (
int i = 0 ; i< numsolves; i++ )
502 if ( residual[i] > max_resid ) max_resid = residual[i] ;
507 std::vector <double> bnorm(numsolves);
508 passb->Norm2( &bnorm[0] ) ;
511 std::vector <double> xnorm(numsolves);
512 passx->Norm2( &xnorm[0] ) ;
516 if (
false && iam == 0 ) {
518 std::cout <<
" Amesos_TestMutliSolver.cpp " << std::endl ;
519 for (
int i = 0 ; i< numsolves && i < 10 ; i++ ) {
520 std::cout <<
"i=" << i
521 <<
" error = " << error[i]
522 <<
" xnorm = " << xnorm[i]
523 <<
" residual = " << residual[i]
524 <<
" bnorm = " << bnorm[i]
529 std::cout << std::endl <<
" max_resid = " << max_resid ;
530 std::cout <<
" max_error = " << max_error << std::endl ;
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)
Amesos_Dscpack: An object-oriented wrapper for Dscpack.
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.
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.
int Amesos_TestMultiSolver(Epetra_Comm &Comm, char *matrix_file, int numsolves, SparseSolverType SparseSolver, bool transpose, int special, AMESOS_MatrixType matrix_type)
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.
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
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.