Skip to content

GitLab

Commit 4abb074a authored by Igor Konshin's avatar Igor Konshin
Browse files

Merge pull request #18 from bvdmitri/master 

mat solve example update, k3biilu2 and fcbiilu2 interfaces
parents 0002d429 1a4bb4bd
Expand all Hide whitespace changes
Inline Side-by-side
.gitignore
View file @ 4abb074a
.svn
Tests/solver_test001/matrices
SyncToy*
build*/
Examples/CMakeLists.txt
View file @ 4abb074a
......@@ -7,4 +7,4 @@ add_subdirectory(FVDiscr)
add_subdirectory(ADFVDiscr)
add_subdirectory(ADMFD)
#add_subdirectory(OctreeCutcell)
add_subdirectory(Solver)
add_subdirectory(Solver)
\ No newline at end of file
Examples/MatSolve/CMakeLists.txt
View file @ 4abb074a
project(MatSolve)
set(SOURCE main.cpp)
set(SOURCE main.cpp inner_parser.cpp)
add_executable(MatSolve ${SOURCE})
......@@ -7,10 +7,10 @@ target_link_libraries(MatSolve inmost)
if(USE_MPI)
message("linking MatSolve with MPI")
target_link_libraries(MatSolve ${MPI_LIBRARIES})
target_link_libraries(MatSolve ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
set_target_properties(MatSolve PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif()
endif(USE_MPI)
if(USE_SOLVER)
......
Examples/MatSolve/inner_options.txt 0 → 100644
View file @ 4abb074a
enum maximum_iterations 300
enum gmres_substeps 4
real relative_tolerance 1.0e-5
real absolute_tolerance 1.0e-10
real divergence_tolerance 1e+200
enum reorder_nonzeros 0
enum rescale_iterations 8
enum adapt_ddpq_tolerance 0
real drop_tolerance 3.0e-5
real reuse_tolerance 1.0e-5
real ddpq_tolerance 0.0
enum condition_estimation 1
enum schwartz_overlap 6
Examples/MatSolve/inner_parser.cpp 0 → 100644
View file @ 4abb074a
#include "inner_parser.h"
void removeSpaces(char* source) {
char* i = source;
char* j = source;
while(*j != 0)
{
*i = *j++;
if(*i != ' ')
i++;
}
*i = 0;
}
InnerOptions *parseInnerDatabaseOptions(char *optionsFile) {
FILE *databaseFile = fopen(optionsFile, "r");
if (!databaseFile) {
std::cout << "Inner options file not found" << std::endl;
return nullptr;
}
InnerOptions *options = new InnerOptions();
char *tmp = (char *) calloc(256, sizeof(char));
char *parameterName = (char *) calloc(128, sizeof(char));
char *parameterValue = (char *) calloc(128, sizeof(char));
char *type = (char *) calloc(36, sizeof(char));
while (!feof(databaseFile) && fgets(tmp, 256, databaseFile)) {
//removeSpaces(tmp);
char *line = tmp;
//Comment str
if (line[0] == '#') continue;
//First 4 chars is 'real' or 'enum'
bool isReal = false, isEnum = false;
sscanf(line, "%s %s %s", type, parameterName, parameterValue);
if (strncmp(type, "real", 4) == 0) {
//fprintf(stdout, "Real parameter:");
isReal = true;
} else if (strncmp(type, "enum", 4) == 0) {
//fprintf(stdout, "Enum parameter:");
isEnum = true;
} else {
fprintf(stderr, "Skipping bad line: %s", line);
continue;
}
options->options.push_back(new InnerOption(std::string(parameterName), std::string(parameterValue), isReal ? REAL : ENUM));
}
free(parameterValue);
free(parameterName);
free(tmp);
free(type);
return options;
}
char *findInnerOptions(const char *databaseFilePath) {
FILE *databaseFile = fopen(databaseFilePath, "r");
if (databaseFile == NULL) {
fprintf(stderr, "Database file not found\n");
#if defined(USE_MPI)
MPI_Finalize();
#endif
exit(1);
}
char *tmp = (char *) calloc(256, sizeof(char));
char *parameterName = (char *) calloc(64, sizeof(char));
char *parameterValue = (char *) calloc(64, sizeof(char));
char *fileName = NULL;
while (!feof(databaseFile) && fgets(tmp, 256, databaseFile)) {
removeSpaces(tmp);
char *line = tmp;
if (strncmp(line, "inner", 5) != 0) continue;
line += 6;
size_t fileNameLength = 0;
while (!isspace(*(line + fileNameLength))) fileNameLength += 1;
fileName = (char *) calloc(fileNameLength + 1, sizeof(char));
memcpy(fileName, line, fileNameLength);
fileName[fileNameLength] = 0;
}
free(tmp);
return fileName;
}
Examples/MatSolve/inner_parser.h 0 → 100644
View file @ 4abb074a
#ifndef INMOST_INNER_PARSER_H_H
#define INMOST_INNER_PARSER_H_H
#include <string>
#include <iostream>
#include <vector>
enum OptionType {
REAL,
ENUM
};
class InnerOption {
public:
InnerOption(const std::string &name, const std::string &value, OptionType type) : name(name), value(value),
type(type) { }
std::string name;
std::string value;
OptionType type;
};
class InnerOptions {
public:
std::vector<InnerOption *> options;
};
char *findInnerOptions(const char *databaseFilePath);
InnerOptions *parseInnerDatabaseOptions(char *optionsFile);
#endif //INMOST_INNER_PARSER_H_H
Examples/MatSolve/main.cpp
View file @ 4abb074a
This diff is collapsed.
Source/Solver/CMakeLists.txt
View file @ 4abb074a
......@@ -42,23 +42,25 @@ if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/solver_fcbiilu2.h" AND EXISTS "${CMAKE_C
set(SOLVER_DEFINITIONS ${SOLVER_DEFINITIONS} -DHAVE_SOLVER_FCBIILU2)
set(HAVE_SOLVER_FCBIILU2 TRUE PARENT_SCOPE)
list(APPEND HEADER ${CMAKE_CURRENT_SOURCE_DIR}/solver_fcbiilu2.h)
list(APPEND SOURCE ${CMAKE_CURRENT_SOURCE_DIR}/solver_fcbiilu2.cpp fcbiilu2.cpp)
list(APPEND SOURCE ${CMAKE_CURRENT_SOURCE_DIR}/solver_fcbiilu2.cpp ${CMAKE_CURRENT_SOURCE_DIR}/fcbiilu2.cpp)
else()
set(HAVE_SOLVER_FCBIILU2 FALSE)
endif()
if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/solver_k3biilu2.h" AND EXISTS
"${CMAKE_CURRENT_SOURCE_DIR}/solver_k3biilu2.cpp" AND EXISTS
if (USE_MPI)
if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/solver_k3biilu2.h" AND EXISTS
"${CMAKE_CURRENT_SOURCE_DIR}/solver_k3biilu2.cpp" AND EXISTS
"${CMAKE_CURRENT_SOURCE_DIR}/k3d.h" AND EXISTS
"${CMAKE_CURRENT_SOURCE_DIR}/k3d.cpp" )
set(SOLVER_DEFINITIONS ${SOLVER_DEFINITIONS} -DHAVE_SOLVER_K3BIILU2)
set(HAVE_SOLVER_K3BIILU2 TRUE PARENT_SCOPE)
list(APPEND HEADER ${CMAKE_CURRENT_SOURCE_DIR}/solver_k3biilu2.h ${CMAKE_CURRENT_SOURCE_DIR}/k3d.h)
list(APPEND SOURCE ${CMAKE_CURRENT_SOURCE_DIR}/solver_k3biilu2.cpp ${CMAKE_CURRENT_SOURCE_DIR}/k3d.cpp)
else()
set(HAVE_SOLVER_K3BIILU2 FALSE PARENT_SCOPE)
set(SOLVER_DEFINITIONS ${SOLVER_DEFINITIONS} -DHAVE_SOLVER_K3BIILU2)
set(HAVE_SOLVER_K3BIILU2 TRUE PARENT_SCOPE)
list(APPEND HEADER ${CMAKE_CURRENT_SOURCE_DIR}/solver_k3biilu2.h ${CMAKE_CURRENT_SOURCE_DIR}/k3d.h)
list(APPEND SOURCE ${CMAKE_CURRENT_SOURCE_DIR}/solver_k3biilu2.cpp ${CMAKE_CURRENT_SOURCE_DIR}/k3d.cpp)
else()
set(HAVE_SOLVER_K3BIILU2 FALSE PARENT_SCOPE)
endif()
endif()
set(HEADER ${HEADER} PARENT_SCOPE)
set(SOURCE ${SOURCE} PARENT_SCOPE)
set(SOLVER_DEFINITIONS ${SOLVER_DEFINITIONS} PARENT_SCOPE)
\ No newline at end of file
Source/Solver/solver.cpp
View file @ 4abb074a
......@@ -1644,7 +1644,7 @@ namespace INMOST
ia[0] = shift;
for(Sparse::Matrix::iterator it = A.Begin(); it != A.End(); it++)
{
for(Row::iterator jt = it->Begin(); jt != it->End(); jt++)
for(Sparse::Row::iterator jt = it->Begin(); jt != it->End(); jt++)
{
ja[k] = jt->first + 1;
values[k] = jt->second;
......
Source/Solver/solver_fcbiilu2.cpp 0 → 100644
View file @ 4abb074a
//#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;
matrix * A = (matrix *) pA;
free(A->A);
A->A = values;
}
void MatrixFinalizeFcbiilu2(void * 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) comm;
(void) name;
}
void VectorCopyDataFcbiilu2(void ** ppA, void * pB)
{
//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 VectorAssignDataFcbiilu2(void * pA, void * pB)
{
//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 VectorPreallocateFcbiilu2(void * pA, int size)
{
vector * A = (vector *)pA;
if( A == NULL ) throw INMOST::DataCorruptedInSolver;
A->n = size;
A->v = (double *)malloc(sizeof(double)*size);
}
void VectorFillFcbiilu2(void * pA, double * values)
{
vector * A = (vector *)pA;
if( A == NULL ) throw INMOST::DataCorruptedInSolver;
memcpy(A->v,values,sizeof(double)*A->n);
}
void VectorLoadFcbiilu2(void * pA, double * values)
{
vector * A = (vector *)pA;
if( A == NULL ) throw INMOST::DataCorruptedInSolver;
memcpy(values,A->v,sizeof(double)*A->n);
}