Commit b2d17c6f authored by Igor Konshin's avatar Igor Konshin
Browse files

interface to new linear solvers was added

parent 2d42caac
...@@ -60,6 +60,27 @@ else() ...@@ -60,6 +60,27 @@ else()
set(HAVE_SOLVER_MPTILU2 FALSE) set(HAVE_SOLVER_MPTILU2 FALSE)
endif() endif()
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" )
add_definitions(-DHAVE_SOLVER_FCBIILU2)
set(HAVE_SOLVER_FCBIILU2 TRUE)
list(APPEND HEADER solver_fcbiilu2.h)
list(APPEND SOURCE solver_fcbiilu2.cpp 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
"${CMAKE_CURRENT_SOURCE_DIR}/k3d.h" AND EXISTS
"${CMAKE_CURRENT_SOURCE_DIR}/k3d.cpp" )
add_definitions(-DHAVE_SOLVER_K3BIILU2)
set(HAVE_SOLVER_K3BIILU2 TRUE)
list(APPEND HEADER solver_k3biilu2.h k3d.h)
list(APPEND SOURCE solver_k3biilu2.cpp k3d.cpp)
else()
set(HAVE_SOLVER_K3BIILU2 FALSE)
endif()
add_library(inmost STATIC ${SOURCE} ${HEADER}) add_library(inmost STATIC ${SOURCE} ${HEADER})
option(USE_MPI "Compile with MPI support" ON) option(USE_MPI "Compile with MPI support" ON)
......
...@@ -4,5 +4,5 @@ add_subdirectory(DrawMatrix) ...@@ -4,5 +4,5 @@ add_subdirectory(DrawMatrix)
add_subdirectory(MatSolve) add_subdirectory(MatSolve)
add_subdirectory(GridGen) add_subdirectory(GridGen)
add_subdirectory(FVDiscr) add_subdirectory(FVDiscr)
add_subdirectory(OctreeCutcell) #add_subdirectory(OctreeCutcell)
add_subdirectory(Solver) add_subdirectory(Solver)
...@@ -42,15 +42,17 @@ namespace INMOST ...@@ -42,15 +42,17 @@ namespace INMOST
enum Type enum Type
{ {
INNER_ILU2, ///< inner Solver based on BiCGStab(L) solver with second order ILU factorization as preconditioner. INNER_ILU2, ///< inner Solver based on BiCGStab(L) solver with second order ILU factorization as preconditioner.
INNER_DDPQILUC, ///< inner Solver based on BiCGStab(L) solver with second order Crout-ILU with inversed-based condition estimation and unsymmetric reordering for diagonal dominance as preconditioner. INNER_DDPQILUC, ///< inner Solver based on BiCGStab(L) solver with second order Crout-ILU with inversed-based condition estimation and unsymmetric reordering for diagonal dominance as preconditioner.
INNER_MPTILUC, ///< inner Solver based on BiCGStab(L) solver with second order Crout-ILU with inversed-based condition estimation and maximum product transversal reordering as preconditioner. INNER_MPTILUC, ///< inner Solver based on BiCGStab(L) solver with second order Crout-ILU with inversed-based condition estimation and maximum product transversal reordering as preconditioner.
INNER_MPTILU2, ///< inner Solver based on BiCGStab(L) solver with second order ILU and maximum product transversal reordering as preconditioner. INNER_MPTILU2, ///< inner Solver based on BiCGStab(L) solver with second order ILU and maximum product transversal reordering as preconditioner.
Trilinos_Aztec, ///< external Solver AztecOO from Trilinos package Trilinos_Aztec, ///< external Solver AztecOO from Trilinos package.
Trilinos_Belos, ///< external Solver Belos from Trilinos package, currently without preconditioner Trilinos_Belos, ///< external Solver Belos from Trilinos package, currently without preconditioner.
Trilinos_ML, ///< external Solver AztecOO with ML preconditioner Trilinos_ML, ///< external Solver AztecOO with ML preconditioner.
Trilinos_Ifpack,///< external Solver AztecOO with Ifpack preconditioner Trilinos_Ifpack,///< external Solver AztecOO with Ifpack preconditioner.
PETSc, ///< external Solver PETSc, @see http://www.mcs.anl.gov/petsc/. PETSc, ///< external Solver PETSc, @see http://www.mcs.anl.gov/petsc/
ANI ///< external Solver from ANI3D based on ILU2 (sequential Fortran version). ANI, ///< external Solver from ANI3D based on ILU2 (sequential Fortran version), @see http://ani3d.sourceforge.net/
FCBIILU2, // external FCBIILU2 Solver (BIILU2 parallel F2C version).
K3BIILU2 // inner K3BIILU2 Solver (BIILU2 parallel version).
}; };
static INMOST_MPI_Type & GetRowEntryType() {return RowEntryType;} static INMOST_MPI_Type & GetRowEntryType() {return RowEntryType;}
...@@ -464,14 +466,14 @@ namespace INMOST ...@@ -464,14 +466,14 @@ namespace INMOST
INMOST_DATA_ENUM_TYPE Nonzeros; ///< Number of nonzero in linked list. INMOST_DATA_ENUM_TYPE Nonzeros; ///< Number of nonzero in linked list.
interval< INMOST_DATA_ENUM_TYPE, Row::entry > LinkedList; ///< Storage for linked list. interval< INMOST_DATA_ENUM_TYPE, Row::entry > LinkedList; ///< Storage for linked list.
public: public:
/// Default constructor without size specfied. /// Default constructor without size specified.
RowMerger(); RowMerger();
/// Constructor with size specified. /// Constructor with size specified.
/// @param interval_begin First index in linked list. /// @param interval_begin First index in linked list.
/// @param interval_end Last index in linked list. /// @param interval_end Last index in linked list.
/// @param Sorted Result should be sorted. /// @param Sorted Result should be sorted or not.
RowMerger(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool Sorted = true); RowMerger(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool Sorted = true);
/// Constructor that gets sizes from matrix /// Constructor that gets sizes from the matrix.
/// @param A Matrix to get sizes from. /// @param A Matrix to get sizes from.
/// @param Sorted Result should be sorted. /// @param Sorted Result should be sorted.
RowMerger(Matrix & A, bool Sorted = true); RowMerger(Matrix & A, bool Sorted = true);
...@@ -483,14 +485,14 @@ namespace INMOST ...@@ -483,14 +485,14 @@ namespace INMOST
/// This behavior may be changed in future. /// This behavior may be changed in future.
/// @param interval_begin First index in linked list. /// @param interval_begin First index in linked list.
/// @param interval_end Last index in linked list. /// @param interval_end Last index in linked list.
/// @param Sorted Result should be sorted. /// @param Sorted Result should be sorted or not.
void Resize(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool Sorted = true); void Resize(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool Sorted = true);
/// Resize linked list for new matrix. /// Resize linked list for new matrix.
/// \warning /// \warning
/// All contents of linked list will be lost after resize. /// All contents of linked list will be lost after resize.
/// This behavior may be changed in future. /// This behavior may be changed in future.
/// @param A Matrix to get sizes from. /// @param A Matrix to get sizes from.
/// @param Sorted Result should be sorted. /// @param Sorted Result should be sorted or not.
void Resize(Matrix & A, bool Sorted = true); void Resize(Matrix & A, bool Sorted = true);
/// Clear linked list. /// Clear linked list.
void Clear(); void Clear();
...@@ -587,7 +589,7 @@ namespace INMOST ...@@ -587,7 +589,7 @@ namespace INMOST
/// works for: /// works for:
/// INNER_MLILUC /// INNER_MLILUC
/// - "adapt_ddpq_tolerance" - adapt ddpq tolerance depending from the complexity /// - "adapt_ddpq_tolerance" - adapt ddpq tolerance depending from the complexity
// of calculation of Schur complement, /// of calculation of Schur complement,
/// works for: /// works for:
/// INNER_MLILUC /// INNER_MLILUC
void SetParameterEnum(std::string name, INMOST_DATA_ENUM_TYPE value); void SetParameterEnum(std::string name, INMOST_DATA_ENUM_TYPE value);
...@@ -632,12 +634,6 @@ namespace INMOST ...@@ -632,12 +634,6 @@ namespace INMOST
/// Get the package Type. /// Get the package Type.
Type GetPackage() const {return _pack;} Type GetPackage() const {return _pack;}
/// Return the number of iterations performed by the last solution.
/// @see Solver::Solve
INMOST_DATA_ENUM_TYPE Iterations();
/// Return the final residual achieved by the last solution.
/// @see Solver::Solve
INMOST_DATA_REAL_TYPE Residual();
/// Set the matrix and construct the preconditioner. /// Set the matrix and construct the preconditioner.
/// @param A Matrix A in linear problem Ax = b /// @param A Matrix A in linear problem Ax = b
/// @param OldPreconditioner If this parameter is set to true, /// @param OldPreconditioner If this parameter is set to true,
...@@ -662,6 +658,12 @@ namespace INMOST ...@@ -662,6 +658,12 @@ namespace INMOST
/// ///
/// @see Solver::SetMatrix /// @see Solver::SetMatrix
bool Solve(Vector & RHS, Vector & SOL); bool Solve(Vector & RHS, Vector & SOL);
/// Return the number of iterations performed by the last solution.
/// @see Solver::Solve
INMOST_DATA_ENUM_TYPE Iterations();
/// Return the final residual achieved by the last solution.
/// @see Solver::Solve
INMOST_DATA_REAL_TYPE Residual();
/// Get the reason of convergence or divergence of the last solution. /// Get the reason of convergence or divergence of the last solution.
/// @see Solver::Solve /// @see Solver::Solve
std::string GetReason(); std::string GetReason();
......
...@@ -10,6 +10,12 @@ ...@@ -10,6 +10,12 @@
#if defined(HAVE_SOLVER_MPTILU2) #if defined(HAVE_SOLVER_MPTILU2)
#include "solver_mtilu2.hpp" #include "solver_mtilu2.hpp"
#endif #endif
#if defined(HAVE_SOLVER_FCBIILU2)
#include "solver_fcbiilu2.h"
#endif
#if defined(HAVE_SOLVER_K3BIILU2)
#include "solver_k3biilu2.h"
#endif
#include "solver_bcgsl.hpp" #include "solver_bcgsl.hpp"
#include <fstream> #include <fstream>
#include <sstream> #include <sstream>
...@@ -52,6 +58,8 @@ namespace INMOST ...@@ -52,6 +58,8 @@ namespace INMOST
static std::string trilinos_belos_database_file = ""; static std::string trilinos_belos_database_file = "";
static std::string trilinos_ml_database_file = ""; static std::string trilinos_ml_database_file = "";
static std::string ani_database_file = ""; static std::string ani_database_file = "";
static std::string fcbiilu2_database_file = "";
static std::string k3biilu2_database_file = "";
#if defined(USE_MPI) #if defined(USE_MPI)
...@@ -1319,9 +1327,9 @@ namespace INMOST ...@@ -1319,9 +1327,9 @@ namespace INMOST
#endif #endif
} }
//////////////////////////////////////////////////////////////////////// //======================================================================
//SOLVER SOURCE CODE //SOLVER SOURCE CODE
//////////////////////////////////////////////////////////////////////// //======================================================================
void Solver::Initialize(int * argc, char *** argv, const char * database) void Solver::Initialize(int * argc, char *** argv, const char * database)
{ {
(void)database; (void)database;
...@@ -1355,6 +1363,10 @@ namespace INMOST ...@@ -1355,6 +1363,10 @@ namespace INMOST
trilinos_belos_database_file = std::string(str+l); trilinos_belos_database_file = std::string(str+l);
else if( !strncmp(str,"ani",k) ) else if( !strncmp(str,"ani",k) )
ani_database_file = std::string(str+l); ani_database_file = std::string(str+l);
else if( !strncmp(str,"fcbiilu2",k) )
fcbiilu2_database_file = std::string(str+l);
else if( !strncmp(str,"k3biilu2",k) )
k3biilu2_database_file = std::string(str+l);
} }
} }
//std::cout << "PETSc \"" << petsc_database_file << "\"" << std::endl; //std::cout << "PETSc \"" << petsc_database_file << "\"" << std::endl;
...@@ -1365,6 +1377,12 @@ namespace INMOST ...@@ -1365,6 +1377,12 @@ namespace INMOST
#if defined(USE_SOLVER_ANI) #if defined(USE_SOLVER_ANI)
SolverInitializeAni(argc,argv,ani_database_file.c_str()); SolverInitializeAni(argc,argv,ani_database_file.c_str());
#endif #endif
#if defined(HAVE_SOLVER_FCBIILU2)
SolverInitializeFcbiilu2(argc,argv,fcbiilu2_database_file.c_str());
#endif
#if defined(HAVE_SOLVER_K3BIILU2)
SolverInitializeK3biilu2(argc,argv,k3biilu2_database_file.c_str());
#endif
#if defined(USE_MPI) #if defined(USE_MPI)
{ {
int flag = 0; int flag = 0;
...@@ -1410,6 +1428,12 @@ namespace INMOST ...@@ -1410,6 +1428,12 @@ namespace INMOST
#if defined(USE_SOLVER_ANI) #if defined(USE_SOLVER_ANI)
SolverFinalizeAni(); SolverFinalizeAni();
#endif #endif
#if defined(HAVE_SOLVER_FCBIILU2)
SolverFinalizeFcbiilu2();
#endif
#if defined(HAVE_SOLVER_K3BIILU2)
SolverFinalizeK3biilu2();
#endif
#if defined(USE_MPI) #if defined(USE_MPI)
{ {
int flag = 0; int flag = 0;
...@@ -1526,6 +1550,18 @@ namespace INMOST ...@@ -1526,6 +1550,18 @@ namespace INMOST
SolverInitDataAni(&solver_data,_comm,name.c_str()); SolverInitDataAni(&solver_data,_comm,name.c_str());
} }
#endif #endif
#if defined(HAVE_SOLVER_FCBIILU2)
if( _pack == FCBIILU2 )
{
SolverInitDataFcbiilu2(&solver_data,_comm,name.c_str());
}
#endif
#if defined(HAVE_SOLVER_K3BIILU2)
if( _pack == K3BIILU2 )
{
SolverInitDataK3biilu2(&solver_data,_comm,name.c_str());
}
#endif
#if defined(USE_SOLVER_TRILINOS) #if defined(USE_SOLVER_TRILINOS)
if( _pack == Trilinos_Aztec || if( _pack == Trilinos_Aztec ||
_pack == Trilinos_Belos || _pack == Trilinos_Belos ||
...@@ -1613,6 +1649,36 @@ namespace INMOST ...@@ -1613,6 +1649,36 @@ namespace INMOST
VectorCopyDataAni(&solution_data,other.solution_data); VectorCopyDataAni(&solution_data,other.solution_data);
} }
#endif #endif
#if defined(HAVE_SOLVER_FCBIILU2)
if( _pack == FCBIILU2 )
{
SolverCopyDataFcbiilu2(&solver_data,other.solver_data,comm);
if( other.matrix_data != NULL )
{
MatrixCopyDataFcbiilu2(&matrix_data,other.matrix_data);
SolverSetMatrixFcbiilu2(solver_data,matrix_data,false,false);
}
if( other.rhs_data != NULL )
VectorCopyDataFcbiilu2(&rhs_data,other.rhs_data);
if( other.solution_data != NULL )
VectorCopyDataFcbiilu2(&solution_data,other.solution_data);
}
#endif
#if defined(HAVE_SOLVER_K3BIILU2)
if( _pack == K3BIILU2 )
{
SolverCopyDataK3biilu2(&solver_data,other.solver_data,comm);
if( other.matrix_data != NULL )
{
MatrixCopyDataK3biilu2(&matrix_data,other.matrix_data);
SolverSetMatrixK3biilu2(solver_data,matrix_data,false,false);
}
if( other.rhs_data != NULL )
VectorCopyDataK3biilu2(&rhs_data,other.rhs_data);
if( other.solution_data != NULL )
VectorCopyDataK3biilu2(&solution_data,other.solution_data);
}
#endif
#if defined(USE_SOLVER_TRILINOS) #if defined(USE_SOLVER_TRILINOS)
if( _pack == Trilinos_Aztec || if( _pack == Trilinos_Aztec ||
_pack == Trilinos_Belos || _pack == Trilinos_Belos ||
...@@ -1655,6 +1721,30 @@ namespace INMOST ...@@ -1655,6 +1721,30 @@ namespace INMOST
SolverDestroyDataAni(&solver_data); SolverDestroyDataAni(&solver_data);
} }
#endif #endif
#if defined(HAVE_SOLVER_FCBIILU2)
if( _pack == FCBIILU2 )
{
if( matrix_data != NULL )
MatrixDestroyDataFcbiilu2(&matrix_data);
if( rhs_data != NULL )
VectorDestroyDataFcbiilu2(&rhs_data);
if( solution_data != NULL )
VectorDestroyDataFcbiilu2(&solution_data);
SolverDestroyDataFcbiilu2(&solver_data);
}
#endif
#if defined(HAVE_SOLVER_K3BIILU2)
if( _pack == K3BIILU2 )
{
if( matrix_data != NULL )
MatrixDestroyDataK3biilu2(&matrix_data);
if( rhs_data != NULL )
VectorDestroyDataK3biilu2(&rhs_data);
if( solution_data != NULL )
VectorDestroyDataK3biilu2(&solution_data);
SolverDestroyDataK3biilu2(&solver_data);
}
#endif
#if defined(USE_SOLVER_TRILINOS) #if defined(USE_SOLVER_TRILINOS)
if( _pack == Trilinos_Aztec || if( _pack == Trilinos_Aztec ||
_pack == Trilinos_Belos || _pack == Trilinos_Belos ||
...@@ -1752,6 +1842,62 @@ namespace INMOST ...@@ -1752,6 +1842,62 @@ namespace INMOST
} }
} }
#endif #endif
#if defined(HAVE_SOLVER_FCBIILU2)
if( _pack == FCBIILU2 )
{
SolverAssignDataFcbiilu2(solver_data,other.solver_data);
if( other.matrix_data != NULL )
{
if( matrix_data != NULL )
MatrixAssignDataFcbiilu2(matrix_data,other.matrix_data);
else
MatrixCopyDataFcbiilu2(&matrix_data,other.matrix_data);
SolverSetMatrixFcbiilu2(solver_data,matrix_data,false,false);
}
if( other.rhs_data != NULL )
{
if( rhs_data != NULL )
VectorAssignDataFcbiilu2(rhs_data,other.rhs_data);
else
VectorCopyDataFcbiilu2(&rhs_data,other.rhs_data);
}
if( other.solution_data != NULL )
{
if( solution_data != NULL )
VectorAssignDataFcbiilu2(solution_data,other.solution_data);
else
VectorCopyDataFcbiilu2(&solution_data,other.solution_data);
}
}
#endif
#if defined(HAVE_SOLVER_K3BIILU2)
if( _pack == K3BIILU2 )
{
SolverAssignDataK3biilu2(solver_data,other.solver_data);
if( other.matrix_data != NULL )
{
if( matrix_data != NULL )
MatrixAssignDataK3biilu2(matrix_data,other.matrix_data);
else
MatrixCopyDataK3biilu2(&matrix_data,other.matrix_data);
SolverSetMatrixK3biilu2(solver_data,matrix_data,false,false);
}
if( other.rhs_data != NULL )
{
if( rhs_data != NULL )
VectorAssignDataK3biilu2(rhs_data,other.rhs_data);
else
VectorCopyDataK3biilu2(&rhs_data,other.rhs_data);
}
if( other.solution_data != NULL )
{
if( solution_data != NULL )
VectorAssignDataK3biilu2(solution_data,other.solution_data);
else
VectorCopyDataK3biilu2(&solution_data,other.solution_data);
}
}
#endif
#if defined(USE_SOLVER_TRILINOS) #if defined(USE_SOLVER_TRILINOS)
if( _pack == Trilinos_Aztec || if( _pack == Trilinos_Aztec ||
_pack == Trilinos_Belos || _pack == Trilinos_Belos ||
...@@ -1772,7 +1918,14 @@ namespace INMOST ...@@ -1772,7 +1918,14 @@ namespace INMOST
void Solver::SetMatrix(Matrix & A, bool OldPreconditioner) void Solver::SetMatrix(Matrix & A, bool OldPreconditioner)
{ {
(void) OldPreconditioner; (void) OldPreconditioner;
std::cout<<"##### SetMatrix: bef.: the very-very beginning... _pack="<<_pack<<std::endl;//db!
bool ok = false; bool ok = false;
#if defined(HAVE_SOLVER_FCBIILU2)
std::cout<<"##### SetMatrix: HAVE_SOLVER_FCBIILU2 is set, OK!!!!!!!!!!!!! _pack="<<_pack<<std::endl;//db!
#endif
#if defined(HAVE_SOLVER_K3BIILU2)
std::cout<<"##### SetMatrix: HAVE_SOLVER_K3BIILU2 is set, OK!!!!!!!!!!!!! _pack="<<_pack<<std::endl;//db!
#endif
#if defined(USE_SOLVER_PETSC) #if defined(USE_SOLVER_PETSC)
if( _pack == PETSc ) if( _pack == PETSc )
{ {
...@@ -1927,6 +2080,164 @@ namespace INMOST ...@@ -1927,6 +2080,164 @@ namespace INMOST
ok = true; ok = true;
} }
#endif #endif
#if defined(HAVE_SOLVER_FCBIILU2)
if( _pack == FCBIILU2 )
{
std::cout<<"##### bef. FCBIILU2: bool modified_pattern = ... "<<std::endl;//db!
bool modified_pattern = false;
for(Matrix::iterator it = A.Begin(); it != A.End() && !modified_pattern; ++it)
modified_pattern |= it->modified_pattern;
//~ if( A.comm != comm ) throw DifferentCommunicatorInSolver;
if( matrix_data == NULL )
{
MatrixInitDataFcbiilu2(&matrix_data,A.GetCommunicator(),A.GetName().c_str());
modified_pattern = true;
}
std::cout<<"##### bef. if( modified_pattern )... "<<modified_pattern<<std::endl;//db!
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));
std::cout<<"##### bef. for(...) myid="<<myid<<" A.Size()="<<A.Size()<<" n="<<n<<" nproc="<<nproc<<" ibl: "<<ibl[0]<<" "<<ibl[1]<<" ... "<<ibl[nproc]<<std::endl;//db!
for(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(Matrix::iterator it = A.Begin(); it != A.End(); it++)
{
for(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;
}
std::cout<<"##### bef. MatrixFillFcbiilu2 myid="<<myid<<" n="<<n<<" nnz="<<nnz<<" shift="<<shift<<" ia: "<<ia[0]<<" "<<ia[1]<<" "<<ia[2]<<" ... "<<ia[n]<<" ja: "<<ja[0]<<" "<<ja[1]<<" "<<ja[2]<<" ... "<<ja[nnz-1]<<" a: "<<values[0]<<" "<<values[1]<<" "<<values[2]<<" ... "<<values[nnz-1]<<std::endl;//db!
MatrixFillFcbiilu2(matrix_data,n,nproc,ibl,ia,ja,values);
}
else
{
INMOST_DATA_ENUM_TYPE nnz = 0, k = 0;
for(Matrix::iterator it = A.Begin(); it != A.End(); it++) nnz += it->Size();
double * values = (double *)malloc(sizeof(double)*nnz);
k = 0;
for(Matrix::iterator it = A.Begin(); it != A.End(); it++)
for(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);
std::cout<<"##### bef. ok=true"<<std::endl;//db!
ok = true;
}
#endif
#if defined(HAVE_SOLVER_K3BIILU2)
if( _pack == K3BIILU2 )
{
std::cout<<"##### bef. K3BIILU2: bool modified_pattern = ... "<<std::endl;//db!
bool modified_pattern = false;
for(Matrix::iterator it = A.Begin(); it != A.End() && !modified_pattern; ++it)
modified_pattern |= it->modified_pattern;
//~ if( A.comm != comm ) throw DifferentCommunicatorInSolver;
if( matrix_data == NULL )
{
MatrixInitDataK3biilu2(&matrix_data,A.GetCommunicator(),A.GetName().c_str());
modified_pattern = true;
}
std::cout<<"##### bef. if( modified_pattern )... "<<modified_pattern<<std::endl;//db!
if( modified_pattern )
{
global_size = local_size = A.Size();
MatrixDestroyDataK3biilu2(&matrix_data);
MatrixInitDataK3biilu2(&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));
std::cout<<"##### bef. for(...) myid="<<myid<<" A.Size()="<<A.Size()<<" n="<<n<<" nproc="<<nproc<<" ibl: "<<ibl[0]<<" "<<ibl[1]<<" ... "<<ibl[nproc]<<std::endl;//db!
for(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(Matrix::iterator it = A.Begin(); it != A.End(); it++)
{
for(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;
}
std::cout<<"##### bef. MatrixFillK3biilu2 myid="<<myid<<" n="<<n<<" nnz="<<nnz<<" shift="<<shift<<" ia: "<<ia[0]<<" "<<ia[1]<<" "<<ia[2]<<" ... "<<ia[n]<<" ja: "<<ja[0]<<" "<<ja[1]<<" "<<ja[2]<<" ... "<<ja[nnz-1]<<" a: "<<values[0]<<" "<<values[1]<<" "<<values[2]<<" ... "<<values[nnz-1]<<std::endl;//db!
MatrixFillK3biilu2(matrix_data,n,nproc,ibl,ia,ja,values);
}
else
{
INMOST_DATA_ENUM_TYPE nnz = 0, k = 0;
for(Matrix::iterator it = A.Begin(); it != A.End(); it++) nnz += it->Size();
double * values = (double *)malloc(sizeof(double)*nnz);
k = 0;
for(Matrix::iterator it = A.Begin(); it != A.End(); it++)
for(Row::iterator jt = it->Begin(); jt != it->End(); jt++)
values[k++] = jt->second;
MatrixFillValuesK3biilu2(matrix_data,values);
}
MatrixFinalizeK3biilu2(matrix_data);
SolverSetMatrixK3biilu2(solver_data,matrix_data,modified_pattern,OldPreconditioner);