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 ) {
304 superludist.SetTrans( transpose ) ;
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) );
442 EPETRA_CHK_ERR( superluserial.Solve(
true,
false, factor, 2, -1,
true, transpose ) );
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 ) ;
687 spooles.SetTrans( transpose ) ;
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] ) ;