33 #include "Epetra_MpiComm.h" 35 #include "Epetra_SerialComm.h" 37 #include "Epetra_Vector.h" 38 #include "Epetra_Time.h" 39 #include "Epetra_RowMatrix.h" 40 #include "Epetra_CrsMatrix.h" 43 #include "Teuchos_ParameterList.hpp" 44 #include "Galeri_Maps.h" 45 #include "Galeri_CrsMatrices.h" 50 #include "Trilinos_Util.h" 51 #include "Trilinos_Util_ReadMatrixMarket2Epetra.h" 53 #include "Teuchos_RCP.hpp" 54 #include "Epetra_Export.h" 57 Epetra_Map *& readMap,
58 const bool transpose,
const bool distribute,
59 bool& symmetric, Epetra_CrsMatrix *& Matrix ) {
61 Epetra_CrsMatrix * readA = 0;
62 Epetra_Vector * readx = 0;
63 Epetra_Vector * readb = 0;
64 Epetra_Vector * readxexact = 0;
70 FILE *in_file = fopen( in_filename,
"r");
74 filename = &in_filename[1] ;
77 filename = in_filename ;
82 std::string FileName (filename);
84 int FN_Size = FileName.size() ;
85 std::string LastFiveBytes = FileName.substr( EPETRA_MAX(0,FN_Size-5), FN_Size );
86 std::string LastFourBytes = FileName.substr( EPETRA_MAX(0,FN_Size-4), FN_Size );
88 if ( LastFiveBytes ==
".TimD" ) {
90 EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename,
false, Comm, readMap, readA, readx,
91 readb, readxexact,
false,
true,
true ) );
94 if ( LastFiveBytes ==
".triU" ) {
96 EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename,
false, Comm, readMap, readA, readx,
100 if ( LastFiveBytes ==
".triS" ) {
102 EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename,
true, Comm, readMap, readA, readx,
103 readb, readxexact) );
106 if ( LastFourBytes ==
".mtx" ) {
107 EPETRA_CHK_ERR( Trilinos_Util_ReadMatrixMarket2Epetra( filename, Comm, readMap,
108 readA, readx, readb, readxexact) );
109 FILE* in_file = fopen( filename,
"r");
110 assert (in_file != NULL) ;
111 const int BUFSIZE = 800 ;
112 char buffer[BUFSIZE] ;
113 fgets( buffer, BUFSIZE, in_file ) ;
114 std::string headerline1 = buffer;
116 if ( headerline1.find(
"symmetric") < BUFSIZE ) symmetric =
true;
118 if ( headerline1.find(
"symmetric") != std::string::npos) symmetric =
true;
125 Trilinos_Util_ReadHb2Epetra( filename, Comm, readMap, readA, readx,
127 if ( LastFourBytes ==
".rsa" ) symmetric = true ;
134 if ( readb )
delete readb;
135 if ( readx )
delete readx;
136 if ( readxexact )
delete readxexact;
138 Epetra_CrsMatrix *serialA ;
139 Epetra_CrsMatrix *transposeA;
142 transposeA =
new Epetra_CrsMatrix( Copy, *readMap, 0 );
144 serialA = transposeA ;
151 assert( (
void *) &serialA->Graph() ) ;
152 assert( (
void *) &serialA->RowMap() ) ;
153 assert( serialA->RowMap().SameAs(*readMap) ) ;
157 Epetra_Map DistMap(readMap->NumGlobalElements(), 0, Comm);
160 Epetra_Export exporter( *readMap, DistMap );
162 Epetra_CrsMatrix *Amat =
new Epetra_CrsMatrix( Copy, DistMap, 0 );
163 Amat->Export(*serialA, exporter, Add);
164 assert(Amat->FillComplete()==0);
202 int main(
int argc,
char *argv[])
205 MPI_Init(&argc, &argv);
206 Epetra_MpiComm Comm(MPI_COMM_WORLD);
208 Epetra_SerialComm Comm;
211 bool verbose = (Comm.MyPID() == 0);
212 double TotalResidual = 0.0;
218 #ifndef FILENAME_SPECIFIED_ON_COMMAND_LINE 219 ParameterList GaleriList;
220 GaleriList.set(
"nx", 4);
221 GaleriList.set(
"ny", 4);
222 GaleriList.set(
"nz", 4 * Comm.NumProc());
223 GaleriList.set(
"mx", 1);
224 GaleriList.set(
"my", 1);
225 GaleriList.set(
"mz", Comm.NumProc());
226 Epetra_Map* Map = CreateMap(
"Cartesian3D", Comm, GaleriList);
234 Epetra_CrsMatrix* Matrix =
CreateCrsMatrix(
"Laplace3D", Map, GaleriList);
236 bool transpose = false ;
237 bool distribute = false ;
239 Epetra_CrsMatrix *Matrix = 0 ;
240 Epetra_Map *Map = 0 ;
241 MyCreateCrsMatrix( argv[1], Comm, Map, transpose, distribute, symmetric, Matrix ) ;
246 Epetra_MultiVector LHS(*Map, 1);
247 Epetra_MultiVector RHS(*Map, 1);
250 Epetra_LinearProblem Problem(Matrix, &LHS, &RHS);
257 List.set(
"PrintTiming",
true);
258 List.set(
"PrintStatus",
true);
259 List.set(
"MaxProcs",Comm.NumProc());
261 std::vector<std::string> SolverType;
262 SolverType.push_back(
"Amesos_Paraklete");
263 SolverType.push_back(
"Amesos_Klu");
266 SolverType.push_back(
"Amesos_Lapack");
267 SolverType.push_back(
"Amesos_Umfpack");
268 SolverType.push_back(
"Amesos_Pardiso");
269 SolverType.push_back(
"Amesos_Taucs");
270 SolverType.push_back(
"Amesos_Superlu");
271 SolverType.push_back(
"Amesos_Superludist");
272 SolverType.push_back(
"Amesos_Mumps");
273 SolverType.push_back(
"Amesos_Dscpack");
274 SolverType.push_back(
"Amesos_Scalapack");
276 Epetra_Time Time(Comm);
284 for (
unsigned int i = 0 ; i < SolverType.size() ; ++i)
287 if (Factory.
Query(SolverType[i]))
293 Matrix->Multiply(
false, LHS, RHS);
300 assert (Solver != 0);
310 std::cout << std::endl
311 <<
"Solver " << SolverType[i]
312 <<
", verbose = " <<
verbose << std::endl ;
316 Time.ResetStartTime();
319 std::cout << std::endl
320 <<
"Solver " << SolverType[i]
321 <<
", symbolic factorization time = " 322 << Time.ElapsedTime() << std::endl;
327 std::cout <<
"Solver " << SolverType[i]
328 <<
", numeric factorization time = " 329 << Time.ElapsedTime() << std::endl;
334 std::cout <<
"Solver " << SolverType[i]
336 << Time.ElapsedTime() << std::endl;
342 double d = 0.0, d_tot = 0.0;
343 for (
int j = 0 ; j< LHS.Map().NumMyElements() ; ++j)
344 d += (LHS[0][j] - 1.0) * (LHS[0][j] - 1.0);
346 Comm.SumAll(&d,&d_tot,1);
348 std::cout <<
"Solver " << SolverType[i] <<
", ||x - x_exact||_2 = " 349 << sqrt(d_tot) << std::endl;
354 TotalResidual += d_tot;
361 if (TotalResidual > 1e-9)
368 return(EXIT_SUCCESS);
virtual int Solve()=0
Solves A X = B (or AT x = B)
virtual int NumericFactorization()=0
Performs NumericFactorization on the matrix A.
virtual int SymbolicFactorization()=0
Performs SymbolicFactorization on the matrix A.
virtual int SetParameters(Teuchos::ParameterList &ParameterList)=0
Updates internal variables.
int CrsMatrixTranspose(Epetra_CrsMatrix *In, Epetra_CrsMatrix *Out)
#define AMESOS_CHK_ERR(a)
int main(int argc, char *argv[])
Factory for binding a third party direct solver to an Epetra_LinearProblem.
#define EPETRA_CHK_ERR(xxx)
virtual int SetUseTranspose(bool UseTranspose)=0
If set true, X will be set to the solution of AT X = B (not A X = B)
int CreateCrsMatrix(const char *in_filename, const Epetra_Comm &Comm, Epetra_Map *&readMap, const bool transpose, const bool distribute, bool &symmetric, Epetra_CrsMatrix *&Matrix)
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.
int MyCreateCrsMatrix(const char *in_filename, const Epetra_Comm &Comm, Epetra_Map *&readMap, const bool transpose, const bool distribute, bool &symmetric, Epetra_CrsMatrix *&Matrix)
Amesos_BaseSolver: A pure virtual class for direct solution of real-valued double-precision operators...
bool Query(const char *ClassType)
Queries whether a given interface is available or not.