Commit 2ab39862 by Kirill Terekhov

Solver update

Added multilevel version of crout-ilu solver with maximum product
transversal preordering on each level. The multilevel version is based
on condition number estimation, factorization of rows with large
condition number or small diagonal is postponed for the next level.

Fixes for near-zero values in maximum product transversal algorithm and
negative path-length.

Fixes for output in DrawMatrix, option for text-based output for small
matrices.
parent c7232fa2
......@@ -422,7 +422,7 @@ void draw()
for(Sparse::Row::iterator jt = it->Begin(); jt != it->End(); ++jt)
if( jt->first != it - m->Begin() )
//DrawEntry(ord->position((it - m->Begin())), m->Size() - ord->position(jt->first));//, sqrt((jt->second-min)/(max-min)));
DrawEntry((it - m->Begin()), m->Size() - jt->first);//, sqrt((jt->second-min)/(max-min)));
DrawEntry(jt->first, m->Size() - (it - m->Begin()));//, sqrt((jt->second-min)/(max-min)));
glEnd();
......@@ -488,11 +488,68 @@ int main(int argc, char ** argv)
{
if( argc < 2 )
{
std::cout << "usage: " << argv[0] << " matrix.mtx " << std::endl;
std::cout << "usage: " << argv[0] << " matrix.mtx [text]" << std::endl;
return -1;
}
m = new Sparse::Matrix();
m->Load(argv[1]);
if( argc > 2 )
{
int recw = (argc > 3 ? (atoi(argv[3])+4) : 10);
if( recw <= 4 )
{
std::cout << "Bad record width " << argv[3] << " on input, setting 10" << std::endl;
recw = 10;
}
std::string line(recw,'-');
if( std::string(argv[2]) == "text" )
{
std::cout << std::fixed << std::setprecision(recw-4);
std::cout << " ";
for(int k = 0; k < (int)m->Size(); ++k)
std::cout << std::setw(recw) << k;
std::cout << ' ' << std::endl;
std::cout << " *";
for(int k = 0; k < (int)m->Size(); ++k)
std::cout << line;
std::cout << "*" << std::endl;
for(int k = 0; k < (int)m->Size(); ++k)
{
std::cout << std::setw(3) << k << "|";
int pos = 0;
for(int l = 0; l < (int)(*m)[k].Size(); ++l)
{
int r = (*m)[k].GetIndex(l);
double v = (*m)[k].GetValue(l);
while(pos < r)
{
std::cout << std::setw(recw) << ' ';
pos++;
}
std::cout << std::setw(recw) << v;
pos++;
}
while(pos < (int)m->Size())
{
std::cout << std::setw(recw) << ' ';
pos++;
}
std::cout << "|" << std::endl;
}
std::cout << " *";
for(int k = 0; k < (int)m->Size(); ++k)
std::cout << line;
std::cout << "*" << std::endl;
}
delete m;
return 0;
}
//ord = new Reorder_ARMS(m,0,m->Size());
std::cout << "Matrix size: " << m->Size() << std::endl;
INMOST_DATA_ENUM_TYPE nnz = 0, nnzrow;
......@@ -551,4 +608,4 @@ int main(int argc, char ** argv)
glutPostRedisplay();
glutMainLoop();
}
\ No newline at end of file
}
......@@ -28,6 +28,7 @@ namespace INMOST
static const Type INNER_ILU2; ///< inner Solver based on BiCGStab(L) solver with second order ILU factorization as preconditioner.
static const Type 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.
static const Type 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.
static const Type INNER_MLMPTILUC; ///< inner Solver based on BiCGStab(L) solver with second order Crout-ILU with inversed-based condition estimation and maximum product transversal reordering as preconditioner.
static const Type INNER_MPTILU2; ///< inner Solver based on BiCGStab(L) solver with second order ILU and maximum product transversal reordering as preconditioner.
static const Type Trilinos_Aztec; ///< external Solver AztecOO from Trilinos package.
static const Type Trilinos_Belos; ///< external Solver Belos from Trilinos package, currently without preconditioner.
......@@ -275,26 +276,26 @@ namespace INMOST
/// - "maximum_iterations" - total number of iterations
/// - "schwartz_overlap" - number of overlapping levels for additive schwartz method,
/// works for:
/// INNER_ILU2, INNER_MLILUC
/// INNER_ILU2, INNER_MPTILU2, INNER_MPTILUC, INNER_MLMPTILUC, INNER_DDPQILUC
/// Trilinos_Aztec, Trilinos_Belos, Trilinos_ML, Trilinos_Ifpack
/// PETSc
/// - "gmres_substeps" - number of gmres steps performed after each bicgstab step,
/// works for:
/// INNER_ILU2, INNER_MLILUC
/// INNER_ILU2, INNER_MPTILU2, INNER_MPTILUC, INNER_MLMPTILUC, INNER_DDPQILUC
/// - "reorder_nonzeros" - place sparser rows at the beggining of matrix during reordering,
/// works for:
/// INNER_MLILUC
/// INNER_DDPQILUC
/// - "rescale_iterations" - number of iterations for two-side matrix rescaling,
/// works for:
/// INNER_ILU2, INNER_MLILUC
/// INNER_ILU2, INNER_MPTILU2, INNER_MPTILUC, INNER_MLMPTILUC, INNER_DDPQILUC
/// - "condition_estimation" - exploit condition estimation of inversed factors to adapt
/// drop and reuse tolerances,
/// works for:
/// INNER_MLILUC
/// INNER_MPTILUC, INNER_MLMPTILUC, INNER_DDPQILUC
/// - "adapt_ddpq_tolerance" - adapt ddpq tolerance depending from the complexity
/// of calculation of Schur complement,
/// works for:
/// INNER_MLILUC
/// INNER_DDPQILUC
/// Set the solver parameter of the real type.
///
/// Parameters:
......@@ -306,7 +307,7 @@ namespace INMOST
/// ||A x(i) - b|| > divergence_tolerance
/// - "drop_tolerance" - tolerance for dropping values during incomplete factorization,
/// works for:
/// INNER_ILU2, INNER_MLILUC
/// INNER_ILU2, INNER_MPTILU2, INNER_MPTILUC, INNER_MLMPTILUC, INNER_DDPQILUC
/// Trilinos_Aztec, Trilinos_Ifpack
/// PETSc
/// - "reuse_tolerance" - tolerance for reusing values during incomplete factorization,
......@@ -315,7 +316,7 @@ namespace INMOST
/// value should be less then "drop_tolerance",
/// typical value is drop_tolerance^2,
/// works for:
/// INNER_ILU2, INNER_MLILUC
/// INNER_ILU2, INNER_MPTILU2, INNER_MPTILUC, INNER_MLMPTILUC, INNER_DDPQILUC
/// - "ddpq_tolerance" - by this tolerance most diagonnaly-dominant elements will be selected
/// to form the next level of factorization, the closer the tolerance
/// is to one the smaller will be the level. Actual rule is:
......@@ -323,11 +324,19 @@ namespace INMOST
/// A(imax,jmax)/(sum(A(imax,:))+sum(A(:,jmax))-A(imax,jmax))
/// where on imax, jmax maximum is reached.
/// works for:
/// INNER_MLILUC
/// INNER_DDPQILUC
/// - "fill_level" - level of fill for ILU-type preconditioners,
/// works for:
/// INNER_ILU2 (if LFILL is defined in solver_ilu2.hpp)
/// Trilinos, Trilinos_Ifpack
/// - "pivot_condition" - delay factoriation of the row of the matrix to the next level, if
/// the estimated condition number is above prescribed value.
/// works for:
/// INNER_MLMPTILUC
/// - "pivot_diag" - delay factoriation of the row of the matrix to the next level, if
/// the inverted diagonal value is above the prescribed value.
/// works for:
/// INNER_MLMPTILUC
/// @see Solver::GetParameter
void SetParameter(std::string name, std::string value);
/// Return the number of iterations performed by the last solution.
......
......@@ -3,6 +3,7 @@
#include "solver_inner/solver_ilu2/SolverILU2.h"
#include "solver_inner/solver_ddpqiluc2/SolverDDPQILUC2.h"
#include "solver_inner/solver_mptiluc/SolverMPTILUC.h"
#include "solver_inner/solver_mlmptiluc/SolverMLMPTILUC.h"
#include "solver_inner/solver_mptilu2/SolverMPTILU2.h"
#if defined(USE_SOLVER_PETSC)
......@@ -51,6 +52,7 @@ namespace INMOST {
const Solver::Type Solver::INNER_ILU2 = "inner_ilu2";
const Solver::Type Solver::INNER_DDPQILUC = "inner_ddpqiluc2";
const Solver::Type Solver::INNER_MPTILUC = "inner_mptiluc";
const Solver::Type Solver::INNER_MLMPTILUC = "inner_mlmptiluc";
const Solver::Type Solver::INNER_MPTILU2 = "inner_mptilu2";
const Solver::Type Solver::Trilinos_Aztec = "trilinos_aztec";
const Solver::Type Solver::Trilinos_Belos = "trilinos_belos";
......@@ -62,10 +64,12 @@ namespace INMOST {
const Solver::Type Solver::K3BIILU2 = "k3biilu2";
const Solver::Type Solver::SUPERLU = "superlu";
SolverInterface *SolverMaster::getSolver(std::string name) {
SolverInterface *SolverMaster::getSolver(std::string name)
{
if (name == "inner_ilu2") return new SolverILU2();
if (name == "inner_ddpqiluc2") return new SolverDDPQILUC2();
if (name == "inner_mptiluc") return new SolverMPTILUC();
if (name == "inner_mlmptiluc") return new SolverMLMPTILUC();
if (name == "inner_mptilu2") return new SolverMPTILU2();
#if defined(USE_SOLVER_PETSC)
if (name == "petsc") return new SolverPETSc();
......@@ -91,11 +95,13 @@ namespace INMOST {
throw INMOST::SolverNotFound;
}
std::vector<std::string> SolverMaster::getAvailableSolvers() {
std::vector<std::string> SolverMaster::getAvailableSolvers()
{
std::vector<std::string> s;
s.push_back("inner_ilu2");
s.push_back("inner_ddpqiluc2");
s.push_back("inner_mptiluc");
s.push_back("inner_mlmptiluc");
s.push_back("inner_mptilu2");
#if defined(USE_SOLVER_PETSC)
s.push_back("petsc");
......@@ -121,16 +127,20 @@ namespace INMOST {
return s;
}
bool SolverMaster::isSolverAvailable(std::string name) {
try {
bool SolverMaster::isSolverAvailable(std::string name)
{
try
{
SolverInterface *s = SolverMaster::getSolver(name);
delete s;
return true;
} catch (...) {
}
catch (...)
{
return false;
}
}
}
#endif //USE_SOLVER
\ No newline at end of file
#endif //USE_SOLVER
......@@ -9,6 +9,7 @@ set(SOURCE ${SOURCE}
add_subdirectory(solver_ilu2)
add_subdirectory(solver_ddpqiluc2)
add_subdirectory(solver_mptiluc)
add_subdirectory(solver_mlmptiluc)
add_subdirectory(solver_mptilu2)
set(SOURCE ${SOURCE} PARENT_SCOPE)
......
set(SOURCE ${SOURCE}
${CMAKE_CURRENT_SOURCE_DIR}/SolverMLMPTILUC.cpp
${CMAKE_CURRENT_SOURCE_DIR}/solver_mlmtiluc2.cpp)
set(HEADER ${HEADER}
${CMAKE_CURRENT_SOURCE_DIR}/SolverMLMPTILUC.h
${CMAKE_CURRENT_SOURCE_DIR}/solver_mlmtiluc2.hpp)
set(SOURCE ${SOURCE} PARENT_SCOPE)
set(HEADER ${HEADER} PARENT_SCOPE)
\ No newline at end of file
#include "SolverMLMPTILUC.h"
namespace INMOST
{
SolverMLMPTILUC::SolverMLMPTILUC()
{
Method *preconditioner = new MLMTILUC_preconditioner(info);
solver = new KSOLVER(preconditioner, info);
matrix = NULL;
rescale_iterations = 6;
condition_estimation = 1;
schwartz_overlap = 1;
gmres_substeps = 2;
reorder_nnz = 1;
drop_tolerance = 0.005;
reuse_tolerance = 0.00005;
pivot_condition = 1.0e+6;
pivot_diag = 1.0e+6;
fill_level = 3;
}
SolverMLMPTILUC::SolverMLMPTILUC(const SolverInterface *other)
{
//You should not really want to copy solver's information
throw INMOST::SolverUnsupportedOperation;
}
void SolverMLMPTILUC::SetMatrix(Sparse::Matrix &A, bool ModifiedPattern, bool OldPreconditioner)
{
if (matrix != NULL)
delete matrix;
matrix = new Sparse::Matrix(A);
info.PrepareMatrix(*matrix, schwartz_overlap);
solver->ReplaceMAT(*matrix);
solver->RealParameter(":tau") = drop_tolerance;
solver->RealParameter(":tau2") = reuse_tolerance;
solver->EnumParameter(":scale_iters") = rescale_iterations;
solver->EnumParameter(":estimator") = condition_estimation;
solver->RealParameter(":pivot_cond") = pivot_condition;
solver->RealParameter(":pivot_diag") = pivot_diag;
if (sizeof(KSOLVER) == sizeof(BCGSL_solver))
solver->EnumParameter("levels") = gmres_substeps;
if (!solver->isInitialized())
solver->Initialize();
}
void SolverMLMPTILUC::SetParameter(std::string name, std::string value)
{
const char *val = value.c_str();
if (name == "rescale_iterations") rescale_iterations = static_cast<INMOST_DATA_ENUM_TYPE>(atoi(val));
else if (name == "condition_estimation") condition_estimation = static_cast<INMOST_DATA_ENUM_TYPE>(atoi(val));
else if (name == "schwartz_overlap") schwartz_overlap = static_cast<INMOST_DATA_ENUM_TYPE>(atoi(val));
else if (name == "gmres_substeps") gmres_substeps = static_cast<INMOST_DATA_ENUM_TYPE>(atoi(val));
else if (name == "reorder_nonzeros") reorder_nnz = static_cast<INMOST_DATA_ENUM_TYPE>(atoi(val));
else if (name == "fill_level") fill_level = static_cast<INMOST_DATA_ENUM_TYPE>(atoi(val));
else if (name == "drop_tolerance") drop_tolerance = atof(val);
else if (name == "reuse_tolerance") reuse_tolerance = atof(val);
else if (name == "pivot_condition") pivot_condition = atof(val);
else if (name == "pivot_diag") pivot_diag = atof(val);
else SolverInner::SetParameter(name, value);
}
const std::string SolverMLMPTILUC::SolverName() const
{
return "inner_mptiluc";
}
SolverMLMPTILUC::~SolverMLMPTILUC() {}
}
#ifndef INMOST_SOLVERMLMPTILUC_H
#define INMOST_SOLVERMLMPTILUC_H
#include <inmost.h>
#include "solver_mlmtiluc2.hpp"
#include "../SolverInner.h"
namespace INMOST
{
class SolverMLMPTILUC : public SolverInner
{
INMOST_DATA_ENUM_TYPE rescale_iterations, condition_estimation, schwartz_overlap, gmres_substeps, reorder_nnz, fill_level;
INMOST_DATA_REAL_TYPE drop_tolerance, reuse_tolerance, pivot_condition, pivot_diag;
public:
SolverMLMPTILUC();
SolverMLMPTILUC(const SolverInterface *other);
virtual void SetMatrix(Sparse::Matrix &A, bool ModifiedPattern, bool OldPreconditioner);
virtual void SetParameter(std::string name, std::string value);
virtual const std::string SolverName() const;
virtual ~SolverMLMPTILUC();
};
}
#endif //INMOST_SOLVERMLMPTILUC_H
This source diff could not be displayed because it is too large. You can view the blob instead.
#ifndef __SOLVER_MLMTILUC2__
#define __SOLVER_MLMTILUC2__
#include <iomanip>
#include "inmost_solver.h"
#include "../solver_prototypes.hpp"
class MLMTILUC_preconditioner : public Method
{
typedef struct Interval_t
{
INMOST_DATA_ENUM_TYPE first, last;
INMOST_DATA_ENUM_TYPE Size() { return last - first; }
} Interval;
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;
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<interval<INMOST_DATA_ENUM_TYPE, Interval> *> F_Address, E_Address;
std::vector<Sparse::Row::entry> E_Entries, F_Entries;
levels_t level_size; //remember size of each level
std::vector<Interval> level_interval;
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> temp; // temporal place for solve phase
//reordering information
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE > ddP,ddQ;
INMOST_DATA_REAL_TYPE condestL, condestU;
INMOST_DATA_ENUM_TYPE estimator;
INMOST_DATA_REAL_TYPE iluc2_tau;
INMOST_DATA_REAL_TYPE pivot_cond;
INMOST_DATA_REAL_TYPE pivot_diag;
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 CheckOrder(interval<INMOST_DATA_ENUM_TYPE, Interval> & Address,
std::vector<Sparse::Row::entry> & Entries,
INMOST_DATA_ENUM_TYPE rbeg, INMOST_DATA_ENUM_TYPE rend);
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);
void ReorderEF(INMOST_DATA_ENUM_TYPE wbeg,
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);
public:
INMOST_DATA_ENUM_TYPE & EnumParameter(std::string name);
INMOST_DATA_REAL_TYPE & RealParameter(std::string name);
void Copy(const Method * other);
MLMTILUC_preconditioner(const MLMTILUC_preconditioner & other);
MLMTILUC_preconditioner & operator =(MLMTILUC_preconditioner const & other);
MLMTILUC_preconditioner(Solver::OrderInfo & info);
bool isInitialized();
bool isFinalized();
bool Initialize();
bool Finalize();
void ApplyB(double alpha, Sparse::Vector & x, double beta, Sparse::Vector & y);
int Descend(int level, Sparse::Vector & inout);
int Ascend(int level, Sparse::Vector & inout);
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();
~MLMTILUC_preconditioner();
};
#endif //__SOLVER_MLMTILUC2__
......@@ -345,7 +345,7 @@ void MTILU2_preconditioner::DumpMatrix(interval<INMOST_DATA_ENUM_TYPE, INMOST_DA
for (INMOST_DATA_ENUM_TYPE it = C_Address[k]; it < C_Address[k+1]; ++it)
{
INMOST_DATA_ENUM_TYPE i = C_Entries[it].first;
if( Cmax[i] == 0.0 )
if( Cmax[i] == 0.0 || C_Entries[it].second == 0.0 )
C_Entries[it].second = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
else
{
......@@ -549,7 +549,11 @@ void MTILU2_preconditioner::DumpMatrix(interval<INMOST_DATA_ENUM_TYPE, INMOST_DA
for(INMOST_DATA_ENUM_TYPE qt = 0; qt < A[k].Size(); ++qt)
{
INMOST_DATA_ENUM_TYPE i = A[k].GetIndex(qt), j = Perm[i];
INMOST_DATA_REAL_TYPE l = exp(V[k]), u = exp(U[i])/Cmax[i];
INMOST_DATA_REAL_TYPE l,u;
if( V[k] == std::numeric_limits<INMOST_DATA_REAL_TYPE>::max() ) l = 1;
else l = exp(V[k]);
if( U[i] == std::numeric_limits<INMOST_DATA_REAL_TYPE>::max() || Cmax[i] == 0 ) u = 1;
else u = exp(U[i])/Cmax[i];
DL[k] = l;
DR[i] = u;
B_Entries[cnt++] = Sparse::Row::make_entry(j, l*u*A[k].GetValue(qt));
......@@ -1090,4 +1094,4 @@ void MTILU2_preconditioner::DumpMatrix(interval<INMOST_DATA_ENUM_TYPE, INMOST_DA
info->Accumulate(output);
return true;
}
#endif
\ No newline at end of file
#endif
......@@ -5,11 +5,11 @@
#include "solver_mtiluc2.hpp"
#include <sstream>
#include <deque>
#define REPORT_ILU
//#define REPORT_ILU
//#undef REPORT_ILU
//#define REPORT_ILU_PROGRESS
//#define REPORT_ILU_END
#define REPORT_ILU_SUMMARY
//#define REPORT_ILU_SUMMARY
//#undef REPORT_ILU_PROGRESS
//#define USE_OMP
......@@ -93,7 +93,7 @@ public:
};
class BinaryHeap
class BinaryHeap1
{
INMOST_DATA_REAL_TYPE * Base;
......@@ -189,13 +189,13 @@ public:
BalanceHeap(0);
return Ret;
}
BinaryHeap(INMOST_DATA_REAL_TYPE * Base, INMOST_DATA_ENUM_TYPE Size)
BinaryHeap1(INMOST_DATA_REAL_TYPE * Base, INMOST_DATA_ENUM_TYPE Size)
: Base(Base)
{
Position.resize(Size,ENUMUNDEF);
Array.reserve(4096);
}
~BinaryHeap()
~BinaryHeap1()
{
}
};
......@@ -757,7 +757,7 @@ public:
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> AugmentPosition(wbeg,wend,ENUMUNDEF);
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> ColumnPosition(wbeg,wend,ENUMUNDEF);
BinaryHeap Heap(&Dist[wbeg],wend-wbeg);
BinaryHeap1 Heap(&Dist[wbeg],wend-wbeg);
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////// Arrays initialization ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
......@@ -790,7 +790,7 @@ public:
for (INMOST_DATA_ENUM_TYPE it = A_Address[k].first; it != A_Address[k].last; ++it)
{
i = A_Entries[it].first;
if( Cmax[i] == 0.0 )
if( Cmax[i] == 0.0 || C_Entries[it] == 0.0 )
C_Entries[it] = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
else
{
......@@ -965,45 +965,27 @@ public:
//T = Timer();
for (k = cbeg; k < cend; ++k)
{
//B_Address[k].first = static_cast<INMOST_DATA_ENUM_TYPE>(B_Entries.size());
bool flip_sign = false;
if( V[k] == std::numeric_limits<INMOST_DATA_REAL_TYPE>::max() ) l = 1;
else l = exp(V[k]);
if( U[k] == std::numeric_limits<INMOST_DATA_REAL_TYPE>::max() || Cmax[k] == 0 ) u = 1;
else u = exp(U[k])/Cmax[k];
//if( l != l || fabs(l) < 1.0e-12 || isnan(l) || isinf(l) ) std::cout << "k " << k << " l is " << l << " V " << V[k] << std::endl;
//if( u != u || fabs(u) < 1.0e-12 || isnan(u) || isinf(u) ) std::cout << "k " << k << " u is " << u << " U " << U[k] << " Cmax " << Cmax[k] << std::endl;
DL[k] = l;
DR[k] = u;
bool flip_sign = false;
for (INMOST_DATA_ENUM_TYPE jt = A_Address[k].first; jt != A_Address[k].last; ++jt)
{
i = A_Entries[jt].first;
j = Perm[A_Entries[jt].first];
l = exp(V[k]);
u = exp(U[i])/Cmax[i];
DL[k] = l;
DR[i] = u;
//l = l*u;
//l = l*u;//exp(U[i]+V[k])/Cmax[i];
//l = pow(10,U[i]+V[k])/Cmax[i];
//l = exp(U[i]+V[k])/Cmax[i];
//A_Entries[jt].first = j;
//A_Entries[jt].second *= l;
//B_Entries.push_back(Sparse::Row::make_entry(j, l*A_Entries[jt].second));
if( U[i] == std::numeric_limits<INMOST_DATA_REAL_TYPE>::max() || Cmax[i] == 0 ) u = 1;
else u = exp(U[i])/Cmax[i];
if( j == k && (l*A_Entries[jt].second*u) < 0.0 ) flip_sign = true;
}
//B_Address[k].last = static_cast<INMOST_DATA_ENUM_TYPE>(B_Entries.size());
if( flip_sign )
{
DL[k] *= -1;
//for(Li = B_Address[k].first; Li < B_Address[k].last; ++Li)
// B_Entries[Li].second *= -1;
}
//std::sort(B_Entries.begin() + B_Address[k].first, B_Entries.end());
//std::sort(A_Entries.begin() + A_Address[k].first, A_Entries.begin() + A_Address[k].last);
if( flip_sign ) DL[k] *= -1;
}
//printf("Reorder matrix %lf\n",Timer()-T);
......@@ -3349,7 +3331,7 @@ swap_algorithm:
printf(" swap %f (%6.2f%%)\n", tswap, 100.0*tswap / ttotal);
printf(" cond %f (%6.2f%%)\n", testimator, 100.0*testimator / ttotal);
#if defined(ILUC2)
printf("nnz A %d LU %d LU2 %d swaps %d\n",nzA,nzLU,nzLU2, swaps);
printf("nnz A %d LU %d LU2 %d\n",nzA,nzLU,nzLU2);//, swaps);
#else
printf("nnz A %d LU %d swaps %d\n",nzA,nzLU,swaps);
#endif
......
......@@ -4,13 +4,13 @@ endif(USE_MESH)
if(USE_AUTODIFF)
add_subdirectory(autodiff_test000)
add_subdirectory(autodiff_test001)
#add_subdirectory(autodiff_test001)
add_subdirectory(autodiff_test002)
endif(USE_AUTODIFF)
if(USE_SOLVER)
add_subdirectory(solver_test000)
#add_subdirectory(solver_test001)
add_subdirectory(solver_test001)
add_subdirectory(solver_test002)
add_subdirectory(solver_test003)
endif(USE_SOLVER)
......
......@@ -39,9 +39,10 @@ if(USE_SOLVER_SUPERLU)
target_link_libraries(solver_test001 ${SUPERLU_LIBRARIES})
endif()
set(SOLVERS inner_ilu2
inner_mptilu2
inner_mptiluc)
set(SOLVERS #inner_ilu2
#inner_mptilu2
inner_mptiluc
inner_mlmptiluc)
foreach(solver ${SOLVERS})
......@@ -179,4 +180,4 @@ foreach(test ${FIDAP_TESTS})
add_test(NAME solver_test001_sparsekit_fidap_${solver}_${test} COMMAND $<TARGET_FILE:solver_test001> ${CMAKE_CURRENT_SOURCE_DIR}/matrices/sparskit/fidap/${test}.mtx ${solver} ${CMAKE_CURRENT_SOURCE_DIR}/matrices/sparskit/fidap/${test}_rhs1.mtx)
endforeach()
endforeach(solver ${SOLVERS})
\ No newline at end of file
endforeach(solver ${SOLVERS})
......@@ -33,12 +33,17 @@ int main(int argc, char ** argv)
{
std::string stype(argv[2]);
for(size_t k = 0; k < stype.size(); ++k) stype[k] = tolower(stype[k]);
/*
if( stype == "inner_mptilu2" )
type = Solver::INNER_MPTILU2;
else if( stype == "inner_mptiluc" )
type = Solver::INNER_MPTILUC;
else if( stype == "inner_ddpqiluc" )
type = Solver::INNER_DDPQILUC;
else
*/
if( Solver::isSolverAvailable(stype) )
type = stype;
else if( stype != "inner_ilu2" && stype != "*" )
{
std::cout << "Unknown type of the solver " << stype << std::endl;
......@@ -54,6 +59,17 @@ int main(int argc, char ** argv)
mat.Load(std::string(argv[1])); //if interval parameters not set, matrix will be divided automatically
BARRIER
if( !rank ) std::cout << "load matrix: " << Timer() - t << std::endl;
/*
b.SetInterval(0,mat.Size());
x.SetInterval(0,mat.Size());
for(int k = 0; k < (int)mat.Size(); ++k)
x[k] = (1+k)*10;
mat.MatVec(1.0,x,0.0,b);
for(int k = 0; k < (int)mat.Size(); ++k)
x[k] = 0;
*/
//mat.Save("test.mtx");
t = Timer();
if( argc > 3 )
......@@ -61,7 +77,7 @@ int main(int argc, char ** argv)
//std::cout << rank << " load vector from " << std::string(argv[3]) << std::endl;
b.Load(std::string(argv[3])); // Load RHS vector
}
else // Set local RHS to 1 if it was not specified
else if( b.Empty() )// Set local RHS to 1 if it was not specified
{
INMOST_DATA_ENUM_TYPE mbeg,mend,k;
mat.GetInterval(mbeg,mend);
......@@ -82,10 +98,11 @@ int main(int argc, char ** argv)
s.SetParameterEnum("condition_estimation",1);
s.SetParameterReal("relative_tolerance",1.0e-9);
s.SetParameterReal("absolute_tolerance",1.0e-16);