#include "Amesos_ConfigDefs.h"
#ifdef HAVE_MPI
#include "mpi.h"
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Vector.h"
#include "Epetra_Time.h"
#include "Epetra_RowMatrix.h"
#include "Epetra_CrsMatrix.h"
#include "Amesos.h"
#include "Teuchos_ParameterList.hpp"
#include "Galeri_Maps.h"
#include "Galeri_CrsMatrices.h"
using namespace Teuchos;
using namespace Galeri;
#include "Trilinos_Util.h"
#include "Trilinos_Util_ReadMatrixMarket2Epetra.h"
#include "CrsMatrixTranspose.h"
#include "Teuchos_RCP.hpp"
#include "Epetra_Export.h"
int MyCreateCrsMatrix( const char *in_filename, const Epetra_Comm &Comm,
Epetra_Map *& readMap,
const bool transpose, const bool distribute,
bool& symmetric, Epetra_CrsMatrix *& Matrix ) {
Epetra_CrsMatrix * readA = 0;
Epetra_Vector * readx = 0;
Epetra_Vector * readb = 0;
Epetra_Vector * readxexact = 0;
FILE *in_file = fopen( in_filename, "r");
const char *filename;
if (in_file == NULL )
filename = &in_filename[1] ;
else {
filename = in_filename ;
fclose( in_file );
}
symmetric = false ;
std::string FileName (filename);
int FN_Size = FileName.size() ;
std::string LastFiveBytes = FileName.substr( EPETRA_MAX(0,FN_Size-5), FN_Size );
std::string LastFourBytes = FileName.substr( EPETRA_MAX(0,FN_Size-4), FN_Size );
if ( LastFiveBytes == ".TimD" ) {
EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename, false, Comm, readMap, readA, readx,
readb, readxexact, false, true, true ) );
symmetric = false;
} else {
if ( LastFiveBytes == ".triU" ) {
EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename, false, Comm, readMap, readA, readx,
readb, readxexact) );
symmetric = false;
} else {
if ( LastFiveBytes == ".triS" ) {
EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename, true, Comm, readMap, readA, readx,
readb, readxexact) );
symmetric = true;
} else {
if ( LastFourBytes == ".mtx" ) {
EPETRA_CHK_ERR( Trilinos_Util_ReadMatrixMarket2Epetra( filename, Comm, readMap,
readA, readx, readb, readxexact) );
FILE* in_file = fopen( filename, "r");
assert (in_file != NULL) ;
const int BUFSIZE = 800 ;
char buffer[BUFSIZE] ;
fgets( buffer, BUFSIZE, in_file ) ;
std::string headerline1 = buffer;
#ifdef TFLOP
if ( headerline1.find("symmetric") < BUFSIZE ) symmetric = true;
#else
if ( headerline1.find("symmetric") != std::string::npos) symmetric = true;
#endif
fclose(in_file);
} else {
Trilinos_Util_ReadHb2Epetra( filename, Comm, readMap, readA, readx,
readb, readxexact) ;
if ( LastFourBytes == ".rsa" ) symmetric = true ;
}
}
}
}
if ( readb ) delete readb;
if ( readx ) delete readx;
if ( readxexact ) delete readxexact;
Epetra_CrsMatrix *serialA ;
Epetra_CrsMatrix *transposeA;
if ( transpose ) {
transposeA = new Epetra_CrsMatrix( Copy, *readMap, 0 );
assert( CrsMatrixTranspose( readA, transposeA ) == 0 );
serialA = transposeA ;
delete readA;
readA = 0 ;
} else {
serialA = readA ;
}
assert( (void *) &serialA->Graph() ) ;
assert( (void *) &serialA->RowMap() ) ;
assert( serialA->RowMap().SameAs(*readMap) ) ;
if ( distribute ) {
Epetra_Map DistMap(readMap->NumGlobalElements(), 0, Comm);
Epetra_Export exporter( *readMap, DistMap );
Epetra_CrsMatrix *Amat = new Epetra_CrsMatrix( Copy, DistMap, 0 );
Amat->Export(*serialA, exporter, Add);
assert(Amat->FillComplete()==0);
Matrix = Amat;
delete readMap;
readMap = 0 ;
delete serialA;
} else {
Matrix = serialA;
}
return 0;
}
int main(int argc, char *argv[])
{
#ifdef HAVE_MPI
MPI_Init(&argc, &argv);
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
bool verbose = (Comm.MyPID() == 0);
double TotalResidual = 0.0;
#ifndef FILENAME_SPECIFIED_ON_COMMAND_LINE
ParameterList GaleriList;
GaleriList.set("nx", 4);
GaleriList.set("ny", 4);
GaleriList.set("nz", 4 * Comm.NumProc());
GaleriList.set("mx", 1);
GaleriList.set("my", 1);
GaleriList.set("mz", Comm.NumProc());
Epetra_Map* Map = CreateMap("Cartesian3D", Comm, GaleriList);
Epetra_CrsMatrix* Matrix = CreateCrsMatrix("Laplace3D", Map, GaleriList);
#else
bool transpose = false ;
bool distribute = false ;
bool symmetric ;
Epetra_CrsMatrix *Matrix = 0 ;
Epetra_Map *Map = 0 ;
MyCreateCrsMatrix( argv[1], Comm, Map, transpose, distribute, symmetric, Matrix ) ;
#endif
Epetra_MultiVector LHS(*Map, 1);
Epetra_MultiVector RHS(*Map, 1);
Epetra_LinearProblem Problem(Matrix, &LHS, &RHS);
ParameterList List;
List.set("PrintTiming",true);
List.set("PrintStatus",true);
List.set("MaxProcs",Comm.NumProc());
std::vector<std::string> SolverType;
SolverType.push_back("Amesos_Paraklete");
SolverType.push_back("Amesos_Klu");
Comm.Barrier() ;
#if 1
SolverType.push_back("Amesos_Lapack");
SolverType.push_back("Amesos_Umfpack");
SolverType.push_back("Amesos_Pardiso");
SolverType.push_back("Amesos_Taucs");
SolverType.push_back("Amesos_Superlu");
SolverType.push_back("Amesos_Superludist");
SolverType.push_back("Amesos_Mumps");
SolverType.push_back("Amesos_Dscpack");
SolverType.push_back("Amesos_Scalapack");
#endif
Epetra_Time Time(Comm);
for (unsigned int i = 0 ; i < SolverType.size() ; ++i)
{
if (Factory.
Query(SolverType[i]))
{
LHS.PutScalar(1.0);
Matrix->Multiply(false, LHS, RHS);
LHS.Random();
assert (Solver != 0);
Comm.Barrier() ;
if (verbose)
std::cout << std::endl
<< "Solver " << SolverType[i]
<< ", verbose = " << verbose << std::endl ;
Comm.Barrier() ;
Time.ResetStartTime();
if (verbose)
std::cout << std::endl
<< "Solver " << SolverType[i]
<< ", symbolic factorization time = "
<< Time.ElapsedTime() << std::endl;
Comm.Barrier() ;
if (verbose)
std::cout << "Solver " << SolverType[i]
<< ", numeric factorization time = "
<< Time.ElapsedTime() << std::endl;
Comm.Barrier() ;
AMESOS_CHK_ERR(Solver->
Solve());
if (verbose)
std::cout << "Solver " << SolverType[i]
<< ", solve time = "
<< Time.ElapsedTime() << std::endl;
Comm.Barrier() ;
double d = 0.0, d_tot = 0.0;
for (int j = 0 ; j< LHS.Map().NumMyElements() ; ++j)
d += (LHS[0][j] - 1.0) * (LHS[0][j] - 1.0);
Comm.SumAll(&d,&d_tot,1);
if (verbose)
std::cout << "Solver " << SolverType[i] << ", ||x - x_exact||_2 = "
<< sqrt(d_tot) << std::endl;
delete Solver;
TotalResidual += d_tot;
}
}
delete Matrix;
delete Map;
if (TotalResidual > 1e-9)
exit(EXIT_FAILURE);
#ifdef HAVE_MPI
MPI_Finalize();
#endif
return(EXIT_SUCCESS);
}
Amesos_BaseSolver: A pure virtual class for direct solution of real-valued double-precision operators...
Definition: Amesos_BaseSolver.h:225
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 SetUseTranspose(bool UseTranspose)=0
If set true, X will be set to the solution of AT X = B (not A X = B)
virtual int SetParameters(Teuchos::ParameterList &ParameterList)=0
Updates internal variables.
Factory for binding a third party direct solver to an Epetra_LinearProblem.
Definition: Amesos.h:44
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.
Definition: Amesos.cpp:69
bool Query(const char *ClassType)
Queries whether a given interface is available or not.
Definition: Amesos.cpp:193