Skip to content

GitLab

Commit 0af669bb authored by Dmitry Bagaev's avatar Dmitry Bagaev
Browse files

FCBIILU2 interface improvments 

1)SetParameter for fcbiilu2
2) deleted void*
parent 7034199e
Hide whitespace changes
Inline Side-by-side
Source/Solver/solver_fcbiilu2/SolverFCBIILU2.cpp
View file @ 0af669bb
#include "SolverFCBIILU2.h"
#include "solver_fcbiilu2.h"
namespace INMOST {
SolverFCBIILU2::SolverFCBIILU2() {
solver_data = NULL;
matrix_data = NULL;
}
SolverInterface *SolverFCBIILU2::Copy(const SolverInterface *other) {
......@@ -47,7 +49,12 @@ namespace INMOST {
void SolverFCBIILU2::Setup(int *argc, char ***argv, SolverParameters &p) {
SolverInitDataFcbiilu2(&solver_data, communicator, p.solverPrefix.c_str());
SolverInitializeFcbiilu2(argc, argv, p.internalFile.c_str());
solver_data->kovl = 0; // number of overlap layers: kovl=0,1,2,...
solver_data->tau = 3e-3; // the ILU2 precision (for the submatrix factorization); tau=3e-3
solver_data->eps = 1e-5; // the residual precision: ||r|| < eps * ||b||; eps=1e-6
solver_data->nit = 999; // number of iterations permitted; nit=999
solver_data->msglev = 2; // messages level; msglev=0 for silent; msglev=1 to output solution statistics
SolverInitializeFcbiilu2(solver_data, argc, argv, p.internalFile.c_str());
}
void SolverFCBIILU2::SetMatrix(Sparse::Matrix &A, bool ModifiedPattern, bool OldPreconditioner) {
......@@ -117,11 +124,11 @@ namespace INMOST {
INMOST_DATA_ENUM_TYPE vbeg, vend;
RHS.GetInterval(vbeg, vend);
void *rhs_data = NULL;
vector *rhs_data = NULL;
VectorInitDataFcbiilu2(&rhs_data, RHS.GetCommunicator(), RHS.GetName().c_str());
VectorPreallocateFcbiilu2(rhs_data, local_size);
void *solution_data = NULL;
vector *solution_data = NULL;
VectorInitDataFcbiilu2(&solution_data, SOL.GetCommunicator(), SOL.GetName().c_str());
VectorPreallocateFcbiilu2(solution_data, local_size);
VectorFillFcbiilu2(rhs_data, &RHS[vbeg]);
......@@ -153,16 +160,22 @@ namespace INMOST {
}
void SolverFCBIILU2::SetParameter(std::string name, std::string value) {
std::cout << "SolverFCBIILU2::SetParameter unsupported operation" << std::endl;
const char *val = value.c_str();
if (name == "kovl") solver_data->kovl = atoi(val);
else if (name == "tau") solver_data->tau = atof(val);
else if (name == "eps") solver_data->eps = atof(val);
else if (name == "nit") solver_data->nit = atoi(val);
else if (name == "msglev") solver_data->msglev = atoi(val);
else std::cout << "Parameter " << name << " is unknown" << std::endl;
//throw INMOST::SolverUnsupportedOperation;
}
const INMOST_DATA_ENUM_TYPE SolverFCBIILU2::Iterations() const {
return static_cast<INMOST_DATA_ENUM_TYPE>(SolverIterationNumberFcbiilu2(solver_data));
return static_cast<INMOST_DATA_ENUM_TYPE>(solver_data->ITER);
}
const INMOST_DATA_REAL_TYPE SolverFCBIILU2::Residual() const {
return SolverResidualNormFcbiilu2(solver_data);
return solver_data->RESID;
}
const std::string SolverFCBIILU2::ReturnReason() const {
......
Source/Solver/solver_fcbiilu2/SolverFCBIILU2.h
View file @ 0af669bb
......@@ -8,8 +8,8 @@ namespace INMOST {
class SolverFCBIILU2: public SolverInterface {
private:
void *solver_data;
void *matrix_data;
bcg *solver_data;
matrix *matrix_data;
INMOST_DATA_ENUM_TYPE local_size, global_size;
public:
SolverFCBIILU2();
......
Source/Solver/solver_fcbiilu2/solver_fcbiilu2.cpp
View file @ 0af669bb
//#include "inmost_solver.h"
//#if defined(USE_SOLVER)
//#if defined(USE_SOLVER_FCBIILU2)
#include "solver_fcbiilu2.h"
#include <stdio.h>
......@@ -12,59 +8,12 @@
#define T(x) // x // Trace of function calls. Use: "T(x) x" for trace and "T(x)" for silence
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 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]
......@@ -80,7 +29,7 @@ int biilu2_bcg (
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
int *istat, //[16], // integer statistics array on return
double *dstat); //[16]); // double statistics array on return
//
// Here, in notation:
......@@ -101,69 +50,67 @@ int biilu2_bcg (
/* 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;
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...
T(std::cout<<"##### inside newmatrixbcg bef. renewbcg \n";)//db!
return renewbcg(s, A->A);
}
else return 0;
s->ibl = A->ibl;
s->ia = A->ia;
s->ja = A->ja;
if (!same_precond) {
//do nothing...
T(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;
T(std::cout<<"##### inside initbcg bef. newmatrixbcg eps="<<eps<<" \n";)//db!
int initbcg(bcg *s, matrix *A) {
//s->eps = eps;
T(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)
{
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;
//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 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;
T(std::cout<<"##### inside renewbcg bef. biilu2_bcg\n";)//db!
double *B = (double*) malloc(sizeof(double)*s->n); //db!!!!!!!!!!!!!!
T(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);
T(std::cout<<"##### inside renewbcg aft. biilu2_bcg\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);
T(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);
......@@ -171,364 +118,309 @@ int renewbcg(bcg *s, double *A)
}
/* Solve linear system */
int solvebcg(bcg *s, vector *b, vector *x)
{
s->kovl = set_kovl;
s->tau = set_tau;
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;
//s->nit = set_nit;
//s->msglev = set_msglev;
int job = 1;
int maxit = s->nit;
int job = 1;
int maxit = s->nit;
int ierr = 0;
T(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);
T(std::cout<<"##### inside solvebcg aft. biilu2_bcg\n";)//db!
s->ITER = s->istat[2];
T(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);
T(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 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 MatrixCopyDataFcbiilu2(matrix **pA, matrix *B) {
if (pA == NULL || B == NULL) throw INMOST::DataCorruptedInSolver;
*pA = (matrix *) malloc(sizeof(matrix));
matrix *A = *pA;
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 MatrixAssignDataFcbiilu2(matrix *A, matrix *B) {
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)
{
T(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;
T(std::cout<<"##### ins. MatrixInitDataFcbiilu2 n=nproc=0 \n";)//db!
}
void MatrixInitDataFcbiilu2(matrix **ppA, INMOST_MPI_Comm comm, const char *name) {
T(std::cout << "##### ins. MatrixInitDataFcbiilu2 \n";)//db!
if (ppA == NULL) throw INMOST::DataCorruptedInSolver;
if (*ppA == NULL) {
*ppA = (matrix *) malloc(sizeof(matrix));
matrix *A = (matrix *) *ppA;
A->n = 0;
A->nproc = 0;
T(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);
T(std::cout<<"##### ins. MatrixDestroyDataFcbiilu2 ...free \n";)//db!
}
free(*pA);
*pA = NULL;
}
void MatrixDestroyDataFcbiilu2(matrix **pA) {
matrix *A = (matrix *) (*pA);
if (A != NULL) {
if (A->n != 0) {
free(A->ibl);
free(A->ia);
free(A->ja);
free(A->A);
T(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)
{
T(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 MatrixFillFcbiilu2(matrix *pA, int size, int nproc, int *ibl, int *ia, int *ja, double *values) {
T(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)
{
T(std::cout<<"##### ins. MatrixFillValuesFcbiilu2 \n";)//db!
if( pA == NULL ) throw INMOST::DataCorruptedInSolver;
matrix * A = (matrix *) pA;
free(A->A);
A->A = values;
void MatrixFillValuesFcbiilu2(matrix *pA, double *values) {
T(std::cout << "##### ins. MatrixFillValuesFcbiilu2 \n";)//db!
if (pA == NULL) throw INMOST::DataCorruptedInSolver;
matrix *A = (matrix *) pA;
free(A->A);
A->A = values;
}
void MatrixFinalizeFcbiilu2(void * data)
{
//don't need to do anything
void MatrixFinalizeFcbiilu2(matrix *data) {
//don't need to do anything
(void) data;
}
void VectorInitDataFcbiilu2(void ** ppA, INMOST_MPI_Comm comm, const char * name)
{
if( ppA == NULL ) throw INMOST::DataCorruptedInSolver;
*ppA = malloc(sizeof(vector));
vector * A = (vector *)*ppA;
A->n = 0;
void VectorInitDataFcbiilu2(vector **ppA, INMOST_MPI_Comm comm, const char *name) {
if (ppA == NULL) throw INMOST::DataCorruptedInSolver;
*ppA = (vector *) malloc(sizeof(vector));
vector *A = (vector *) *ppA;
A->n = 0;
(void) comm;
(void) name;
}
void VectorCopyDataFcbiilu2(void ** ppA, void * pB)
{
T(std::cout<<"##### ins. VectorCopyDataFcbiilu2 \n";)//db!
if( ppA == NULL || pB == NULL ) throw INMOST::DataCorruptedInSolver;
*ppA = malloc(sizeof(vector));
vector * A = (vector *)*ppA;
vector * B = (vector *)pB;
A->n = B->n;
if( B->n != 0 )
{
A->v = (double *)malloc(sizeof(double)*A->n);
memcpy(A->v,B->v,sizeof(double)*A->n);
}
void VectorCopyDataFcbiilu2(vector **ppA, vector *pB) {
T(std::cout << "##### ins. VectorCopyDataFcbiilu2 \n";)//db!
if (ppA == NULL || pB == NULL) throw INMOST::DataCorruptedInSolver;
*ppA = (vector *) malloc(sizeof(vector));
vector *A = (vector *) *ppA;
vector *B = (vector *) pB;
A->n = B->n;
if (B->n != 0) {
A->v = (double *) malloc(sizeof(double) * A->n);
memcpy(A->v, B->v, sizeof(double) * A->n);
}
}
void VectorAssignDataFcbiilu2(void * pA, void * pB)
{
T(std::cout<<"##### ins. VectorAssignDataFcbiilu2 \n";)//db!
vector * A = (vector *)pA;
vector * B = (vector *)pB;
if( A == NULL || B == NULL ) throw INMOST::DataCorruptedInSolver;
if( A != B )
{
if( A->n != 0 ) free(A->v);
A->n = B->n;
if( B->n != 0 )
{
A->v = (double *) malloc(sizeof(double)*A->n);
memcpy(A->v,B->v,sizeof(double)*A->n);
}
}
void VectorAssignDataFcbiilu2(vector *pA, vector *pB) {
T(std::cout << "##### ins. VectorAssignDataFcbiilu2 \n";)//db!
vector *A = (vector *) pA;
vector *B = (vector *) pB;
if (A ==