Amesos Package Browser (Single Doxygen Collection) Development
Amesos_TestMrhsSolver.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //                Amesos: Direct Sparse Solver Package
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "Amesos_ConfigDefs.h"
00030 #include "Teuchos_ParameterList.hpp"
00031 //#include "Trilinos_Util_ReadTriples2Epetra.h"
00032 //#include "Trilinos_Util_ReadMatrixMarket2Epetra.h"
00033 #include "Trilinos_Util.h"
00034 #include "Epetra_LocalMap.h"
00035 #include "Epetra_Map.h"
00036 #include "Epetra_Vector.h"
00037 #include "Epetra_MultiVector.h"
00038 #include "Epetra_Export.h"
00039 #include "Epetra_CrsMatrix.h"
00040 #include "Epetra_LinearProblem.h"
00041 #include "Epetra_Time.h"
00042 #ifdef HAVE_AMESOS_DSCPACK
00043 #include "Amesos_Dscpack.h"
00044 #endif
00045 #ifdef HAVE_AMESOS_LAPACK
00046 #include "Amesos_Lapack.h"
00047 #endif
00048 #ifdef HAVE_AMESOS_SLUS
00049 #include "Epetra_SLU.h"
00050 #endif
00051 #ifdef HAVE_AMESOS_UMFPACK
00052 #include "Amesos_Umfpack.h"
00053 #endif
00054 #ifdef HAVE_AMESOS_KLU
00055 #include "Amesos_Klu.h"
00056 #endif
00057 #ifdef HAVE_AMESOS_TAUCS
00058 #include "Amesos_Taucs.h"
00059 #endif
00060 #ifdef HAVE_AMESOS_PARDISO
00061 #include "Amesos_Pardiso.h"
00062 #endif
00063 #ifdef HAVE_AMESOS_PARAKLETE
00064 #include "Amesos_Paraklete.h"
00065 #endif
00066 #ifdef HAVE_AMESOS_MUMPS
00067 #include "Amesos_Mumps.h"
00068 #endif
00069 #ifdef HAVE_AMESOS_SCALAPACK
00070 #include "Amesos_Scalapack.h"
00071 #endif
00072 #ifdef HAVE_AMESOS_SLUD
00073 #include "SuperludistOO.h"
00074 #endif
00075 #ifdef HAVE_AMESOS_SUPERLU
00076 #include "Amesos_Superlu.h"
00077 #endif
00078 #ifdef HAVE_AMESOS_SUPERLUDIST
00079 #include "Amesos_Superludist.h"
00080 #endif
00081 #ifdef HAVE_AMESOS_SLUD2
00082 #include "Superludist2_OO.h"
00083 #endif
00084 #ifdef TEST_SPOOLES
00085 #include "SpoolesOO.h"
00086 #endif
00087 #include "SparseSolverResult.h"
00088 #include "Amesos_TestSolver.h"
00089 #include "CrsMatrixTranspose.h"
00090 #include "SparseDirectTimingVars.h"
00091 
00092 #include <vector>
00093 //
00094 //  TestMrhsSolver.cpp reads in a matrix in Harwell-Boeing format, 
00095 //  calls one of the sparse direct solvers, using multiple right hand sides
00096 //  (one per solve) and computes the error and residual.  
00097 //
00098 //  TestSolver ignores the Harwell-Boeing right hand sides, creating
00099 //  random right hand sides instead.  
00100 //
00101 //  TestMrhsSolver can test either A x = b or A^T x = b.
00102 //  This can be a bit confusing because sparse direct solvers 
00103 //  use compressed column storage - the transpose of Trilinos'
00104 //  sparse row storage.
00105 //
00106 //  Matrices:
00107 //    readA - Serial.  As read from the file.
00108 //    transposeA - Serial.  The transpose of readA.
00109 //    serialA - if (transpose) then transposeA else readA 
00110 //    distributedA - readA distributed to all processes
00111 //    passA - if ( distributed ) then distributedA else serialA
00112 //
00113 //
00114 int Amesos_TestMrhsSolver( Epetra_Comm &Comm, char *matrix_file, int numsolves, 
00115          SparseSolverType SparseSolver, bool transpose, 
00116          int special, AMESOS_MatrixType matrix_type ) {
00117 
00118 
00119   Comm.Barrier();
00120 
00121   Epetra_Map * readMap;
00122   Epetra_CrsMatrix * readA; 
00123   Epetra_Vector * readx; 
00124   Epetra_Vector * readb;
00125   Epetra_Vector * readxexact;
00126 
00127   std::string FileName = matrix_file ;
00128   int FN_Size = FileName.size() ; 
00129   std::string LastFiveBytes = FileName.substr( EPETRA_MAX(0,FN_Size-5), FN_Size );
00130   std::string LastFourBytes = FileName.substr( EPETRA_MAX(0,FN_Size-4), FN_Size );
00131   bool NonContiguousMap = false; 
00132 
00133   if ( LastFiveBytes == ".triU" ) { 
00134     // Call routine to read in unsymmetric Triplet matrix
00135     NonContiguousMap = true; 
00136     EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( matrix_file, false, Comm, readMap, readA, readx, 
00137                   readb, readxexact, NonContiguousMap ) );
00138   } else {
00139     if ( LastFiveBytes == ".triS" ) { 
00140       NonContiguousMap = true; 
00141       // Call routine to read in symmetric Triplet matrix
00142       EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( matrix_file, true, Comm, readMap, readA, readx, 
00143               readb, readxexact) );
00144     } else {
00145       if (  LastFourBytes == ".mtx" ) { 
00146   EPETRA_CHK_ERR( Trilinos_Util_ReadMatrixMarket2Epetra( matrix_file, Comm, readMap, 
00147                      readA, readx, readb, readxexact) );
00148       } else {
00149   // Call routine to read in HB problem
00150   Trilinos_Util_ReadHb2Epetra( matrix_file, Comm, readMap, readA, readx, 
00151                  readb, readxexact) ;
00152       }
00153     }
00154   }
00155 
00156 
00157   Epetra_CrsMatrix transposeA(Copy, *readMap, 0);
00158   Epetra_CrsMatrix *serialA ; 
00159 
00160   if ( transpose ) {
00161     assert( CrsMatrixTranspose( readA, &transposeA ) == 0 ); 
00162     serialA = &transposeA ; 
00163   } else {
00164     serialA = readA ; 
00165   }
00166 
00167   
00168   // Create uniform distributed map
00169   Epetra_Map map(readMap->NumGlobalElements(), 0, Comm);
00170   Epetra_Map* map_;
00171 
00172   if( NonContiguousMap ) {
00173     //
00174     //  map gives us NumMyElements and MyFirstElement;
00175     //
00176     int NumGlobalElements =  readMap->NumGlobalElements();
00177     int NumMyElements = map.NumMyElements();
00178     int MyFirstElement = map.MinMyGID();
00179     std::vector<int> MapMap_( NumGlobalElements );
00180     readMap->MyGlobalElements( &MapMap_[0] ) ;
00181     Comm.Broadcast( &MapMap_[0], NumGlobalElements, 0 ) ; 
00182     map_ = new Epetra_Map( NumGlobalElements, NumMyElements, &MapMap_[MyFirstElement], 0, Comm);
00183   } else {
00184     map_ = new Epetra_Map( map ) ; 
00185   }
00186 
00187 
00188   // Create Exporter to distribute read-in matrix and vectors
00189   Epetra_Export exporter(*readMap, *map_);
00190   Epetra_CrsMatrix A(Copy, *map_, 0);
00191 
00192   Epetra_RowMatrix * passA = 0; 
00193   Epetra_MultiVector * passx = 0; 
00194   Epetra_MultiVector * passb = 0;
00195   Epetra_MultiVector * passxexact = 0;
00196   Epetra_MultiVector * passresid = 0;
00197   Epetra_MultiVector * passtmp = 0;
00198 
00199   Epetra_MultiVector x(*map_,numsolves);
00200   Epetra_MultiVector b(*map_,numsolves);
00201   Epetra_MultiVector xexact(*map_,numsolves);
00202   Epetra_MultiVector resid(*map_,numsolves);
00203   Epetra_MultiVector tmp(*map_,numsolves);
00204 
00205 
00206   Epetra_MultiVector serialx(*readMap,numsolves);
00207   Epetra_MultiVector serialb(*readMap,numsolves);
00208   Epetra_MultiVector serialxexact(*readMap,numsolves);
00209   Epetra_MultiVector serialresid(*readMap,numsolves);
00210   Epetra_MultiVector serialtmp(*readMap,numsolves);
00211 
00212   bool distribute_matrix = ( matrix_type == AMESOS_Distributed ) ; 
00213   if ( distribute_matrix ) { 
00214     //
00215     //  Initialize x, b and xexact to the values read in from the file
00216     //
00217 
00218     A.Export(*serialA, exporter, Add);
00219     Comm.Barrier();
00220 
00221     assert(A.FillComplete()==0);    
00222     Comm.Barrier();
00223 
00224     passA = &A; 
00225     passx = &x; 
00226     passb = &b;
00227     passxexact = &xexact;
00228     passresid = &resid;
00229     passtmp = &tmp;
00230   } else { 
00231     passA = serialA; 
00232     passx = &serialx; 
00233     passb = &serialb;
00234     passxexact = &serialxexact;
00235     passresid = &serialresid;
00236     passtmp = &serialtmp;
00237   }
00238 
00239   passxexact->SetSeed(131) ; 
00240   passxexact->Random();
00241   passx->SetSeed(11231) ; 
00242   passx->Random();
00243 
00244   passb->PutScalar( 0.0 );
00245   passA->Multiply( transpose, *passxexact, *passb ) ; 
00246 
00247   Epetra_MultiVector CopyB( *passb ) ;
00248 
00249   double Anorm = passA->NormInf() ; 
00250   SparseDirectTimingVars::SS_Result.Set_Anorm(Anorm) ;
00251 
00252   Epetra_LinearProblem Problem(  (Epetra_RowMatrix *) passA, 
00253          (Epetra_MultiVector *) passx, 
00254          (Epetra_MultiVector *) passb );
00255 
00256   double max_resid = 0.0;
00257   for ( int j = 0 ; j < special+1 ; j++ ) { 
00258     
00259     Epetra_Time TotalTime( Comm ) ; 
00260     if ( false ) { 
00261 #ifdef TEST_UMFPACK
00262 
00263       unused code
00264 
00265     } else if ( SparseSolver == UMFPACK ) { 
00266       UmfpackOO umfpack( (Epetra_RowMatrix *) passA, 
00267        (Epetra_MultiVector *) passx, 
00268        (Epetra_MultiVector *) passb ) ; 
00269       
00270       umfpack.SetTrans( transpose ) ; 
00271       umfpack.Solve() ; 
00272 #endif
00273 #ifdef TEST_SUPERLU
00274     } else if ( SparseSolver == SuperLU ) { 
00275       SuperluserialOO superluserial ; 
00276       superluserial.SetUserMatrix( (Epetra_RowMatrix *) passA) ; 
00277 
00278       superluserial.SetPermc( SuperLU_permc ) ; 
00279       superluserial.SetTrans( transpose ) ; 
00280       superluserial.SetUseDGSSV( special == 0 ) ; 
00281 
00282       for ( int i= 0 ; i < numsolves ; i++ ) { 
00283   //    set up to sovle A X[:,i] = B[:,i]
00284   Epetra_Vector *passb_i = (*passb)(i) ;
00285   Epetra_Vector *passx_i = (*passx)(i) ;
00286   superluserial.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
00287   superluserial.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
00288   //      superluserial.SetRHS( (Epetra_MultiVector *) passb_i ; 
00289   superluserial.Solve() ; 
00290   if ( i == 0 ) {
00291     SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() ); 
00292   } else { 
00293     if ( i < numsolves-1 ) 
00294       SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() ); 
00295     else
00296       SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() ); 
00297   }
00298 
00299       }
00300 #endif
00301 #ifdef HAVE_AMESOS_SLUD
00302     } else if ( SparseSolver == SuperLUdist ) { 
00303       SuperludistOO superludist( Problem ) ; 
00304       superludist.SetTrans( transpose ) ; 
00305 
00306       bool factor = true; 
00307       for ( int i= 0 ; i < numsolves ; i++ ) { 
00308   //    set up to sovle A X[:,i] = B[:,i]
00309   Epetra_Vector *passb_i = (*passb)(i) ;
00310   Epetra_Vector *passx_i = (*passx)(i) ;
00311   Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
00312   Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
00313   EPETRA_CHK_ERR( superludist.Solve( factor ) ); 
00314   factor = false; 
00315   if ( i == 0 ) 
00316     SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() ); 
00317   else { 
00318     if ( i < numsolves-1 ) 
00319       SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() ); 
00320     else
00321       SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() ); 
00322   }
00323 
00324       }
00325 #endif
00326 #ifdef HAVE_AMESOS_SLUD2
00327     } else if ( SparseSolver == SuperLUdist2 ) { 
00328       Superludist2_OO superludist2( Problem ) ; 
00329       superludist2.SetTrans( transpose ) ; 
00330 
00331       bool factor = true; 
00332       for ( int i= 0 ; i < numsolves ; i++ ) { 
00333   //    set up to sovle A X[:,i] = B[:,i]
00334   Epetra_Vector *passb_i = (*passb)(i) ;
00335   Epetra_Vector *passx_i = (*passx)(i) ;
00336   Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
00337   Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
00338   EPETRA_CHK_ERR( superludist2.Solve( factor ) ); 
00339   factor = false; 
00340   if ( i == 0 ) 
00341     SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() ); 
00342   else { 
00343     if ( i < numsolves-1 ) 
00344       SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() ); 
00345     else
00346       SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() ); 
00347   }
00348 
00349       }
00350 #endif
00351 #ifdef HAVE_AMESOS_DSCPACK
00352     } else if ( SparseSolver == DSCPACK ) { 
00353       Teuchos::ParameterList ParamList ;
00354       Amesos_Dscpack dscpack( Problem ) ; 
00355       ParamList.set( "MaxProcs", -3 );
00356       EPETRA_CHK_ERR( dscpack.SetParameters( ParamList ) ); 
00357 
00358       bool factor = true; 
00359       for ( int i= 0 ; i < numsolves ; i++ ) { 
00360   //    set up to sovle A X[:,i] = B[:,i]
00361   Epetra_Vector *passb_i = (*passb)(i) ;
00362   Epetra_Vector *passx_i = (*passx)(i) ;
00363   Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
00364   Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
00365   EPETRA_CHK_ERR( dscpack.Solve( ) ); 
00366   factor = false; 
00367   if ( i == 0 ) 
00368     SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() ); 
00369   else { 
00370     if ( i < numsolves-1 ) 
00371       SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() ); 
00372     else
00373       SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() ); 
00374   }
00375 
00376       }
00377 #endif
00378 #ifdef HAVE_AMESOS_UMFPACK
00379     } else if ( SparseSolver == UMFPACK ) { 
00380       Teuchos::ParameterList ParamList ;
00381       Amesos_Umfpack umfpack( Problem ) ; 
00382       ParamList.set( "MaxProcs", -3 );
00383       EPETRA_CHK_ERR( umfpack.SetParameters( ParamList ) ); 
00384       EPETRA_CHK_ERR( umfpack.SetUseTranspose( transpose ) ); 
00385 
00386       bool factor = true; 
00387       for ( int i= 0 ; i < numsolves ; i++ ) { 
00388   //    set up to sovle A X[:,i] = B[:,i]
00389   Epetra_Vector *passb_i = (*passb)(i) ;
00390   Epetra_Vector *passx_i = (*passx)(i) ;
00391   Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
00392   Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
00393   EPETRA_CHK_ERR( umfpack.Solve( ) ); 
00394   factor = false; 
00395   if ( i == 0 ) 
00396     SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() ); 
00397   else { 
00398     if ( i < numsolves-1 ) 
00399       SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() ); 
00400     else
00401       SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() ); 
00402   }
00403 
00404       }
00405 #endif
00406 #ifdef HAVE_AMESOS_SUPERLU
00407     } else if ( SparseSolver == SUPERLU ) { 
00408       Teuchos::ParameterList ParamList ;
00409       Amesos_Superlu superlu( Problem ) ; 
00410       ParamList.set( "MaxProcs", -3 );
00411       EPETRA_CHK_ERR( superlu.SetParameters( ParamList ) ); 
00412       EPETRA_CHK_ERR( superlu.SetUseTranspose( transpose ) ); 
00413 
00414       bool factor = true; 
00415       EPETRA_CHK_ERR( superlu.SymbolicFactorization(  ) ); 
00416       EPETRA_CHK_ERR( superlu.NumericFactorization(  ) ); 
00417       for ( int i= 0 ; i < numsolves ; i++ ) { 
00418   //    set up to sovle A X[:,i] = B[:,i]
00419   Epetra_Vector *passb_i = (*passb)(i) ;
00420   Epetra_Vector *passx_i = (*passx)(i) ;
00421   Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
00422   Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
00423   EPETRA_CHK_ERR( superlu.Solve( ) ); 
00424   factor = false; 
00425   if ( i == 0 ) 
00426     SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() ); 
00427   else { 
00428     if ( i < numsolves-1 ) 
00429       SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() ); 
00430     else
00431       SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() ); 
00432   }
00433 
00434       }
00435 #endif
00436 #ifdef HAVE_AMESOS_SLUS
00437     } else if ( SparseSolver == SuperLU ) { 
00438       Epetra_SLU superluserial( &Problem ) ;
00439       
00440       bool factor = true; 
00441 
00442       for ( int i= 0 ; i < numsolves ; i++ ) { 
00443   //    set up to sovle A X[:,i] = B[:,i]
00444   Epetra_Vector *passb_i = (*passb)(i) ;
00445   Epetra_Vector *passx_i = (*passx)(i) ;
00446   Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
00447   Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
00448   EPETRA_CHK_ERR( superluserial.Solve( true, false, factor, 2, -1, true, transpose ) ); 
00449   //  factor = false; 
00450   if ( i == 0 ) 
00451     SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() ); 
00452   else { 
00453     if ( i < numsolves-1 ) 
00454       SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() ); 
00455     else
00456       SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() ); 
00457   }
00458 
00459       }
00460 #endif
00461 #ifdef HAVE_AMESOS_KLU
00462     } else if ( SparseSolver == KLU ) { 
00463       Teuchos::ParameterList ParamList ;
00464       //      ParamList.set("OutputLevel",2);
00465       Amesos_Klu klu( Problem ) ; 
00466       // ParamList.set ("ScaleMethod", 0) ;
00467       ParamList.set( "MaxProcs", -3 );
00468       EPETRA_CHK_ERR( klu.SetParameters( ParamList ) ); 
00469       ParamList.set( "MaxProcs", -3 );
00470       EPETRA_CHK_ERR( klu.SetParameters( ParamList ) ); 
00471       EPETRA_CHK_ERR( klu.SetUseTranspose( transpose ) ); 
00472 
00473       bool factor = true; 
00474       EPETRA_CHK_ERR( klu.SymbolicFactorization(  ) ); 
00475       for ( int trials = 0 ; trials <= 1 ; trials++) {
00476     EPETRA_CHK_ERR( klu.NumericFactorization(  ) ); 
00477     for ( int i= 0 ; i < numsolves ; i++ ) {
00478       //    set up to sovle A X[:,i] = B[:,i]
00479       Epetra_Vector *passb_i = (*passb)(i) ;
00480       Epetra_Vector *passx_i = (*passx)(i) ;
00481       Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
00482       Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
00483 
00484       EPETRA_CHK_ERR( klu.Solve( ) ); 
00485       factor = false; 
00486       if ( i == 0 ) {
00487         SparseDirectTimingVars::SS_Result.Set_First_Time(
00488           TotalTime.ElapsedTime() ); 
00489       } else {
00490         if ( i < numsolves-1 ) 
00491     SparseDirectTimingVars::SS_Result.Set_Middle_Time(
00492       TotalTime.ElapsedTime() ); 
00493         else
00494     SparseDirectTimingVars::SS_Result.Set_Last_Time(
00495       TotalTime.ElapsedTime() ); 
00496       }
00497     }
00498       }
00499 
00500 #endif
00501 #ifdef HAVE_AMESOS_LAPACK
00502     } else if ( SparseSolver == LAPACK ) { 
00503       Teuchos::ParameterList ParamList ;
00504       Amesos_Lapack lapack( Problem ) ; 
00505       EPETRA_CHK_ERR( lapack.SetUseTranspose( transpose ) ); 
00506 
00507       bool factor = true; 
00508       EPETRA_CHK_ERR( lapack.SymbolicFactorization(  ) ); 
00509       EPETRA_CHK_ERR( lapack.NumericFactorization(  ) ); 
00510       for ( int i= 0 ; i < numsolves ; i++ ) { 
00511   //    set up to sovle A X[:,i] = B[:,i]
00512   Epetra_Vector *passb_i = (*passb)(i) ;
00513   Epetra_Vector *passx_i = (*passx)(i) ;
00514   Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
00515   Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
00516   EPETRA_CHK_ERR( lapack.Solve( ) ); 
00517   factor = false; 
00518   if ( i == 0 ) 
00519     SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() ); 
00520   else { 
00521     if ( i < numsolves-1 ) 
00522       SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() ); 
00523     else
00524       SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() ); 
00525   }
00526 
00527       }
00528 #endif
00529 #ifdef HAVE_AMESOS_TAUCS
00530     } else if ( SparseSolver == TAUCS ) { 
00531       Teuchos::ParameterList ParamList ;
00532       Amesos_Taucs taucs( Problem ) ; 
00533       ParamList.set( "MaxProcs", -3 );
00534       EPETRA_CHK_ERR( taucs.SetParameters( ParamList ) ); 
00535       EPETRA_CHK_ERR( taucs.SetUseTranspose( transpose ) ); 
00536       EPETRA_CHK_ERR( taucs.SymbolicFactorization( ) ); 
00537       EPETRA_CHK_ERR( taucs.NumericFactorization( ) ); 
00538 
00539       for ( int i= 0 ; i < numsolves ; i++ ) { 
00540   //    set up to sovle A X[:,i] = B[:,i]
00541   Epetra_Vector *passb_i = (*passb)(i) ;
00542   Epetra_Vector *passx_i = (*passx)(i) ;
00543   Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
00544   Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
00545   EPETRA_CHK_ERR( taucs.Solve( ) ); 
00546   //  factor = false; 
00547   if ( i == 0 ) 
00548     SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() ); 
00549   else { 
00550     if ( i < numsolves-1 ) 
00551       SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() ); 
00552     else
00553       SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() ); 
00554   }
00555 
00556       }
00557 #endif
00558 #ifdef HAVE_AMESOS_PARDISO
00559     } else if ( SparseSolver == PARDISO ) { 
00560       Teuchos::ParameterList ParamList ;
00561       Amesos_Pardiso pardiso( Problem ) ; 
00562       ParamList.set( "MaxProcs", -3 );
00563       EPETRA_CHK_ERR( pardiso.SetParameters( ParamList ) ); 
00564       EPETRA_CHK_ERR( pardiso.SetUseTranspose( transpose ) ); 
00565       EPETRA_CHK_ERR( pardiso.SymbolicFactorization( ) ); 
00566       EPETRA_CHK_ERR( pardiso.NumericFactorization( ) ); 
00567 
00568       for ( int i= 0 ; i < numsolves ; i++ ) { 
00569   //    set up to sovle A X[:,i] = B[:,i]
00570   Epetra_Vector *passb_i = (*passb)(i) ;
00571   Epetra_Vector *passx_i = (*passx)(i) ;
00572   Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
00573   Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
00574   EPETRA_CHK_ERR( pardiso.Solve( ) ); 
00575   //  factor = false; 
00576   if ( i == 0 ) 
00577     SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() ); 
00578   else { 
00579     if ( i < numsolves-1 ) 
00580       SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() ); 
00581     else
00582       SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() ); 
00583   }
00584 
00585       }
00586 #endif
00587 #ifdef HAVE_AMESOS_PARAKLETE
00588     } else if ( SparseSolver == PARAKLETE ) { 
00589       Teuchos::ParameterList ParamList ;
00590       Amesos_Paraklete paraklete( Problem ) ; 
00591       ParamList.set( "MaxProcs", -3 );
00592       EPETRA_CHK_ERR( paraklete.SetParameters( ParamList ) ); 
00593       EPETRA_CHK_ERR( paraklete.SetUseTranspose( transpose ) ); 
00594       EPETRA_CHK_ERR( paraklete.SymbolicFactorization( ) ); 
00595       EPETRA_CHK_ERR( paraklete.NumericFactorization( ) ); 
00596 
00597       for ( int i= 0 ; i < numsolves ; i++ ) { 
00598   //    set up to sovle A X[:,i] = B[:,i]
00599   Epetra_Vector *passb_i = (*passb)(i) ;
00600   Epetra_Vector *passx_i = (*passx)(i) ;
00601   Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
00602   Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
00603   EPETRA_CHK_ERR( paraklete.Solve( ) ); 
00604   //  factor = false; 
00605   if ( i == 0 ) 
00606     SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() ); 
00607   else { 
00608     if ( i < numsolves-1 ) 
00609       SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() ); 
00610     else
00611       SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() ); 
00612   }
00613 
00614       }
00615 #endif
00616 #ifdef HAVE_AMESOS_MUMPS
00617     } else if ( SparseSolver == MUMPS ) { 
00618       Teuchos::ParameterList ParamList ;
00619       Amesos_Mumps mumps( Problem ) ; 
00620       ParamList.set( "MaxProcs", -3 );
00621       EPETRA_CHK_ERR( mumps.SetParameters( ParamList ) ); 
00622       EPETRA_CHK_ERR( mumps.SetUseTranspose( transpose ) ); 
00623       EPETRA_CHK_ERR( mumps.SymbolicFactorization( ) ); 
00624       EPETRA_CHK_ERR( mumps.NumericFactorization( ) ); 
00625 
00626       for ( int i= 0 ; i < numsolves ; i++ ) { 
00627   //    set up to sovle A X[:,i] = B[:,i]
00628   Epetra_Vector *passb_i = (*passb)(i) ;
00629   Epetra_Vector *passx_i = (*passx)(i) ;
00630   Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
00631   Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
00632   EPETRA_CHK_ERR( mumps.Solve( ) ); 
00633   //  factor = false; 
00634   if ( i == 0 ) 
00635     SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() ); 
00636   else { 
00637     if ( i < numsolves-1 ) 
00638       SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() ); 
00639     else
00640       SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() ); 
00641   }
00642 
00643       }
00644 #endif
00645 #ifdef HAVE_AMESOS_SCALAPACK
00646     } else if ( SparseSolver == SCALAPACK ) { 
00647       Teuchos::ParameterList ParamList ;
00648       Amesos_Scalapack scalapack( Problem ) ; 
00649       ParamList.set( "MaxProcs", -3 );
00650       EPETRA_CHK_ERR( scalapack.SetParameters( ParamList ) ); 
00651       EPETRA_CHK_ERR( scalapack.SetUseTranspose( transpose ) ); 
00652 
00653       EPETRA_CHK_ERR( scalapack.SymbolicFactorization( ) ); 
00654       EPETRA_CHK_ERR( scalapack.NumericFactorization( ) ); 
00655       bool factor = true; 
00656       for ( int i= 0 ; i < numsolves ; i++ ) { 
00657   //    set up to sovle A X[:,i] = B[:,i]
00658   Epetra_Vector *passb_i = (*passb)(i) ;
00659   Epetra_Vector *passx_i = (*passx)(i) ;
00660   Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
00661   Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
00662   EPETRA_CHK_ERR( scalapack.Solve( ) ); 
00663   //  factor = false; 
00664   if ( i == 0 ) 
00665     SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() ); 
00666   else { 
00667     if ( i < numsolves-1 ) 
00668       SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() ); 
00669     else
00670       SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() ); 
00671   }
00672 
00673       }
00674 #endif
00675 #ifdef HAVE_AMESOS_SUPERLUDIST
00676     } else if ( SparseSolver == SUPERLUDIST ) { 
00677       Teuchos::ParameterList ParamList ;
00678       ParamList.set( "MaxProcs", -3 );
00679       Amesos_Superludist superludist( Problem ) ; 
00680       EPETRA_CHK_ERR( superludist.SetParameters( ParamList ) ); 
00681       EPETRA_CHK_ERR( superludist.SetUseTranspose( transpose ) ); 
00682       EPETRA_CHK_ERR( superludist.SymbolicFactorization(  ) ); 
00683       EPETRA_CHK_ERR( superludist.NumericFactorization(  ) ); 
00684       SparseDirectTimingVars::SS_Result.Set_First_Time( TotalTime.ElapsedTime() ); 
00685 
00686       bool factor = true; 
00687       for ( int i= 0 ; i < numsolves ; i++ ) { 
00688   //    set up to sovle A X[:,i] = B[:,i]
00689   Epetra_Vector *passb_i = (*passb)(i) ;
00690   Epetra_Vector *passx_i = (*passx)(i) ;
00691   Problem.SetLHS( dynamic_cast<Epetra_MultiVector *>(passx_i) ) ;
00692   Problem.SetRHS( dynamic_cast<Epetra_MultiVector *>(passb_i) );
00693   EPETRA_CHK_ERR( superludist.Solve( ) ); 
00694   factor = false; 
00695   if ( i < numsolves-1 ) 
00696     SparseDirectTimingVars::SS_Result.Set_Middle_Time( TotalTime.ElapsedTime() ); 
00697   else
00698     SparseDirectTimingVars::SS_Result.Set_Last_Time( TotalTime.ElapsedTime() ); 
00699       }
00700 #endif
00701 #ifdef TEST_SPOOLES
00702     } else if ( SparseSolver == SPOOLES ) { 
00703       SpoolesOO spooles( (Epetra_RowMatrix *) passA, 
00704        (Epetra_MultiVector *) passx, 
00705        (Epetra_MultiVector *) passb ) ; 
00706     
00707       spooles.SetTrans( transpose ) ; 
00708       spooles.Solve() ;
00709 #endif 
00710 #ifdef TEST_SPOOLESSERIAL
00711     } else if ( SparseSolver == SPOOLESSERIAL ) { 
00712       SpoolesserialOO spoolesserial( (Epetra_RowMatrix *) passA, 
00713              (Epetra_MultiVector *) passx, 
00714              (Epetra_MultiVector *) passb ) ; 
00715     
00716       spoolesserial.Solve() ;
00717 #endif 
00718     } else { 
00719       SparseDirectTimingVars::log_file << "Solver not implemented yet" << std::endl ;
00720       std::cerr << "\n\n####################  Requested solver not available (Or not tested with multiple RHS) on this platform #####################\n" << std::endl ;
00721     }
00722 
00723     SparseDirectTimingVars::SS_Result.Set_Total_Time( TotalTime.ElapsedTime() ); 
00724 
00725     //
00726     //  Compute the error = norm(xcomp - xexact )
00727     //
00728     std::vector <double> error(numsolves) ; 
00729     double max_error = 0.0;
00730   
00731     passresid->Update(1.0, *passx, -1.0, *passxexact, 0.0);
00732 
00733     passresid->Norm2(&error[0]);
00734     for ( int i = 0 ; i< numsolves; i++ ) 
00735       if ( error[i] > max_error ) max_error = error[i] ; 
00736     SparseDirectTimingVars::SS_Result.Set_Error(max_error) ;
00737 
00738     //  passxexact->Norm2(&error[0] ) ; 
00739     //  passx->Norm2(&error ) ; 
00740 
00741     //
00742     //  Compute the residual = norm(Ax - b)
00743     //
00744     std::vector <double> residual(numsolves) ; 
00745   
00746     passtmp->PutScalar(0.0);
00747     passA->Multiply( transpose, *passx, *passtmp);
00748     passresid->Update(1.0, *passtmp, -1.0, *passb, 0.0); 
00749     //    passresid->Update(1.0, *passtmp, -1.0, CopyB, 0.0); 
00750     passresid->Norm2(&residual[0]);
00751 
00752     for ( int i = 0 ; i< numsolves; i++ ) 
00753       if ( residual[i] > max_resid ) max_resid = residual[i] ; 
00754 
00755 
00756     SparseDirectTimingVars::SS_Result.Set_Residual(max_resid) ;
00757     
00758     std::vector <double> bnorm(numsolves); 
00759     passb->Norm2( &bnorm[0] ) ; 
00760     SparseDirectTimingVars::SS_Result.Set_Bnorm(bnorm[0]) ;
00761 
00762     std::vector <double> xnorm(numsolves); 
00763     passx->Norm2( &xnorm[0] ) ; 
00764     SparseDirectTimingVars::SS_Result.Set_Xnorm(xnorm[0]) ;
00765 
00766   }
00767   delete readA;
00768   delete readx;
00769   delete readb;
00770   delete readxexact;
00771   delete readMap;
00772   delete map_;
00773 
00774   Comm.Barrier();
00775    return 0;
00776 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines