Commit 4a1aff21 authored by Dmitry Bagaev's avatar Dmitry Bagaev

Inner solvers refactored using SolverInterface, Trillions on the way

parent efb0bb4a
......@@ -6,6 +6,7 @@
#include "inmost.h"
#include "inner_parser.h"
using namespace INMOST;
#if defined(USE_MPI)
......@@ -14,14 +15,14 @@ using namespace INMOST;
#define BARRIER
#endif
int main(int argc, char ** argv) {
int main(int argc, char **argv) {
int processRank = 0, processorsCount = 1;
#if defined(USE_MPI)
#if defined(USE_MPI)
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &processRank); // Get the rank of the current process
MPI_Comm_size(MPI_COMM_WORLD, &processorsCount ); // Get the total number of processors used
#endif
MPI_Comm_size(MPI_COMM_WORLD, &processorsCount); // Get the total number of processors used
#endif
if (processRank == 0) {
std::cout << "Starting MatSolve2" << std::endl;
......@@ -57,10 +58,14 @@ int main(int argc, char ** argv) {
std::cout << "-x, --xvector <X vector file name>" << std::endl;
std::cout << "-p, --parameters <Solver parameters file name>" << std::endl;
std::cout << "-t, --type <Solver type name>" << std::endl;
std::cout << " Available solvers:" << std::endl;
Solver::Initialize(NULL, NULL, NULL);
std::vector<std::string> availableSolvers = SolverFactory::getAvailableSolvers();
for (auto it = availableSolvers.begin(); it != availableSolvers.end(); it++) {
std::cout << " " << *it << std::endl;
}
Solver::Finalize();
}
#if defined(USE_MPI)
MPI_Finalize();
#endif
return 0;
}
//Matrix file name found with -m or --matrix options
......@@ -127,7 +132,7 @@ int main(int argc, char ** argv) {
if (!matrixFound) {
if (processRank == 0) {
std::cout <<
"Matrix not found, you can specify matrix file name using -m or --matrix options, otherwise specify -h option to see all options, exiting...";
"Matrix not found, you can specify matrix file name using -m or --matrix options, otherwise specify -h option to see all options, exiting...";
}
return -1;
}
......@@ -135,23 +140,33 @@ int main(int argc, char ** argv) {
if (!typeFound) {
if (processRank == 0) {
std::cout <<
"Solver type not found in command line, you can specify solver type with -t or --type option, using INNER_ILU2 solver by default." <<
std::endl;
"Solver type not found in command line, you can specify solver type with -t or --type option, using INNER_ILU2 solver by default."
<<
std::endl;
}
}
if (!vectorBFound) {
if (processRank == 0) {
std::cout <<
"B vector not found, you can specify b vector file name with -b or --bvector option, using identity vector by default." <<
std::endl;
"B vector not found, you can specify b vector file name with -b or --bvector option, using identity vector by default."
<<
std::endl;
}
}
// Initialize the linear solver in accordance with args
Solver::Initialize(&argc, &argv, parametersFound ? parametersFileName.c_str() : NULL);
if (!SolverFactory::isSolverAvailable(solverName)) {
if (processRank == 0) {
std::cout << "Solver " << solverName << " is not available" << std::endl;
}
Solver::Finalize();
exit(1);
}
Solver solver = Solver(solverName);
//solver2.Finalize();
if (processRank == 0) {
std::cout << "Solving with " << solverName << std::endl;
......@@ -195,9 +210,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()));
solver.SetPropertyEnum(option->name, (unsigned int) atoi(option->value.c_str()));
} else {
//s.SetParameterReal(option->name, atof(option->value.c_str()));
solver.SetPropertyReal(option->name, atof(option->value.c_str()));
};
}
delete options;
......@@ -207,24 +222,6 @@ int main(int argc, char ** argv) {
}
BARRIER
//s.SetParameterEnum("maximum_iterations", 1000);
//s.SetParameterEnum("gmres_substeps", 4);
//s.SetParameterReal("relative_tolerance", 1.0e-6);
//s.SetParameterReal("absolute_tolerance", 1.0e-16);
//s.SetParameterReal("divergence_tolerance", 1e+200);
//s.SetParameterEnum("reorder_nonzeros", 0);
//s.SetParameterEnum("rescale_iterations", 8);
//s.SetParameterEnum("adapt_ddpq_tolerance", 0);
//s.SetParameterReal("drop_tolerance", 3.0e-3);
//s.SetParameterReal("reuse_tolerance", 1.0e-5);
//s.SetParameterReal("ddpq_tolerance", 0.0);
//s.SetParameterEnum("condition_estimation", 1);
//s.SetParameterEnum("schwartz_overlap", 3);
tempTimer = Timer();
solver.SetMatrix(mat);
//s.SetMatrix(mat); // Compute the preconditioner for the original matrix
......@@ -266,7 +263,7 @@ int main(int argc, char ** argv) {
info.RestoreVector(x);
if (processRank == 0) {
std::cout << "||Ax-b||=" << sqrt(recv[0]) << " ||b||=" << sqrt(recv[1]) << " ||Ax-b||/||b||=" <<
sqrt(recv[0] / (recv[1] + 1.0e-100)) << std::endl;
sqrt(recv[0] / (recv[1] + 1.0e-100)) << std::endl;
std::cout << "norms: " << Timer() - tempTimer << std::endl;
std::cout << processorsCount << " processors for Solver " << solverName;
if (success) {
......@@ -312,15 +309,15 @@ int main(int argc, char ** argv) {
norm += 1.0e-100;
if (processRank == 0) {
std::cout << "Difference with exact solution \"" << vectorXFileName << "\": " << std::scientific <<
std::setprecision(6) << std::endl;
std::setprecision(6) << std::endl;
std::cout << "dif1 = " << dif1 << " dif2 = " << dif2 << " dif8 = " << dif8 << " ||ex||_1 = " <<
norm << std::endl;
norm << std::endl;
std::cout << "rel1 = " << dif1 / norm << " rel2 = " << dif2 / norm << " rel8 = " << dif8 / norm <<
std::endl;
std::endl;
}
}
}
BARRIER
Solver::Finalize(); // Finalize solver and close MPI activity
return 0;
Solver::Finalize(); // Finalize solver and close MPI activity
return 0;
}
set(SOURCE
${SOURCE}
PARENT_SCOPE
)
${SOURCE}
PARENT_SCOPE
)
set(HEADER
${HEADER}
${CMAKE_CURRENT_SOURCE_DIR}/inmost.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_options_cmake.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_common.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_data.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_dense.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_mesh.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_solver.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_solver_interface.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_partitioner.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_autodiff.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_expression.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_variable.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_sparse.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_nonlinear.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_xml.h
${CMAKE_CURRENT_SOURCE_DIR}/container.hpp
PARENT_SCOPE
)
\ No newline at end of file
${HEADER}
${CMAKE_CURRENT_SOURCE_DIR}/inmost.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_options_cmake.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_common.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_data.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_dense.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_mesh.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_solver.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_solver_interface.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_solver_factory
${CMAKE_CURRENT_SOURCE_DIR}/inmost_partitioner.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_autodiff.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_expression.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_variable.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_sparse.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_nonlinear.h
${CMAKE_CURRENT_SOURCE_DIR}/inmost_xml.h
${CMAKE_CURRENT_SOURCE_DIR}/container.hpp
PARENT_SCOPE
)
\ No newline at end of file
......@@ -7,6 +7,7 @@
#include "inmost_dense.h"
#include "inmost_solver.h"
#include "inmost_solver_interface.h"
#include "inmost_solver_factory.h"
#include "inmost_partitioner.h"
#include "inmost_variable.h"
#include "inmost_nonlinear.h"
......
#ifndef INMOST_SOLVERFACTORY_H
#define INMOST_SOLVERFACTORY_H
#include <inmost.h>
#ifndef INMOST_SOLVER_FACTORY
#define INMOST_SOLVER_FACTORY
namespace INMOST {
struct SolverBaseFactory {
virtual SolverInterface *create() = 0;
virtual SolverInterface *copy(const SolverInterface *other) = 0;
virtual ~SolverBaseFactory() {};
};
......@@ -32,16 +28,12 @@ namespace INMOST {
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
#endif //INMOST_SOLVER_FACTORY
set(SOURCE
${SOURCE}
${CMAKE_CURRENT_SOURCE_DIR}/Solver.cpp
${CMAKE_CURRENT_SOURCE_DIR}/OrderInfo.cpp
${CMAKE_CURRENT_SOURCE_DIR}/SolverFactory.cpp
${CMAKE_CURRENT_SOURCE_DIR}/sparse.cpp
)
set(HEADER
${HEADER}
${CMAKE_CURRENT_SOURCE_DIR}/SolverFactory.h
)
add_subdirectory(solver_inner)
if (USE_SOLVER_PETSC)
......
This diff is collapsed.
......@@ -3,7 +3,6 @@
//
#include "inmost.h"
#include "SolverFactory.h"
namespace INMOST {
......
This diff is collapsed.
set(HEADER ${HEADER} ${CMAKE_CURRENT_SOURCE_DIR}/solver_prototypes.hpp
${CMAKE_CURRENT_SOURCE_DIR}/solver_bcgsl.hpp)
set(HEADER ${HEADER}
${CMAKE_CURRENT_SOURCE_DIR}/solver_prototypes.hpp
${CMAKE_CURRENT_SOURCE_DIR}/solver_bcgsl.hpp)
add_subdirectory(solver_ilu2)
add_subdirectory(solver_ddpqiluc2)
add_subdirectory(solver_mptiluc)
add_subdirectory(solver_mptilu2)
set(SOURCE ${SOURCE} PARENT_SCOPE)
set(HEADER ${HEADER} PARENT_SCOPE)
\ No newline at end of file
set(SOURCE ${SOURCE}
${CMAKE_CURRENT_SOURCE_DIR}/SolverDDPQILUC2.cpp
${CMAKE_CURRENT_SOURCE_DIR}/solver_ddpqiluc2.cpp)
set(HEADER ${HEADER}
${CMAKE_CURRENT_SOURCE_DIR}/SolverDDPQILUC2.h
${CMAKE_CURRENT_SOURCE_DIR}/solver_ddpqiluc2.hpp)
set(SOURCE ${SOURCE} PARENT_SCOPE)
set(HEADER ${HEADER} PARENT_SCOPE)
\ No newline at end of file
#include "SolverDDPQILUC2.h"
namespace INMOST {
SolverDDPQILUC2::SolverDDPQILUC2() {
Method *preconditioner = new ILUC_preconditioner(info);
solver = new KSOLVER(preconditioner, info);
matrix = NULL;
additive_schwartz_overlap = DEFAULT_ADDITIVE_SCHWARTZ_OVERLAP;
maximum_iterations = DEFAULT_MAXIMUM_ITERATIONS;
absolute_tolerance = DEFAULT_ABSOLUTE_TOLERANCE;
relative_tolerance = DEFAULT_RELATIVE_TOLERANCE;
divergence_tolerance = DEFAULT_DIVERGENCE_TOLERANCE;
preconditioner_drop_tolerance = DEFAULT_PRECONDITIONER_DROP_TOLERANCE;
preconditioner_reuse_tolerance = DEFAULT_PRECONDITIONER_REUSE_TOLERANCE;
preconditioner_rescale_iterations = DEFAULT_PRECONDITIONER_RESCALE_ITERS;
preconditioner_ddpq_tolerance = DEFAULT_PRECONDITIONER_DDPQ_TOLERANCE;
preconditioner_reorder_nonzero = DEFAULT_PRECONDITIONER_REORDER_NONZEROS;
preconditioner_condition_estimation = DEFAULT_PRECONDITIONER_CONDITION_ESTIMATION;
preconditioner_adapt_ddpq_tolerance = DEFAULT_PRECONDITIONER_ADAPT_DDPQ_TOLERANCE;
solver_gmres_substeps = DEFAULT_SOLVER_GMRES_SUBSTEPS;
}
SolverDDPQILUC2::SolverDDPQILUC2(const SolverInterface *other) {
throw INMOST::SolverUnsupportedOperation;
}
void SolverDDPQILUC2::Assign(const SolverInterface *other) {
throw INMOST::SolverUnsupportedOperation;
}
void SolverDDPQILUC2::Initialize(int *argc, char ***argv, const char *parameters_file, std::string prefix) {
}
void SolverDDPQILUC2::SetMatrix(Sparse::Matrix &A, bool ModifiedPattern, bool OldPreconditioner) {
if (matrix != NULL) {
delete matrix;
}
matrix = new Sparse::Matrix(A);
info.PrepareMatrix(*matrix, additive_schwartz_overlap);
solver->ReplaceMAT(*matrix);
solver->RealParameter(":tau") = preconditioner_drop_tolerance;
solver->RealParameter(":tau2") = preconditioner_reuse_tolerance;
solver->EnumParameter(":scale_iters") = preconditioner_rescale_iterations;
solver->RealParameter(":ddpq_tau") = preconditioner_ddpq_tolerance;
solver->EnumParameter(":reorder_nnz") = preconditioner_reorder_nonzero;
solver->EnumParameter(":estimator") = preconditioner_condition_estimation;
solver->EnumParameter(":ddpq_tau_adapt") = preconditioner_adapt_ddpq_tolerance;
if (sizeof(KSOLVER) == sizeof(BCGSL_solver)) {
solver->EnumParameter("levels") = solver_gmres_substeps;
}
if (!solver->isInitialized()) {
solver->Initialize();
}
}
bool SolverDDPQILUC2::Solve(Sparse::Vector &RHS, Sparse::Vector &SOL) {
solver->EnumParameter("maxits") = maximum_iterations;
solver->RealParameter("rtol") = relative_tolerance;
solver->RealParameter("atol") = absolute_tolerance;
solver->RealParameter("divtol") = divergence_tolerance;
return solver->Solve(RHS, SOL);
}
bool SolverDDPQILUC2::Clear() {
info.Clear();
if (matrix != NULL) {
delete matrix;
matrix = NULL;
}
if (solver != NULL) {
delete solver;
solver = NULL;
}
return true;
}
bool SolverDDPQILUC2::isMatrixSet() {
return matrix != NULL;
}
INMOST_DATA_REAL_TYPE SolverDDPQILUC2::GetPropertyReal(std::string property) const {
return solver->RealParameter(property);
}
INMOST_DATA_ENUM_TYPE SolverDDPQILUC2::GetPropertyEnum(std::string property) const {
return solver->EnumParameter(property);
}
void SolverDDPQILUC2::SetPropertyReal(std::string property, INMOST_DATA_REAL_TYPE value) {
if (property == "absolute_tolerance")
absolute_tolerance = value;
else if (property == "relative_tolerance")
relative_tolerance = value;
else if (property == "divergence_tolerance")
divergence_tolerance = value;
else if (property == "drop_tolerance")
preconditioner_drop_tolerance = value;
else if (property == "reuse_tolerance")
preconditioner_reuse_tolerance = value;
else if (property == "ddpq_tolerance")
preconditioner_ddpq_tolerance = value;
else std::cout << "Parameter " << property << " of real type is unknown" << std::endl;
}
void SolverDDPQILUC2::SetPropertyEnum(std::string property, INMOST_DATA_ENUM_TYPE value) {
if (property == "maximum_iterations")
maximum_iterations = value;
else if (property == "rescale_iterations")
preconditioner_rescale_iterations = value;
else if (property == "condition_estimation")
preconditioner_condition_estimation = value;
else if (property == "adapt_ddpq_tolerance")
preconditioner_adapt_ddpq_tolerance = value;
else if (property == "schwartz_overlap")
additive_schwartz_overlap = value;
else if (property == "gmres_substeps")
solver_gmres_substeps = value;
else if (property == "reorder_nonzeros")
preconditioner_reorder_nonzero = value;
else std::cout << "Parameter " << property << " of integral type is unknown" << std::endl;
}
const INMOST_DATA_ENUM_TYPE SolverDDPQILUC2::Iterations() const {
return solver->GetIterations();
}
const INMOST_DATA_REAL_TYPE SolverDDPQILUC2::Residual() const {
return solver->GetResidual();
}
const std::string SolverDDPQILUC2::ReturnReason() const {
return solver->GetReason();
}
const std::string SolverDDPQILUC2::SolverName() const {
return "inner_ddpqiluc2";
}
void SolverDDPQILUC2::Finalize() {
}
SolverDDPQILUC2::~SolverDDPQILUC2() {
this->Clear();
}
}
#ifndef INMOST_SOLVERDDPQILUC2_H
#define INMOST_SOLVERDDPQILUC2_H
#include <inmost.h>
#include "solver_ddpqiluc2.hpp"
#include "../solver_bcgsl.hpp"
namespace INMOST {
class SolverDDPQILUC2 : public SolverInterface {
private:
Sparse::Matrix *matrix;
KSOLVER *solver;
Solver::OrderInfo info;
INMOST_DATA_ENUM_TYPE additive_schwartz_overlap;
INMOST_DATA_ENUM_TYPE maximum_iterations;
INMOST_DATA_REAL_TYPE absolute_tolerance;
INMOST_DATA_REAL_TYPE relative_tolerance;
INMOST_DATA_REAL_TYPE divergence_tolerance;
INMOST_DATA_REAL_TYPE preconditioner_drop_tolerance;
INMOST_DATA_REAL_TYPE preconditioner_reuse_tolerance;
INMOST_DATA_ENUM_TYPE preconditioner_rescale_iterations;
INMOST_DATA_REAL_TYPE preconditioner_ddpq_tolerance;
INMOST_DATA_ENUM_TYPE preconditioner_reorder_nonzero;
INMOST_DATA_ENUM_TYPE preconditioner_condition_estimation;
INMOST_DATA_ENUM_TYPE preconditioner_adapt_ddpq_tolerance;
INMOST_DATA_ENUM_TYPE solver_gmres_substeps;
public:
SolverDDPQILUC2();
SolverDDPQILUC2(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(Sparse::Vector &RHS, Sparse::Vector &SOL);
virtual bool Clear();
virtual bool isMatrixSet();
virtual INMOST_DATA_REAL_TYPE GetPropertyReal(std::string property) const;
virtual INMOST_DATA_ENUM_TYPE GetPropertyEnum(std::string property) const;
virtual void SetPropertyReal(std::string property, INMOST_DATA_REAL_TYPE value);
virtual void SetPropertyEnum(std::string property, INMOST_DATA_ENUM_TYPE 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 ~SolverDDPQILUC2();
};
}
#endif //INMOST_SOLVERDDPQILUC2_H
#ifndef __SOLVER_DDPQILUC2__
#define __SOLVER_DDPQILUC2__
#include <iomanip>
#include "inmost_solver.h"
#include "../solver_prototypes.hpp"
class ILUC_preconditioner : public Method
{
typedef struct Interval_t
{
INMOST_DATA_ENUM_TYPE first, last;
INMOST_DATA_ENUM_TYPE Size() { return last - first; }
} Interval;
typedef std::pair<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> coord;
typedef std::pair<INMOST_DATA_REAL_TYPE, coord > wgt_coord;
typedef std::vector< wgt_coord > wgt_coords;
typedef struct row_col_t
{
Sparse::Row row, col;
INMOST_DATA_REAL_TYPE diag;
} row_col;
typedef dynarray<INMOST_DATA_ENUM_TYPE,256> levels_t;
//result of multilevel preconditioner
// |LDU F |
// A = | |
// | E C |
// LU decomposition is stored in skyline format with reversed direction
//
// For the next step of multilevel factorization we have to compute
// S = C - (E U) D (L F)
//
// U is stored by rows
// L is stored by columns
// F is stored by rows
// E is stored by rows
// C is stored by rows
//
// During LF calculation F is traversed transposed (by columns) and
// each column is solved with L.
// LF is obtained by columns.
//
// During EU calculating E is traversed by rows and
// each row is solved with U. Then re
//
// For the faster multiplication of
std::vector<Sparse::Row::entry> LU_Entries, B_Entries;
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> LU_Diag;
interval<INMOST_DATA_ENUM_TYPE, Interval> U_Address, L_Address, B_Address;
std::vector<Sparse::Row::entry> E_Entries, F_Entries;
std::vector<interval<INMOST_DATA_ENUM_TYPE, Interval> *> E_Address;
interval<INMOST_DATA_ENUM_TYPE, Interval> F_Address;
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> temp; // temporal place for solve phase
levels_t level_size; //remember size of each level
std::vector<Interval> level_interval;
//reordering information
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE > ddP,ddQ;
INMOST_DATA_ENUM_TYPE reorder_nnz, ddpq_tau_adapt, estimator;
INMOST_DATA_REAL_TYPE ddpq_tau, iluc2_tau;
INMOST_DATA_REAL_TYPE tau, eps;
INMOST_DATA_ENUM_TYPE sciters;
Sparse::Matrix * Alink;
Solver::OrderInfo * info;
bool init;
void DumpMatrix(interval<INMOST_DATA_ENUM_TYPE, Interval> & Address,
std::vector<Sparse::Row::entry> & Entries,
INMOST_DATA_ENUM_TYPE wmbeg, INMOST_DATA_ENUM_TYPE wmend,
std::string file_name);
void SwapEntries(interval<INMOST_DATA_ENUM_TYPE, Interval> & Address,
std::vector<Sparse::Row::entry> & Entries,
INMOST_DATA_ENUM_TYPE rbeg, INMOST_DATA_ENUM_TYPE rend,
INMOST_DATA_ENUM_TYPE k, INMOST_DATA_ENUM_TYPE j);
void SwapLine(interval<INMOST_DATA_ENUM_TYPE, Interval> & Line, INMOST_DATA_ENUM_TYPE i, INMOST_DATA_ENUM_TYPE j);
void SwapE(INMOST_DATA_ENUM_TYPE i, INMOST_DATA_ENUM_TYPE j);
void ReorderEF(INMOST_DATA_ENUM_TYPE mobeg,
INMOST_DATA_ENUM_TYPE cbeg,
INMOST_DATA_ENUM_TYPE cend,
INMOST_DATA_ENUM_TYPE wend,
interval<INMOST_DATA_ENUM_TYPE, bool> & donePQ,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localP,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localQ);
void inversePQ(INMOST_DATA_ENUM_TYPE wbeg,
INMOST_DATA_ENUM_TYPE wend,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localP,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localQ,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & invP,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & invQ);
void applyPQ(INMOST_DATA_ENUM_TYPE wbeg,
INMOST_DATA_ENUM_TYPE wend,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localP,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localQ,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & invP,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & invQ);
public:
INMOST_DATA_ENUM_TYPE & EnumParameter(std::string name);
INMOST_DATA_REAL_TYPE & RealParameter(std::string name);
void Copy(const Method * other);
ILUC_preconditioner(const ILUC_preconditioner & other);
ILUC_preconditioner & operator =(ILUC_preconditioner const & other);
ILUC_preconditioner(Solver::OrderInfo & info);
bool isInitialized();
bool isFinalized();
bool Initialize();
bool Finalize();
void Multiply(int level, Sparse::Vector & input, Sparse::Vector & output);
void ApplyLU(int level, Sparse::Vector & input,Sparse::Vector & output);
void ApplyB(int level, double alpha, Sparse::Vector & x, double beta, Sparse::Vector & y);
int Descend(int level, Sparse::Vector & input, Sparse::Vector & output);
int Ascend(int level, Sparse::Vector & input, Sparse::Vector & output);
bool Solve(Sparse::Vector & input, Sparse::Vector & output);
bool ReplaceMAT(Sparse::Matrix & A);
bool ReplaceSOL(Sparse::Vector & x);
bool ReplaceRHS(Sparse::Vector & b);
Method * Duplicate();
~ILUC_preconditioner();
};
#endif //__SOLVER_DDPQILUC2__
set(SOURCE ${SOURCE} ${CMAKE_CURRENT_SOURCE_DIR}/SolverILU2.cpp)
set(HEADER ${HEADER} ${CMAKE_CURRENT_SOURCE_DIR}/solver_ilu2.hpp
${CMAKE_CURRENT_SOURCE_DIR}/SolverILU2.h)
set(SOURCE ${SOURCE}
${CMAKE_CURRENT_SOURCE_DIR}/SolverILU2.cpp)
set(HEADER ${HEADER}
${CMAKE_CURRENT_SOURCE_DIR}/solver_ilu2.hpp
${CMAKE_CURRENT_SOURCE_DIR}/SolverILU2.h)
set(SOURCE ${SOURCE} PARENT_SCOPE)
set(HEADER ${HEADER} PARENT_SCOPE)
\ No newline at end of file
......@@ -13,9 +13,10 @@ namespace INMOST {
preconditioner_reuse_tolerance = DEFAULT_PRECONDITIONER_REUSE_TOLERANCE;
preconditioner_rescale_iterations = DEFAULT_PRECONDITIONER_RESCALE_ITERS;
preconditioner_fill_level = DEFAULT_PRECONDITIONER_FILL_LEVEL;
solver_gmres_substeps = DEFAULT_SOLVER_GMRES_SUBSTEPS;
Method *prec = new ILU2_preconditioner(info);
solver = new BCGS_solver(prec, info);
Method *preconditioner = new ILU2_preconditioner(info);
solver = new KSOLVER(preconditioner, info);
matrix = NULL;
}
......@@ -34,19 +35,22 @@ namespace INMOST {
}
void SolverILU2::SetMatrix(Sparse::Matrix &A, bool ModifiedPattern, bool OldPreconditioner) {
Sparse::Matrix *m = new Sparse::Matrix(A);
info.PrepareMatrix(*m, additive_schwartz_overlap);
solver->ReplaceMAT(*m);
if (matrix != NULL) {