#include "Amesos_ConfigDefs.h"
#ifdef HAVE_MPI
#include "mpi.h"
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Amesos.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_Import.h"
#include "Epetra_Export.h"
#include "Epetra_Map.h"
#include "Epetra_MultiVector.h"
#include "Epetra_Vector.h"
#include "Epetra_LinearProblem.h"
#include "Galeri_ReadHB.h"
int main(int argc, char *argv[])
{
#ifdef HAVE_MPI
MPI_Init(&argc, &argv);
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
std::string matrix_file;
std::string solver_type;
if (argc > 1)
matrix_file = argv[1];
else
matrix_file = "662_bus.rsa";
if (argc > 2)
solver_type = argv[2];
else
solver_type = "Klu";
if (Comm.MyPID() == 0)
std::cout << "Reading matrix `" << matrix_file << "'";
Epetra_Map* readMap;
Epetra_CrsMatrix* readA;
Epetra_Vector* readx;
Epetra_Vector* readb;
Epetra_Vector* readxexact;
try
{
Galeri::ReadHB(matrix_file.c_str(), Comm, readMap,
readA, readx, readb, readxexact);
}
catch(...)
{
std::cout << "Caught exception, maybe file name is incorrect" << std::endl;
#ifdef HAVE_MPI
MPI_Finalize();
#else
exit(EXIT_SUCCESS);
#endif
}
Epetra_Map* mapPtr = 0;
if(readMap->GlobalIndicesInt())
mapPtr = new Epetra_Map((int) readMap->NumGlobalElements(), 0, Comm);
else if(readMap->GlobalIndicesLongLong())
mapPtr = new Epetra_Map(readMap->NumGlobalElements(), 0, Comm);
else
assert(false);
Epetra_Map& map = *mapPtr;
Epetra_CrsMatrix A(Copy, map, 0);
const Epetra_Map &OriginalMap = readA->RowMatrixRowMap() ;
assert (OriginalMap.SameAs(*readMap));
Epetra_Export exporter(OriginalMap, map);
Epetra_Vector x(map);
Epetra_Vector b(map);
Epetra_Vector xexact(map);
x.Export(*readx, exporter, Add);
b.Export(*readb, exporter, Add);
xexact.Export(*readxexact, exporter, Add);
A.Export(*readA, exporter, Add);
A.FillComplete();
delete readMap;
delete readA;
delete readx;
delete readb;
delete readxexact;
delete mapPtr;
Epetra_LinearProblem Problem(&A,&x,&b);
if (!Comm.MyPID())
std::cout << "Calling Amesos..." << std::endl;
Solver = Factory.
Create(solver_type,Problem);
if (Solver == 0) {
std::cerr << "Selected solver (" << solver_type << ") is not available" << std::endl;
#ifdef HAVE_MPI
MPI_Finalize();
#endif
return(EXIT_SUCCESS);
}
double residual;
Epetra_Vector Ax(b.Map());
A.Multiply(false, x, Ax);
Ax.Update(1.0, b, -1.0);
Ax.Norm2(&residual);
if (!Comm.MyPID())
std::cout << "After Amesos solution, ||b - A * x||_2 = " << residual << std::endl;
delete Solver;
if (residual > 1e-5)
return(EXIT_FAILURE);
#ifdef HAVE_MPI
MPI_Finalize();
#endif
return(EXIT_SUCCESS);
}