Commit 00395d87 authored by Dmitry Bagaev's avatar Dmitry Bagaev

New version of Solver, just PETSc, others in progress

parent 5156345a
......@@ -2,3 +2,7 @@
Tests/solver_test001/matrices
SyncToy*
build*/
*.iml
*.xml
......@@ -4,6 +4,7 @@
#include <cmath>
#include "inmost.h"
#include "Source/Solver/refactoring/Solver2.h"
#include "inner_parser.h"
using namespace INMOST;
......@@ -177,8 +178,11 @@ int main(int argc, char ** argv) {
break;
}
Solver::Initialize(&argc, &argv, parametersFound ? parametersFileName.c_str()
: NULL); // Initialize the linear solver in accordance with args
// Initialize the linear solver in accordance with args
Solver2::Initialize(&argc, &argv, parametersFound ? parametersFileName.c_str() : NULL);
Solver2 solver = Solver2("petsc");
Solver2 solver2 = solver;
//solver2.Finalize();
if (processRank == 0) {
std::cout << "Solving with " << Solver::TypeName(type) << std::endl;
......@@ -212,7 +216,7 @@ int main(int argc, char ** argv) {
bool success;
double resid, realresid = 0;
std::string reason;
Solver s(type); // Declare the linear solver by specified type
//Solver s(type); // Declare the linear solver by specified type
if (parametersFound) {
char *fileName = findInnerOptions(parametersFileName.c_str());
......@@ -222,9 +226,9 @@ int main(int argc, char ** argv) {
for (unsigned ii = 0; ii < options->options.size(); ii++) {
InnerOption *option = options->options[ii];
if (option->type == ENUM) {
s.SetParameterEnum(option->name, (unsigned int) atoi(option->value.c_str()));
//s.SetParameterEnum(option->name, (unsigned int) atoi(option->value.c_str()));
} else {
s.SetParameterReal(option->name, atof(option->value.c_str()));
//s.SetParameterReal(option->name, atof(option->value.c_str()));
};
}
delete options;
......@@ -253,17 +257,18 @@ int main(int argc, char ** argv) {
tempTimer = Timer();
s.SetMatrix(mat); // Compute the preconditioner for the original matrix
solver2.SetMatrix(mat);
//s.SetMatrix(mat); // Compute the preconditioner for the original matrix
BARRIER
if (processRank == 0) std::cout << "preconditioner time: " << Timer() - tempTimer << std::endl;
tempTimer = Timer();
success = s.Solve(b, x); // Solve the linear system with the previously computted preconditioner
success = solver2.Solve(b, x); // Solve the linear system with the previously computted preconditioner
BARRIER
solvingTimer = Timer() - solvingTimer;
if (processRank == 0) std::cout << "iterations time: " << Timer() - tempTimer << std::endl;
iters = s.Iterations(); // Get the number of iterations performed
resid = s.Residual(); // Get the final residual achieved
reason = s.GetReason(); // Get the convergence reason
iters = solver2.Iterations(); // Get the number of iterations performed
resid = solver2.Residual(); // Get the final residual achieved
reason = solver2.ReturnReason(); // Get the convergence reason
//x.Save("output.sol"); // Save the solution if required
// Compute the true residual
......@@ -347,6 +352,6 @@ int main(int argc, char ** argv) {
}
}
BARRIER
Solver::Finalize(); // Finalize solver and close MPI activity
Solver2::Finalize(); // Finalize solver and close MPI activity
return 0;
}
......@@ -221,6 +221,7 @@ namespace INMOST
IndexesAreDifferentInSolver,
PrepareMatrixFirst,
CannotReusePreconditionerOfDifferentSize,
SolverNotFound,
/// The list of errors may occur in the Partitioner.
ErrorInPartitioner = 500,
......
add_subdirectory(refactoring)
set(SOURCE
${SOURCE}
${CMAKE_CURRENT_SOURCE_DIR}/solver.cpp
......
if(USE_SOLVER_PETSC)
add_subdirectory(solver_petsc)
endif()
set(HEADER
${HEADER}
${CMAKE_CURRENT_SOURCE_DIR}/SolverInterface.h
${CMAKE_CURRENT_SOURCE_DIR}/SolverFactory.h
${CMAKE_CURRENT_SOURCE_DIR}/Solver2.h)
set(SOURCE
${SOURCE}
${CMAKE_CURRENT_SOURCE_DIR}/Solver2.cpp
${CMAKE_CURRENT_SOURCE_DIR}/SolverFactory.cpp)
set(HEADER ${HEADER} PARENT_SCOPE)
set(SOURCE ${SOURCE} PARENT_SCOPE)
\ No newline at end of file
//
// Created by Dmitri Bagaev on 22.09.16.
//
#include <Source/Solver/refactoring/solver_petsc/SolverPETSc.h>
#include <inmost_sparse.h>
#include "Solver2.h"
namespace INMOST {
int *Solver2::argc = NULL;
char ***Solver2::argv = NULL;
const char *Solver2::database = NULL;
bool Solver2::is_initialized = false;
bool Solver2::is_finalized = false;
Solver2::Solver2(std::string solverName, std::string prefix, INMOST_MPI_Comm _comm) {
this->solver = SolverFactory::getSolver(solverName);
this->prefix = prefix;
solver->setCommunicator(_comm);
std::string solverDatabasePath = Solver2::parseDatabase(solverName);
solver->Initialize(argc, argv, solverDatabasePath.c_str(), prefix);
}
Solver2::Solver2(const Solver2 &other) {
this->solver = SolverFactory::copySolver(other.solver);
this->prefix = other.prefix;
solver->setCommunicator(other.solver->getCommunicator());
std::string solverDatabasePath = Solver2::parseDatabase(solver->getSolverName());
solver->Initialize(argc, argv, solverDatabasePath.c_str(), prefix);
}
Solver2& Solver2::operator=(const Solver2& other) {
if( this != &other ) {
this->solver->setCommunicator(other.solver->getCommunicator());
this->prefix = other.prefix;
this->solver->Assign(other.solver);
}
return *this;
}
void Solver2::Initialize(int *argc, char ***argv, const char *database) {
Solver2::argc = argc;
Solver2::argv = argv;
Solver2::database = database;
Solver2::is_initialized = true;
Solver2::is_finalized = false;
Sparse::CreateRowEntryType();
//Register all available solvers
#if defined(USE_SOLVER_PETSC)
SolverFactory::registerSolver<SolverPETSc>("petsc");
#endif
}
void Solver2::SetMatrix(Sparse::Matrix & A, bool ModifiedPattern, bool OldPreconditioner) {
solver->SetMatrix(A, ModifiedPattern, OldPreconditioner);
}
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;
INMOST_DATA_ENUM_TYPE vbeg,vend;
RHS.GetInterval(vbeg,vend);
if( RHS.Size() != SOL.Size() )
{
if( SOL.Size() == 0 )
{
SOL.SetInterval(vbeg,vend);
for(Sparse::Vector::iterator ri = SOL.Begin(); ri != SOL.End(); ++ri) *ri = 0.0;
}
else throw InconsistentSizesInSolver;
}
return solver->Solve(RHS, SOL);
}
void Solver2::Finalize() {
Sparse::DestroyRowEntryType();
Solver2::is_finalized = true;
Solver2::is_initialized = false;
}
bool Solver2::isInitialized() {
return is_initialized;
}
bool Solver2::isFinalized() {
return is_finalized;
}
const INMOST_DATA_ENUM_TYPE Solver2::Iterations() const {
return solver->Iterations();
}
const INMOST_DATA_REAL_TYPE Solver2::Residual() const {
return solver->Residual();
}
const std::string Solver2::ReturnReason() const {
return solver->ReturnReason();
}
std::string Solver2::getSolverName() const {
return solver->getSolverName();
}
std::string Solver2::getSolverPrefix() const {
return prefix;
}
Solver2::~Solver2() {
solver->Finalize();
delete solver;
}
std::string Solver2::parseDatabase(std::string solverName) {
const char *name = solverName.c_str();
if( database != NULL ) {
FILE * f = fopen(database, "r");
if (f != NULL) {
char str[4096];
while( !feof(f) && fgets(str,4096,f)) {
int k = 0, l;
for(k = 0; k < (int)strlen(str); ++k) {
if( str[k] == ':' ) break;
}
if( k == strlen(str) ) continue; //invalid line
for(l = 0; l < k; ++l) str[l] = tolower(str[l]);
l = (int)strlen(str)-1; // Right-trim string
while(l > 0 && isspace(str[l]) ) --l;
str[l+1] = 0;
l = k+1;
while(l < (int)strlen(str) && isspace(str[l]) ) ++l;
if( l == strlen(str) ) continue; //skip empty entry
if( !strncmp(str, name, k) ) {
return std::string(str+l);
}
}
fclose(f);
}
}
return std::string("");
}
}
//
// Created by Dmitri Bagaev on 22.09.16.
//
#ifndef INMOST_SOLVERCONTAINER_H
#define INMOST_SOLVERCONTAINER_H
#include <string>
#include "SolverInterface.h"
#define DEFAULT_ADDITIVE_SCHWARTZ_OVERLAP 1
#define DEFAULT_ABSOLUTE_TOLERANCE 1.0e-5
#define DEFAULT_RELATIVE_TOLERANCE 1.0e-12
#define DEFAULT_DIVERGENCE_TOLERANCE 1.0e+100
#define DEFAULT_MAXIMUM_ITERATIONS 2500
#define DEFAULT_SOLVER_GMRES_SUBSTEPS 2
#define DEFAULT_PRECONDITIONER_DROP_TOLERANCE 0.005
#define DEFAULT_PRECONDITIONER_REUSE_TOLERANCE 0.00005
#define DEFAULT_PRECONDITIONER_FILL_LEVEL 3
#define DEFAULT_PRECONDITIONER_DDPQ_TOLERANCE 0.75
#define DEFAULT_PRECONDITIONER_REORDER_NONZEROS 1
#define DEFAULT_PRECONDITIONER_RESCALE_ITERS 6
#define DEFAULT_PRECONDITIONER_CONDITION_ESTIMATION 1
#define DEFAULT_PRECONDITIONER_ADAPT_DDPQ_TOLERANCE 1
namespace INMOST {
class Solver2 {
private:
static int *argc;
static char ***argv;
static const char *database;
static bool is_initialized;
static bool is_finalized;
//Actual solver using for solving system
SolverInterface *solver;
std::string prefix;
public:
Solver2(std::string solverName, std::string prefix = "", INMOST_MPI_Comm _comm = INMOST_MPI_COMM_WORLD);
Solver2(const Solver2& other);
Solver2& operator =(const Solver2& other);
std::string getSolverName() const;
std::string getSolverPrefix() const;
static void Initialize(int *argc, char ***argv, const char *database);
static bool isInitialized();
static bool isFinalized();
static void Finalize();
void SetMatrix(Sparse::Matrix & A, bool ModifiedPattern = true, bool OldPreconditioner = false);
bool Solve(INMOST::Sparse::Vector & RHS, INMOST::Sparse::Vector & SOL);
const INMOST_DATA_ENUM_TYPE Iterations() const;
const INMOST_DATA_REAL_TYPE Residual() const;
const std::string ReturnReason() const;
~Solver2();
static std::string parseDatabase(std::string solverName);
};
}
#endif //INMOST_SOLVERCONTAINER_H
//
// Created by Dmitri Bagaev on 28.09.16.
//
#include "SolverFactory.h"
#include "SolverInterface.h"
namespace INMOST {
std::map<std::string, SolverBaseFactory *> SolverFactory::solvers = std::map<std::string, SolverBaseFactory *>();
SolverInterface *SolverFactory::getSolver(std::string name) {
auto iterator = SolverFactory::solvers.find(name);
if (iterator != SolverFactory::solvers.end()) {
return iterator->second->create();
} else {
throw INMOST::SolverNotFound;
}
}
SolverInterface *SolverFactory::copySolver(const SolverInterface *other) {
auto iterator = SolverFactory::solvers.find(other->getSolverName());
if (iterator != SolverFactory::solvers.end()) {
return iterator->second->copy(other);
} else {
throw INMOST::SolverNotFound;
}
}
std::vector<std::string> SolverFactory::getAvailableSolvers() {
std::vector<std::string> s;
for (auto bi = SolverFactory::solvers.begin(); bi != SolverFactory::solvers.end(); bi++) {
s.push_back(bi->first);
}
return s;
}
bool SolverFactory::isSolverAvailable(std::string name) {
return SolverFactory::solvers.find(name) != SolverFactory::solvers.end();
}
}
//
// Created by Dmitri Bagaev on 28.09.16.
//
#ifndef INMOST_SOLVERFACTORY_H
#define INMOST_SOLVERFACTORY_H
#include "SolverInterface.h"
namespace INMOST {
struct SolverBaseFactory {
virtual SolverInterface* create() = 0;
virtual SolverInterface* copy(const SolverInterface* other) = 0;
virtual ~SolverBaseFactory() {};
};
template <class C> struct SolverCreateFactory : SolverBaseFactory {
SolverInterface *create() {
return new C;
};
SolverInterface *copy(const SolverInterface* other) {
return new C(other);
};
};
class SolverFactory {
private:
static std::map<std::string, SolverBaseFactory *> solvers;
public:
template<class T> static void registerSolver(std::string name) {
solvers.insert(std::make_pair(name, new SolverCreateFactory<T>));
};
static SolverInterface *getSolver(std::string name);
static SolverInterface *copySolver(const SolverInterface *other);
static std::vector<std::string> getAvailableSolvers();
static bool isSolverAvailable(std::string name);
};
}
#endif //INMOST_SOLVERFACTORY_H
//
// Created by Dmitri Bagaev on 22.09.16.
//
#ifndef INMOST_SOLVERINTERFACE_H
#define INMOST_SOLVERINTERFACE_H
#include <string>
#include "inmost_sparse.h"
namespace INMOST {
class SolverInterface {
private:
INMOST_MPI_Comm communicator;
public:
SolverInterface() {};
SolverInterface(const SolverInterface* other) {};
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 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 void Finalize() = 0;
virtual ~SolverInterface() {};
void setCommunicator(INMOST_MPI_Comm _communicator) {
communicator = _communicator;
}
INMOST_MPI_Comm getCommunicator() {
return communicator;
}
};
}
#endif //INMOST_SOLVERINTERFACE_H
set(SOURCE ${SOURCE} ${CMAKE_CURRENT_SOURCE_DIR}/new_solver_petsc.cpp)
set(SOURCE ${SOURCE} ${CMAKE_CURRENT_SOURCE_DIR}/new_solver_petsc.h)
set(SOURCE ${SOURCE} ${CMAKE_CURRENT_SOURCE_DIR}/SolverPETSC.cpp)
set(HEADER ${HEADER} ${CMAKE_CURRENT_SOURCE_DIR}/SolverPETSc.h)
set(SOURCE ${SOURCE} PARENT_SCOPE)
set(HEADER ${HEADER} PARENT_SCOPE)
\ No newline at end of file
//
// Created by Dmitri Bagaev on 22.09.16.
//
#include "SolverPETSc.h"
namespace INMOST {
SolverPETSc::SolverPETSc() {
this->ksp = NULL;
this->matrix = NULL;
this->rhs = NULL;
this->solution = NULL;
}
SolverPETSc::SolverPETSc(const SolverInterface *other) {
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());
if (solver->matrix != NULL) {
MatrixCopyDataPetsc(&matrix, solver->matrix);
SolverSetMatrixPetsc(ksp, matrix,false,false);
}
if (solver->rhs != NULL)
VectorCopyDataPetsc(&rhs, solver->rhs);
if (solver->solution != NULL)
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";
}
bool SolverPETSc::isMatrixSet() {
return matrix != NULL;
}
void SolverPETSc::SetMatrix(Sparse::Matrix & A, bool ModifiedPattern, bool OldPreconditioner) {
bool modified_pattern = ModifiedPattern;
if( !isMatrixSet()) {
MatrixInitDataPetsc(&matrix, A.GetCommunicator(),A.GetName().c_str());
modified_pattern = true;
}
INMOST_DATA_ENUM_TYPE mbeg,mend;
A.GetInterval(mbeg,mend);
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());
#else
global_size = local_size;
#endif
int max = 0;
{
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 {
diag_nonzeroes[k]++;
}
}
if( diag_nonzeroes[k] + off_diag_nonzeroes[k] > max ) {
max = diag_nonzeroes[k] + off_diag_nonzeroes[k];
}
k++;
}
MatrixPreallocatePetsc(matrix, local_size, global_size, diag_nonzeroes, off_diag_nonzeroes);
delete [] diag_nonzeroes;
delete [] off_diag_nonzeroes;
}
if( max > 0 ) {
int * col_positions = new int[max];
double * col_values = new double[max];
unsigned k = 0, m;
for(Sparse::Matrix::iterator it = A.Begin(); it != A.End(); ++it) {
m = 0;
for(INMOST_DATA_ENUM_TYPE i = 0; i < it->Size(); i++) {
col_positions[m] = it->GetIndex(i);
col_values[m] = it->GetValue(i);
m++;
}
MatrixFillPetsc(matrix, mbeg + k, m, col_positions, col_values);
k++;
}
delete [] col_positions;
delete [] col_values;
}
} else {
unsigned max = 0;
for(Sparse::Matrix::iterator it = A.Begin(); it != A.End(); ++it) {
if( it->Size() > max ) {
max = it->Size();
}
}
int * col_positions = new int[max];
double * col_values = new double[max];
unsigned k = 0, m;
for(Sparse::Matrix::iterator it = A.Begin(); it != A.End(); ++it) {
m = 0;
for(INMOST_DATA_ENUM_TYPE i = 0; i < it->Size(); i++) {
col_positions[m] = it->GetIndex(i);
col_values[m] = it->GetValue(i);
m++;
}
MatrixFillPetsc(matrix, mbeg + k, m, col_positions, col_values);
k++;
}
delete [] col_positions;
delete [] col_values;
}
MatrixFinalizePetsc(matrix);
if (parameters_file == "") {
SolverSetDropTolerancePetsc(ksp, DEFAULT_PRECONDITIONER_DROP_TOLERANCE);
SolverSetFillLevelPetsc(ksp, DEFAULT_PRECONDITIONER_FILL_LEVEL);
SolverSetOverlapPetsc(ksp, DEFAULT_ADDITIVE_SCHWARTZ_OVERLAP);
}
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());
}
VectorPreallocatePetsc(rhs,local_size,global_size);
if (solution == NULL) {
VectorInitDataPetsc(&solution, SOL.GetCommunicator(), SOL.GetName().c_str());
}
VectorPreallocatePetsc(solution, local_size, global_size);
INMOST_DATA_ENUM_TYPE vbeg,vend;
RHS.GetInterval(vbeg,vend);
int * positions = new int[local_size];
double * values = new double[local_size];
{
unsigned k = 0;
for(Sparse::Vector::iterator it = RHS.Begin(); it != RHS.End(); ++it) {
positions[k] = vbeg+k;
values[k] = *it;
k++;
}
VectorFillPetsc(rhs, local_size, positions, values);
VectorFinalizePetsc(rhs);
k = 0;
for(Sparse::Vector::iterator it = SOL.Begin(); it != SOL.End(); ++it) {
values[k] = *it;
k++;
}
VectorFillPetsc(solution, local_size, positions, values);
VectorFinalizePetsc(solution);
}
if(parameters_file == "") {
SolverSetTolerancesPetsc(ksp,
DEFAULT_RELATIVE_TOLERANCE,
DEFAULT_ABSOLUTE_TOLERANCE,
DEFAULT_DIVERGENCE_TOLERANCE,
DEFAULT_MAXIMUM_ITERATIONS);
}
bool result = SolverSolvePetsc(ksp, rhs, solution);
if (result) {
VectorLoadPetsc(solution, local_size, positions, values);
unsigned k = 0;
for(Sparse::Vector::iterator it = SOL.Begin(); it != SOL.End(); ++it) {
*it = values[k];
k++;
}
}
delete [] positions;
delete [] values;
return result;
}
const INMOST_DATA_ENUM_TYPE SolverPETSc::Iterations() const {
return SolverIterationNumberPetsc(ksp);