Skip to content

GitLab

Commit ed7ce804 authored by Dmitry Bagaev's avatar Dmitry Bagaev
Browse files

FCBIILU2 refactored using SolverInterface

parent e6fbb7d2
Hide whitespace changes
Inline Side-by-side
Source/Solver/CMakeLists.txt
View file @ ed7ce804
......@@ -8,7 +8,7 @@ set(SOURCE
)
add_subdirectory(solver_inner)
#add_subdirectory(solver_fcbiilu2)
add_subdirectory(solver_fcbiilu2)
add_subdirectory(solver_k3biilu2)
if (USE_SOLVER_PETSC)
......
Source/Solver/solver.cpp
View file @ ed7ce804
......@@ -25,10 +25,14 @@
#include "solver_superlu/SolverSUPERLU.h"
#endif
#if defined (HAVE_SOLVER_K3BIILU2)
#if defined(HAVE_SOLVER_K3BIILU2)
#include "solver_k3biilu2/SolverK3BIILU2.h"
#endif
#if defined(HAVE_SOLVER_FCBIILU2)
#include "solver_fcbiilu2/SolverFCBIILU2.h"
#endif
namespace INMOST {
int *Solver::argc = NULL;
......@@ -116,6 +120,9 @@ namespace INMOST {
#endif
#if defined(HAVE_SOLVER_K3BIILU2)
SolverMaster::registerSolver<SolverK3BIILU2>("k3biilu2");
#endif
#if defined(HAVE_SOLVER_FCBIILU2)
SolverMaster::registerSolver<SolverFCBIILU2>("fcbiilu2");
#endif
Sparse::CreateRowEntryType();
}
......
Source/Solver/solver_fcbiilu2/CMakeLists.txt 0 → 100644
View file @ ed7ce804
if (EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/solver_fcbiilu2.h" AND EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/solver_fcbiilu2.cpp" AND EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/fcbiilu2.cpp")
set(SOLVER_DEFINITIONS ${SOLVER_DEFINITIONS} -DHAVE_SOLVER_FCBIILU2 PARENT_SCOPE)
set(HAVE_SOLVER_FCBIILU2 TRUE PARENT_SCOPE)
set(SOURCE ${SOURCE}
${CMAKE_CURRENT_SOURCE_DIR}/solver_fcbiilu2.cpp
${CMAKE_CURRENT_SOURCE_DIR}/fcbiilu2.cpp
${CMAKE_CURRENT_SOURCE_DIR}/SolverFCBIILU2.cpp
PARENT_SCOPE)
set(HEADER ${HEADER}
${CMAKE_CURRENT_SOURCE_DIR}/solver_fcbiilu2.h
${CMAKE_CURRENT_SOURCE_DIR}/SolverFCBIILU2.h
PARENT_SCOPE)
else ()
set(HAVE_SOLVER_FCBIILU2 FALSE PARENT_SCOPE)
endif ()
\ No newline at end of file
Source/Solver/solver_fcbiilu2/SolverFCBIILU2.cpp 0 → 100644
View file @ ed7ce804
#include "SolverFCBIILU2.h"
namespace INMOST {
SolverFCBIILU2::SolverFCBIILU2() {
}
SolverFCBIILU2::SolverFCBIILU2(const SolverInterface *other) {
const SolverFCBIILU2 *fcother = static_cast<const SolverFCBIILU2 *>(other);
SolverCopyDataFcbiilu2(&solver_data, fcother->solver_data, communicator);
if (fcother->matrix_data != NULL) {
MatrixCopyDataFcbiilu2(&matrix_data, fcother->matrix_data);
SolverSetMatrixFcbiilu2(solver_data, matrix_data, false, false);
}
}
void SolverFCBIILU2::Assign(const SolverInterface *other) {
const SolverFCBIILU2 *fcother = static_cast<const SolverFCBIILU2 *>(other);
SolverAssignDataFcbiilu2(solver_data, fcother->solver_data);
if (fcother->matrix_data != NULL) {
if (matrix_data != NULL) {
MatrixAssignDataFcbiilu2(matrix_data, fcother->matrix_data);
} else {
MatrixCopyDataFcbiilu2(&matrix_data, fcother->matrix_data);
}
SolverSetMatrixFcbiilu2(solver_data, matrix_data, false, false);
}
}
void SolverFCBIILU2::Initialize(int *argc, char ***argv, const char *parameters_file, std::string prefix) {
SolverInitDataFcbiilu2(&solver_data, communicator, prefix.c_str());
SolverInitializeFcbiilu2(argc, argv, parameters_file);
}
void SolverFCBIILU2::SetMatrix(Sparse::Matrix &A, bool ModifiedPattern, bool OldPreconditioner) {
bool modified_pattern = ModifiedPattern;
//~ if( A.comm != comm ) throw DifferentCommunicatorInSolver;
if (matrix_data == NULL) {
MatrixInitDataFcbiilu2(&matrix_data, A.GetCommunicator(), A.GetName().c_str());
modified_pattern = true;
}
if (modified_pattern) {
global_size = local_size = A.Size();
MatrixDestroyDataFcbiilu2(&matrix_data);
MatrixInitDataFcbiilu2(&matrix_data, A.GetCommunicator(), A.GetName().c_str());
INMOST_DATA_ENUM_TYPE nnz = 0, k = 0, q = 1, shift = 1;
int nproc, myid;
#if defined(USE_MPI)
MPI_Comm_size(A.GetCommunicator(), &nproc);
MPI_Comm_rank(A.GetCommunicator(), &myid);
#else
nproc = 1;
myid = 0;
#endif
int *ibl = (int *) malloc(sizeof(int) * (nproc + 1));
ibl[0] = 0;
int n = A.Size();
#if defined(USE_MPI)
INMOST_DATA_ENUM_TYPE mbeg, mend;
A.GetInterval(mbeg, mend);
n = mend - mbeg;
int block_end = mend;
MPI_Allgather(&block_end, 1, MPI_INT, &ibl[1], 1, MPI_INT, A.GetCommunicator());
#else
ibl[1] = n;
#endif
int *ia = (int *) malloc(sizeof(int) * (n + 1));
for (Sparse::Matrix::iterator it = A.Begin(); it != A.End(); it++) nnz += it->Size();
int *ja = (int *) malloc(sizeof(int) * nnz);
double *values = (double *) malloc(sizeof(double) * nnz);
ia[0] = shift;
for (Sparse::Matrix::iterator it = A.Begin(); it != A.End(); it++) {
for (Sparse::Row::iterator jt = it->Begin(); jt != it->End(); jt++) {
ja[k] = jt->first + 1;
values[k] = jt->second;
//std::cout<<"# q="<<q<<" k="<<k<<" ja="<<ja[k]<<" a="<<values[k]<<std::endl;//db!
k++;
}
shift += it->Size();
ia[q++] = shift;
}
MatrixFillFcbiilu2(matrix_data, n, nproc, ibl, ia, ja, values);
} else {
INMOST_DATA_ENUM_TYPE nnz = 0, k = 0;
for (Sparse::Matrix::iterator it = A.Begin(); it != A.End(); it++) nnz += it->Size();
double *values = (double *) malloc(sizeof(double) * nnz);
k = 0;
for (Sparse::Matrix::iterator it = A.Begin(); it != A.End(); it++)
for (Sparse::Row::iterator jt = it->Begin(); jt != it->End(); jt++)
values[k++] = jt->second;
MatrixFillValuesFcbiilu2(matrix_data, values);
}
MatrixFinalizeFcbiilu2(matrix_data);
SolverSetMatrixFcbiilu2(solver_data, matrix_data, modified_pattern, OldPreconditioner);
}
bool SolverFCBIILU2::Solve(INMOST::Sparse::Vector &RHS, INMOST::Sparse::Vector &SOL) {
INMOST_DATA_ENUM_TYPE vbeg, vend;
RHS.GetInterval(vbeg, vend);
void *rhs_data = NULL;
VectorInitDataFcbiilu2(&rhs_data, RHS.GetCommunicator(), RHS.GetName().c_str());
VectorPreallocateFcbiilu2(rhs_data, local_size);
void *solution_data = NULL;
VectorInitDataFcbiilu2(&solution_data, SOL.GetCommunicator(), SOL.GetName().c_str());
VectorPreallocateFcbiilu2(solution_data, local_size);
VectorFillFcbiilu2(rhs_data, &RHS[vbeg]);
VectorFinalizeFcbiilu2(rhs_data);
VectorFillFcbiilu2(solution_data, &SOL[vbeg]);
VectorFinalizeFcbiilu2(solution_data);
bool result = SolverSolveFcbiilu2(solver_data, rhs_data, solution_data);
if (result) VectorLoadFcbiilu2(solution_data, &SOL[vbeg]);
return result;
}
bool SolverFCBIILU2::Clear() {
if (matrix_data != NULL) {
MatrixDestroyDataFcbiilu2(&matrix_data);
}
SolverDestroyDataFcbiilu2(&solver_data);
return true;
}
bool SolverFCBIILU2::isMatrixSet() {
return matrix_data != NULL;
}
void SolverFCBIILU2::SetDefaultParameters() {
}
SolverParameter SolverFCBIILU2::GetParameter(std::string name) const {
throw INMOST::SolverUnsupportedOperation;
}
void SolverFCBIILU2::SetParameter(std::string name, std::string value) {
throw INMOST::SolverUnsupportedOperation;
}
const INMOST_DATA_ENUM_TYPE SolverFCBIILU2::Iterations() const {
return static_cast<INMOST_DATA_ENUM_TYPE>(SolverIterationNumberFcbiilu2(solver_data));
}
const INMOST_DATA_REAL_TYPE SolverFCBIILU2::Residual() const {
return SolverResidualNormFcbiilu2(solver_data);
}
const std::string SolverFCBIILU2::ReturnReason() const {
return "Unspecified for FCBIILU2";
}
const std::string SolverFCBIILU2::SolverName() const {
return "fcbiilu2";
}
void SolverFCBIILU2::Finalize() {
SolverFinalizeFcbiilu2();
}
SolverFCBIILU2::~SolverFCBIILU2() {
this->Clear();
}
}
\ No newline at end of file
Source/Solver/solver_fcbiilu2/SolverFCBIILU2.h 0 → 100644
View file @ ed7ce804
#ifndef INMOST_SOLVERFCBIILU2_H
#define INMOST_SOLVERFCBIILU2_H
#include "inmost_solver_interface.h"
#include "solver_fcbiilu2.h"
namespace INMOST {
class SolverFCBIILU2: public SolverInterface {
private:
void *solver_data;
void *matrix_data;
INMOST_DATA_ENUM_TYPE local_size, global_size;
public:
SolverFCBIILU2();
SolverFCBIILU2(const SolverInterface *other);
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 bool Clear();
virtual bool isMatrixSet();
virtual void SetDefaultParameters();
virtual SolverParameter GetParameter(std::string name) const;
virtual void SetParameter(std::string name, std::string value);
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 ~SolverFCBIILU2();
};
}
#endif //INMOST_SOLVERFCBIILU2_H
Source/Solver/solver_fcbiilu2/solver_fcbiilu2.cpp 0 → 100644
View file @ ed7ce804
//#include "inmost_solver.h"
//#if defined(USE_SOLVER)
//#if defined(USE_SOLVER_FCBIILU2)
#include "solver_fcbiilu2.h"
#include <stdio.h>
#include <stdlib.h>
#include <memory.h>
#include <math.h>
#include <string>
static int set_kovl = 0; // number of overlap layers: kovl=0,1,2,...
static double set_tau = 3e-3; // the ILU2 precision (for the submatrix factorization); tau=3e-3
static double set_eps = 1e-5; // the residual precision: ||r|| < eps * ||b||; eps=1e-6
static int set_nit = 999; // number of iterations permitted; nit=999
//static int set_kgmr = 1; // iterative scheme (BiCGStab: kgmr=1; GMR[kdeg]: kgmr>3*kdeg+2)
//static int set_kdeg = 1; // degree of polynomial for GMR[d] (ignored if kgmr=1)
static int set_msglev = 2; // messages level; msglev=0 for silent; msglev=1 to output solution statistics
/* BiCGStab solver structure */
typedef struct
{
int n; // local number of unknowns at the local processor
int nproc; // number of processors
int * ibl; // block splitting: ibl[0]=0; ibl[nproc]=nglob
int * ia; // row pointers: ia[0]=0; ia[nloc]=nzloc
int * ja; // column numbers (NOTE: starting from 0 or 1); ja[nzloc]
double * a; // matrix A coefficients; a[nzloc]
int len_r8; // size of the working memory, set len_r8=nbl for the first call
double * W; // poiter to the working memory W[len_r8]
int kovl; // number of overlap layers: kovl=0,1,2,...
double tau; // the ILU2 precision (for the submatrix factorization); tau=3e-3
double eps; // the residual precision: ||r|| < eps * ||b||; eps=1e-6
int nit; // number of iterations permitted; nit=999
int msglev; // messages level; msglev=0 for silent; msglev=1 to output solution statistics
int ierr; // error flag on return; ierr=0 for success
int istat[16]; // integer statistics array on return
double dstat[16]; // double statistics array on return
double RESID; // residual norm
int ITER; // number of BiCGStab iterations performed
} bcg;
typedef struct
{
int n; // local number of unknowns at the local processor
int nproc; // number of processors
int * ibl; // block splitting: ibl[0]=0; ibl[nproc]=nglob
int * ia;
int * ja;
double * A;
} matrix;
typedef struct
{
int n; // local number of unknowns at the local processor
double * v;
} vector;
/*****************************************************************************/
//----- File: biilu2.h if available
// Solve the linear system A X = B by
// biilu2_bcg solver with working memory allocations and statistics output
int biilu2_bcg (
int *ibl, // block splitting: ibl[0]=0; ibl[nproc]=n
int *ia, // row pointers: ia[0]=0; ia[nloc]=nzloc
int *ja, // column numbers (NOTE: starting from 0 or 1); ja[nzloc]
double *a, // matrix A coefficients; a[nzloc]
double *b, // right-hand side B; b[nloc]
double *x, // initial guess to the solution X on entry; solution X on return; x[nloc]
int job, // job number: 0 - construct preconditioner; 1 - use the previously constructed one
int *len_r8,// size of the working memory, set len_r8=nbl for the first call
double **S, // poiter to the working memory S[len_r8]
int kovl, // number of overlap layers: kovl=0,1,2,...
double tau, // the ILU2 precision (for the submatrix factorization); tau=3e-3
double eps, // the residual precision: ||r|| < eps * ||b||; eps=1e-6
int nit, // number of iterations permitted; nit=999; if(nit==0) preconditioner construction only
int msglev, // messages level; msglev=0 for silent; msglev=1 to output solution statistics
int *ierr, // error flag on return; ierr=0 for success
int *istat, //[16], // integer statistics array on return
double *dstat); //[16]); // double statistics array on return
//
// Here, in notation:
// nproc - the number of processors (or blocks), nproc is equal to the MPI communicator size
// nzloc - the local number of nonzero elements, nzloc=ia[[myid+1]]-ia[ibl[myid]]
// nloc - the local number of unknowns at the current processor, nloc=ibl[myid+1]-ibl[myid]
// n - the total number of unknowns in matrix A, n=ibl[nproc]
//
// ON ENTRY:
// ibl, ia, ja, a, b, x, job, len_r8, S, kovl, tau, eps, nit, msglev
// ON RETURN:
// x, len_r8, S, ierr, istat, dstat
//
// Run selftest by
// mpicc -lm -DBIILU2_BICGST_SELFTEST biilu2.c && mpirun -np 2 a.out
/*****************************************************************************/
/* Initialize bcg solver */
static int initbcg(bcg *s, matrix *A, double eps);
/* Reinitialize solver preconditioner with new matrix A */
static int renewbcg(bcg *s, double *A);
/* Solve linear system */
/*static*/ int solvebcg(bcg *s, vector *b, vector *x);
/* Free memory used by solver */
static void freebcg(bcg *s);
/*****************************************************************************/
/* Initialize solver with new matrix A */
static int newmatrixbcg(bcg *s, matrix *A, bool same_precond)
{
if( s->n != A->n && same_precond ) throw INMOST::CannotReusePreconditionerOfDifferentSize;
s->n = A->n;
s->nproc = A->nproc;
s->ibl = A->ibl;
s->ia = A->ia;
s->ja = A->ja;
if( !same_precond )
{
//do nothing...
//std::cout<<"##### inside newmatrixbcg bef. renewbcg \n";//db!
return renewbcg(s, A->A);
}
else return 0;
}
/* solver */
/* Initialize bcg solver */
int initbcg(bcg *s, matrix *A, double eps)
{
s->eps = eps;
//std::cout<<"##### inside initbcg bef. newmatrixbcg eps="<<eps<<" \n";//db!
return newmatrixbcg(s, A, false);
}
/* Reinitialize solver preconditioner with new matrix A */
int renewbcg(bcg *s, double *A)
{
//reinitialize matrix values
s->a = A;
//BIILU2 preconditioner construction...
s->kovl = set_kovl;
s->tau = set_tau;
s->msglev = set_msglev;
int job = 0;
int maxit = 0;
s->len_r8 = 0; //to be set by nbl inside biilu2_bcg
if(s->W) free(s->W); s->W=NULL;
int ierr = 0;
//std::cout<<"##### inside renewbcg bef. biilu2_bcg\n";//db!
double *B = (double*) malloc(sizeof(double)*s->n); //db!!!!!!!!!!!!!!
//double *X = (double*) malloc(sizeof(double)*s->n); //db!!!!!!!!!!!!!!
biilu2_bcg (s->ibl, s->ia, s->ja, s->a,
//B, X, //!!!!!KAPORIN!!!!!!!!!!!!!
B, NULL, //!!!!!KAPORIN!!!!!!!!!!!!
//NULL, NULL, //!!!!!INK!!!!!!!!!!!
job, &s->len_r8, &s->W,
s->kovl, s->tau, s->eps, maxit, s->msglev,
&ierr, s->istat, s->dstat);
//std::cout<<"##### inside renewbcg aft. biilu2_bcg\n";//db!
free(B); //free(X);//db!!!!!!!!!!!!!!!
if (ierr) printf("initialization of biilu2 failed, ierr=%d\n", ierr);
return ierr;
}
/* Solve linear system */
int solvebcg(bcg *s, vector *b, vector *x)
{
s->kovl = set_kovl;
s->tau = set_tau;
// s->eps = set_eps;
s->nit = set_nit;
s->msglev = set_msglev;
int job = 1;
int maxit = s->nit;
int ierr = 0;
//std::cout<<"##### inside solvebcg bef. biilu2_bcg\n";//db!
biilu2_bcg (s->ibl, s->ia, s->ja, s->a, b->v, x->v,
job, &s->len_r8, &s->W,
s->kovl, s->tau, s->eps, maxit, s->msglev,
&ierr, s->istat, s->dstat);
//std::cout<<"##### inside solvebcg aft. biilu2_bcg\n";//db!
s->ITER = s->istat[2];
s->RESID = s->dstat[2];
return ierr;
}
/* Free memory used by solver */
void freebcg(bcg *s)
{
if(s->W) free(s->W); s->W=NULL;
}
/*****************************************************************************/
void MatrixCopyDataFcbiilu2(void ** ppA, void * pB)
{
matrix * B = (matrix *)pB;
if( ppA == NULL || pB == NULL ) throw INMOST::DataCorruptedInSolver;
*ppA = malloc(sizeof(matrix));
matrix * A = (matrix *)ppA;
A->n = B->n;
if( B->n != 0 )
{
int nnz = B->ia[B->n] - B->ia[0];
A->nproc = B->nproc;
A->ibl = (int *) malloc(sizeof(int)*(A->nproc+1));
memcpy(A->ibl,B->ibl,sizeof(int)*(A->nproc+1));
A->ia = (int *) malloc(sizeof(int)*(A->n+1));
memcpy(A->ia,B->ia,sizeof(int)*(A->n+1));
A->ja = (int *) malloc(sizeof(int)*nnz);
memcpy(A->ja,B->ja,sizeof(int)*nnz);
A->A = (double *) malloc(sizeof(double)*nnz);
memcpy(A->A,B->A,sizeof(double)*nnz);
}
}
void MatrixAssignDataFcbiilu2(void * pA, void * pB)
{
matrix * A = (matrix *)pA;
matrix * B = (matrix *)pB;
if( A == NULL || B == NULL ) throw INMOST::DataCorruptedInSolver;
if( A != B )
{
if( A->n != 0 )
{
free(A->ibl);
free(A->ia);
free(A->ja);
free(A->A);
}
if( B->n != 0 )
{
int nnz = B->ia[B->n] - B->ia[0];
A->n = B->n;
A->nproc = B->nproc;
A->ibl = (int *) malloc(sizeof(int)*(A->nproc+1));
memcpy(A->ibl,B->ibl,sizeof(int)*(A->nproc+1));
A->ia = (int *) malloc(sizeof(int)*(A->n+1));
memcpy(A->ia,B->ia,sizeof(int)*(A->n+1));
A->ja = (int *) malloc(sizeof(int)*nnz);
memcpy(A->ja,B->ja,sizeof(int)*nnz);
A->A = (double *) malloc(sizeof(double)*nnz);
memcpy(A->A,B->A,sizeof(double)*nnz);
}
}
}
void MatrixInitDataFcbiilu2(void ** ppA, INMOST_MPI_Comm comm, const char * name)
{
//std::cout<<"##### ins. MatrixInitDataFcbiilu2 \n";//db!
if( ppA == NULL ) throw INMOST::DataCorruptedInSolver;
if( *ppA == NULL )
{
*ppA = malloc(sizeof(matrix));
matrix * A = (matrix *)*ppA;
A->n = 0;
A->nproc = 0;
//std::cout<<"##### ins. MatrixInitDataFcbiilu2 n=nproc=0 \n";//db!
}
(void) comm;
(void) name;
}
void MatrixDestroyDataFcbiilu2(void ** pA)
{
matrix * A = (matrix *)(*pA);
if( A != NULL )
{
if( A->n != 0 )
{
free(A->ibl);
free(A->ia);
free(A->ja);
free(A->A);
//std::cout<<"##### ins. MatrixDestroyDataFcbiilu2 ...free \n";//db!
}
free(*pA);
*pA = NULL;
}
}
void MatrixFillFcbiilu2(void * pA, int size, int nproc, int * ibl, int * ia, int * ja, double * values)
{
//std::cout<<"##### ins. MatrixFillFcbiilu2 n="<<size<<" nproc="<<nproc<<" \n";//db!
if( pA == NULL ) throw INMOST::DataCorruptedInSolver;
matrix * A = (matrix *) pA;
A->n = size;
A->nproc = nproc;
A->ibl = ibl;
A->ia = ia;
A->ja = ja;
A->A = values;
}
void MatrixFillValuesFcbiilu2(void * pA, double * values)
{
//std::cout<<"##### ins. MatrixFillValuesFcbiilu2 \n";//db!
if( pA == NULL ) throw INMOST::DataCorruptedInSolver;