Commit 42bc9a04 authored by Dmitry Bagaev's avatar Dmitry Bagaev

Some fixes for petsc and MatSolve example

parent 00395d87
......@@ -2,6 +2,7 @@
#include <iostream>
#include <iomanip>
#include <cmath>
#include <cstdio>
#include "inmost.h"
#include "Source/Solver/refactoring/Solver2.h"
......@@ -79,11 +80,20 @@ int main(int argc, char ** argv) {
}
//Matrix file name found with -m or --matrix options
if (strcmp(argv[i], "-m") == 0 || strcmp(argv[i], "--matrix") == 0) {
if (processRank == 0) {
std::cout << "Matrix file found: " << argv[i + 1] << std::endl;
}
matrixFound = true;
matrixFileName = std::string(argv[i + 1]);
FILE *matrixFile = fopen(matrixFileName.c_str(), "r");
if (matrixFile == NULL) {
if (processRank == 0) {
std::cout << "Matrix file not found: " << argv[i + 1] << std::endl;
exit(1);
}
} else {
if (processRank == 0) {
std::cout << "Matrix file found: " << argv[i + 1] << std::endl;
}
}
fclose(matrixFile);
i++;
continue;
}
......
......@@ -18,7 +18,7 @@ namespace INMOST {
this->solver = SolverFactory::getSolver(solverName);
this->prefix = prefix;
solver->setCommunicator(_comm);
solver->SetCommunicator(_comm);
std::string solverDatabasePath = Solver2::parseDatabase(solverName);
solver->Initialize(argc, argv, solverDatabasePath.c_str(), prefix);
}
......@@ -27,14 +27,14 @@ namespace INMOST {
this->solver = SolverFactory::copySolver(other.solver);
this->prefix = other.prefix;
solver->setCommunicator(other.solver->getCommunicator());
std::string solverDatabasePath = Solver2::parseDatabase(solver->getSolverName());
solver->SetCommunicator(other.solver->GetCommunicator());
std::string solverDatabasePath = Solver2::parseDatabase(solver->SolverName());
solver->Initialize(argc, argv, solverDatabasePath.c_str(), prefix);
}
Solver2& Solver2::operator=(const Solver2& other) {
if( this != &other ) {
this->solver->setCommunicator(other.solver->getCommunicator());
this->solver->SetCommunicator(other.solver->GetCommunicator());
this->prefix = other.prefix;
this->solver->Assign(other.solver);
}
......@@ -60,7 +60,7 @@ namespace INMOST {
bool Solver2::Solve(INMOST::Sparse::Vector & RHS, INMOST::Sparse::Vector & SOL) {
if( !solver->isMatrixSet()) throw MatrixNotSetInSolver;
if( RHS.GetCommunicator() != solver->getCommunicator() || SOL.GetCommunicator() != solver->getCommunicator()) throw DifferentCommunicatorInSolver;
if( RHS.GetCommunicator() != solver->GetCommunicator() || SOL.GetCommunicator() != solver->GetCommunicator()) throw DifferentCommunicatorInSolver;
INMOST_DATA_ENUM_TYPE vbeg,vend;
RHS.GetInterval(vbeg,vend);
if( RHS.Size() != SOL.Size() )
......@@ -79,7 +79,6 @@ namespace INMOST {
Sparse::DestroyRowEntryType();
Solver2::is_finalized = true;
Solver2::is_initialized = false;
}
bool Solver2::isInitialized() {
......@@ -102,11 +101,11 @@ namespace INMOST {
return solver->ReturnReason();
}
std::string Solver2::getSolverName() const {
return solver->getSolverName();
std::string Solver2::SolverName() const {
return solver->SolverName();
}
std::string Solver2::getSolverPrefix() const {
std::string Solver2::SolverPrefix() const {
return prefix;
}
......
......@@ -42,8 +42,8 @@ namespace INMOST {
Solver2(const Solver2& other);
Solver2& operator =(const Solver2& other);
std::string getSolverName() const;
std::string getSolverPrefix() const;
std::string SolverName() const;
std::string SolverPrefix() const;
static void Initialize(int *argc, char ***argv, const char *database);
static bool isInitialized();
......
......@@ -19,7 +19,7 @@ namespace INMOST {
}
SolverInterface *SolverFactory::copySolver(const SolverInterface *other) {
auto iterator = SolverFactory::solvers.find(other->getSolverName());
auto iterator = SolverFactory::solvers.find(other->SolverName());
if (iterator != SolverFactory::solvers.end()) {
return iterator->second->copy(other);
} else {
......
......@@ -11,29 +11,38 @@
namespace INMOST {
class SolverInterface {
private:
protected:
INMOST_MPI_Comm communicator;
public:
SolverInterface() {};
SolverInterface(const SolverInterface* other) {};
virtual void Assign(const SolverInterface* other) = 0;
virtual const std::string getSolverName() const = 0;
virtual void Initialize(int *argc, char ***argv, const char *parameters_file, std::string prefix) = 0;
virtual bool Solve(INMOST::Sparse::Vector & RHS, INMOST::Sparse::Vector & SOL) = 0;
virtual void SetMatrix(Sparse::Matrix & A, bool ModifiedPattern, bool OldPreconditioner) = 0;
virtual void Assign(const SolverInterface* other) = 0;
virtual bool Solve(INMOST::Sparse::Vector & RHS, INMOST::Sparse::Vector & SOL) = 0;
virtual bool isMatrixSet() = 0;
virtual INMOST_DATA_REAL_TYPE GetPropertyReal(std::string property) const = 0;
virtual INMOST_DATA_ENUM_TYPE GetPropertyEnum(std::string property) const = 0;
virtual const INMOST_DATA_ENUM_TYPE Iterations() const = 0;
virtual const INMOST_DATA_REAL_TYPE Residual() const = 0;
virtual const std::string ReturnReason() const = 0;
virtual bool isMatrixSet() = 0;
virtual const std::string SolverName() const = 0;
virtual void Finalize() = 0;
virtual ~SolverInterface() {};
void setCommunicator(INMOST_MPI_Comm _communicator) {
void SetCommunicator(INMOST_MPI_Comm _communicator) {
communicator = _communicator;
}
INMOST_MPI_Comm getCommunicator() {
INMOST_MPI_Comm GetCommunicator() {
return communicator;
}
};
......
......@@ -6,7 +6,10 @@
namespace INMOST {
unsigned int SolverPETSc::petscSolversCount = 0;
SolverPETSc::SolverPETSc() {
petscSolversCount++;
this->ksp = NULL;
this->matrix = NULL;
this->rhs = NULL;
......@@ -14,12 +17,13 @@ namespace INMOST {
}
SolverPETSc::SolverPETSc(const SolverInterface *other) {
petscSolversCount++;
const SolverPETSc *solver = static_cast<const SolverPETSc*>(other);
this->ksp = NULL;
this->matrix = NULL;
this->rhs = NULL;
this->solution = NULL;
SolverCopyDataPetsc(&ksp, solver->ksp, this->getCommunicator());
SolverCopyDataPetsc(&ksp, solver->ksp, solver->communicator);
if (solver->matrix != NULL) {
MatrixCopyDataPetsc(&matrix, solver->matrix);
SolverSetMatrixPetsc(ksp, matrix,false,false);
......@@ -30,22 +34,38 @@ namespace INMOST {
VectorCopyDataPetsc(&solution, solver->solution);
}
void SolverPETSc::Initialize(int *argc, char ***argv, const char *parameters_file, std::string prefix) {
this->parameters_file = std::string(parameters_file);
SolverInitializePetsc(argc, argv, parameters_file);
SolverInitDataPetsc(&ksp, this->getCommunicator(), prefix.c_str());
}
void SolverPETSc::Finalize() {
SolverFinalizePetsc();
}
const std::string SolverPETSc::getSolverName() const {
return "petsc";
void SolverPETSc::Assign(const SolverInterface *other) {
const SolverPETSc *other_solver = static_cast<const SolverPETSc*>(other);
this->parametersFile = other_solver->parametersFile;
SolverAssignDataPetsc(ksp, other_solver->ksp);
if (other_solver->matrix != NULL) {
if (matrix != NULL) {
MatrixAssignDataPetsc(matrix, other_solver->matrix);
} else {
MatrixCopyDataPetsc(&matrix, other_solver->matrix);
}
SolverSetMatrixPetsc(ksp, matrix, false, false);
}
if (other_solver->rhs != NULL) {
if (rhs != NULL) {
VectorAssignDataPetsc(rhs, other_solver->rhs);
} else {
VectorCopyDataPetsc(&rhs, other_solver->rhs);
}
}
if (other_solver->solution!= NULL) {
if (solution != NULL) {
VectorAssignDataPetsc(solution, other_solver->solution);
} else {
VectorCopyDataPetsc(&solution, other_solver->solution);
}
}
}
bool SolverPETSc::isMatrixSet() {
return matrix != NULL;
void SolverPETSc::Initialize(int *argc, char ***argv, const char *parameters_file, std::string prefix) {
this->parametersFile = std::string(parameters_file);
SolverInitializePetsc(argc, argv, parameters_file);
SolverInitDataPetsc(&ksp, this->communicator, prefix.c_str());
}
void SolverPETSc::SetMatrix(Sparse::Matrix & A, bool ModifiedPattern, bool OldPreconditioner) {
......@@ -56,10 +76,10 @@ namespace INMOST {
}
INMOST_DATA_ENUM_TYPE mbeg,mend;
A.GetInterval(mbeg,mend);
if( modified_pattern ) {
if( modified_pattern ) {
local_size = A.Size();
#if defined(USE_MPI)
MPI_Allreduce(&local_size,&global_size,1,INMOST_MPI_DATA_ENUM_TYPE,MPI_SUM, this->getCommunicator());
MPI_Allreduce(&local_size,&global_size,1,INMOST_MPI_DATA_ENUM_TYPE,MPI_SUM, this->communicator);
#else
global_size = local_size;
#endif
......@@ -68,18 +88,18 @@ namespace INMOST {
int * diag_nonzeroes = new int[local_size];
int * off_diag_nonzeroes = new int[local_size];
unsigned k = 0;
for(Sparse::Matrix::iterator it = A.Begin(); it != A.End(); ++it) {
diag_nonzeroes[k] = off_diag_nonzeroes[k] = 0;
for(INMOST_DATA_ENUM_TYPE i = 0; i < it->Size(); i++) {
INMOST_DATA_ENUM_TYPE index = it->GetIndex(i);
if( index < mbeg || index > mend-1 ) {
off_diag_nonzeroes[k]++;
} else {
} else {
diag_nonzeroes[k]++;
}
}
if( diag_nonzeroes[k] + off_diag_nonzeroes[k] > max ) {
if( diag_nonzeroes[k] + off_diag_nonzeroes[k] > max ) {
max = diag_nonzeroes[k] + off_diag_nonzeroes[k];
}
k++;
......@@ -130,7 +150,7 @@ namespace INMOST {
}
MatrixFinalizePetsc(matrix);
if (parameters_file == "") {
if (parametersFile == "") {
SolverSetDropTolerancePetsc(ksp, DEFAULT_PRECONDITIONER_DROP_TOLERANCE);
SolverSetFillLevelPetsc(ksp, DEFAULT_PRECONDITIONER_FILL_LEVEL);
SolverSetOverlapPetsc(ksp, DEFAULT_ADDITIVE_SCHWARTZ_OVERLAP);
......@@ -139,33 +159,6 @@ namespace INMOST {
SolverSetMatrixPetsc(ksp, matrix, modified_pattern, OldPreconditioner);
}
void SolverPETSc::Assign(const SolverInterface *other) {
const SolverPETSc *other_solver = static_cast<const SolverPETSc*>(other);
SolverAssignDataPetsc(ksp, other_solver->ksp);
if (other_solver->matrix != NULL) {
if (matrix != NULL) {
MatrixAssignDataPetsc(matrix, other_solver->matrix);
} else {
MatrixCopyDataPetsc(&matrix, other_solver->matrix);
}
SolverSetMatrixPetsc(ksp, matrix, false, false);
}
if (other_solver->rhs != NULL) {
if (rhs != NULL) {
VectorAssignDataPetsc(rhs, other_solver->rhs);
} else {
VectorCopyDataPetsc(&rhs, other_solver->rhs);
}
}
if (other_solver->solution!= NULL) {
if (solution != NULL) {
VectorAssignDataPetsc(solution, other_solver->solution);
} else {
VectorCopyDataPetsc(&solution, other_solver->solution);
}
}
}
bool SolverPETSc::Solve(INMOST::Sparse::Vector &RHS, INMOST::Sparse::Vector &SOL) {
if (rhs == NULL) {
VectorInitDataPetsc(&rhs, RHS.GetCommunicator(), RHS.GetName().c_str());
......@@ -188,7 +181,7 @@ namespace INMOST {
}
VectorFillPetsc(rhs, local_size, positions, values);
VectorFinalizePetsc(rhs);
k = 0;
for(Sparse::Vector::iterator it = SOL.Begin(); it != SOL.End(); ++it) {
values[k] = *it;
......@@ -198,7 +191,7 @@ namespace INMOST {
VectorFinalizePetsc(solution);
}
if(parameters_file == "") {
if(parametersFile == "") {
SolverSetTolerancesPetsc(ksp,
DEFAULT_RELATIVE_TOLERANCE,
DEFAULT_ABSOLUTE_TOLERANCE,
......@@ -217,8 +210,20 @@ namespace INMOST {
}
delete [] positions;
delete [] values;
return result;
return result;
}
bool SolverPETSc::isMatrixSet() {
return matrix != NULL;
}
INMOST_DATA_REAL_TYPE SolverPETSc::GetPropertyReal(std::string property) const {
throw NotImplemented;
}
INMOST_DATA_ENUM_TYPE SolverPETSc::GetPropertyEnum(std::string property) const {
throw NotImplemented;
}
const INMOST_DATA_ENUM_TYPE SolverPETSc::Iterations() const {
......@@ -233,8 +238,18 @@ namespace INMOST {
return std::string(SolverConvergedReasonPetsc(ksp));
}
SolverPETSc::~SolverPETSc() {
const std::string SolverPETSc::SolverName() const {
return "petsc";
}
void SolverPETSc::Finalize() {
if (petscSolversCount == 1) {
SolverFinalizePetsc();
}
}
SolverPETSc::~SolverPETSc() {
petscSolversCount--;
}
}
......@@ -17,7 +17,10 @@ namespace INMOST {
class SolverPETSc : public SolverInterface {
private:
std::string parameters_file;
//TODO
//may be should find another way to count petsc solvers and finalize them
static unsigned int petscSolversCount;
std::string parametersFile;
KSP* ksp;
Mat* matrix;
Vec* rhs;
......@@ -27,16 +30,25 @@ namespace INMOST {
public:
SolverPETSc();
SolverPETSc(const SolverInterface* other);
virtual const std::string getSolverName() const;
virtual void Assign(const SolverInterface* other);
virtual void Initialize(int *argc, char ***argv, const char *parameters_file, std::string prefix);
virtual void SetMatrix(Sparse::Matrix & A, bool ModifiedPattern, bool OldPreconditioner);
virtual bool Solve(INMOST::Sparse::Vector &RHS, INMOST::Sparse::Vector &SOL);
virtual void SetMatrix(Sparse::Matrix & A, bool ModifiedPattern, bool OldPreconditioner);
virtual void Assign(const SolverInterface* other);
virtual const INMOST_DATA_ENUM_TYPE Iterations() const;
virtual const INMOST_DATA_REAL_TYPE Residual() const;
virtual const std::string ReturnReason() const;
virtual bool isMatrixSet();
virtual void Finalize();
virtual INMOST_DATA_REAL_TYPE GetPropertyReal(std::string property) const;
virtual INMOST_DATA_ENUM_TYPE GetPropertyEnum(std::string property) const;
virtual const INMOST_DATA_ENUM_TYPE Iterations() const;
virtual const INMOST_DATA_REAL_TYPE Residual() const;
virtual const std::string ReturnReason() const;
virtual const std::string SolverName() const;
virtual void Finalize();
virtual ~SolverPETSc();
};
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment