Commit e6f3aa27 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

add support for verbosity parameter ininmost solvers

parent 23bbbe64
......@@ -14,10 +14,6 @@
#define PSEUDOINVERSE // same trick as in petsc with pseudoinverse
//#define USE_LAPACK_SVD // use lapack's dgesvd routine instead of built-in svdnxn
//#if !defined(NDEBUG)
#define REPORT_RESIDUAL
//#endif
//#define USE_OMP
namespace INMOST
......@@ -437,7 +433,7 @@ namespace INMOST
{
INMOST_DATA_REAL_TYPE length;
INMOST_DATA_REAL_TYPE rtol, atol, divtol, last_resid;
INMOST_DATA_ENUM_TYPE iters, maxits, l, last_it;
INMOST_DATA_ENUM_TYPE iters, maxits, l, last_it, verbosity;
INMOST_DATA_REAL_TYPE resid;
INMOST_DATA_REAL_TYPE * tau, * sigma, * gamma, *theta1, * theta2, * theta3;
Sparse::Vector r_tilde, x0, t, * u, * r;
......@@ -446,6 +442,7 @@ namespace INMOST
std::string reason;
Solver::OrderInfo * info;
bool init;
public:
INMOST_DATA_ENUM_TYPE GetIterations() {return last_it;}
INMOST_DATA_REAL_TYPE GetResidual() {return last_resid;}
......@@ -468,6 +465,7 @@ namespace INMOST
if (prec != NULL) return prec->EnumParameter(name.substr(1, name.size() - 1));
}
if (name == "maxits") return maxits;
else if (name == "verbosity" ) return verbosity;
else if (name == "levels")
{
if( init ) throw - 1; //solver was already initialized, value should not be changed
......@@ -477,7 +475,7 @@ namespace INMOST
throw - 1;
}
BCGSL_solver(Method * prec, Solver::OrderInfo & info)
:rtol(1e-8), atol(1e-9), divtol(1e+40), maxits(1500),l(2),prec(prec),info(&info)
:rtol(1e-8), atol(1e-9), divtol(1e+40), maxits(1500),l(2),prec(prec),info(&info), verbosity(0)
{
Alink = NULL;
init = false;
......@@ -533,6 +531,7 @@ namespace INMOST
resid = b->resid;
Alink = b->Alink;
info = b->info;
verbosity = b->verbosity;
if (init) Finalize();
if (b->prec != NULL)
{
......@@ -618,16 +617,18 @@ namespace INMOST
for(INMOST_DATA_INTEGER_TYPE k = ivbeg; k < ivend; ++k) // r_tilde = r[0] / dot(r[0],r[0])
r_tilde[k] /= resid;
last_it = 0;
#if defined(REPORT_RESIDUAL)
if( info->GetRank() == 0 )
if( verbosity )
{
//std::cout << "iter " << last_it << " residual " << resid << std::endl;
//std::cout << "iter " << last_it << " resid " << resid << "\r";
//printf("iter %3d resid %12g | %12g relative %12g | %12g\r", last_it, resid, atol, resid / resid0, rtol);
printf("iter %3d resid %12g | %g\r", last_it, resid, atol);
fflush(stdout);
if( info->GetRank() == 0 )
{
//std::cout << "iter " << last_it << " residual " << resid << std::endl;
//std::cout << "iter " << last_it << " resid " << resid << "\r";
//printf("iter %3d resid %12g | %12g relative %12g | %12g\r", last_it, resid, atol, resid / resid0, rtol);
printf("iter %3d resid %12g | %g\r", last_it, resid, atol);
fflush(stdout);
}
}
#endif
bool halt = false;
if( last_resid < atol || last_resid < rtol*resid0 )
......@@ -1212,21 +1213,23 @@ namespace INMOST
resid = sqrt(resid/rhs_norm);
}
#if defined(REPORT_RESIDUAL)
if( info->GetRank() == 0 )
if( verbosity )
{
//std::cout << "iter " << last_it << " residual " << resid << " time " << tt << " matvec " << ts*0.5/l << " precond " << tp*0.5/l << std::endl;
//std::cout << "iter " << last_it << " resid " << resid << "\r";
//printf("iter %3d resid %12g | %12g relative %12g | %12g\r", last_it, resid, atol, resid / resid0, rtol);
if( info->GetRank() == 0 )
{
//std::cout << "iter " << last_it << " residual " << resid << " time " << tt << " matvec " << ts*0.5/l << " precond " << tp*0.5/l << std::endl;
//std::cout << "iter " << last_it << " resid " << resid << "\r";
//printf("iter %3d resid %12g | %12g relative %12g | %12g\r", last_it, resid, atol, resid / resid0, rtol);
#if defined(USE_OMP)
#pragma omp single
#endif
{
printf("iter %3d resid %12g | %g\r", last_it, resid, atol);
fflush(stdout);
{
printf("iter %3d resid %12g | %g\r", last_it, resid, atol);
fflush(stdout);
}
}
}
#endif
#if defined(USE_OMP)
#pragma omp single
#endif
......@@ -1301,7 +1304,7 @@ namespace INMOST
class BCGS_solver : public IterativeMethod
{
INMOST_DATA_REAL_TYPE rtol, atol, divtol, last_resid;
INMOST_DATA_ENUM_TYPE iters, maxits, last_it;
INMOST_DATA_ENUM_TYPE iters, maxits, last_it,verbosity;
INMOST_DATA_REAL_TYPE resid;
Sparse::Vector r0, p, y, s, t, z, r, v;
Sparse::Matrix * Alink;
......@@ -1331,11 +1334,12 @@ namespace INMOST
if (prec != NULL) return prec->EnumParameter(name.substr(1, name.size() - 1));
}
if (name == "maxits") return maxits;
else if (name == "verbosity") return verbosity;
else if (prec != NULL) return prec->EnumParameter(name);
throw - 1;
}
BCGS_solver(Method * prec, Solver::OrderInfo & info)
:rtol(1e-8), atol(1e-11), divtol(1e+40), iters(0), maxits(1500),prec(prec),info(&info)
:rtol(1e-8), atol(1e-11), divtol(1e+40), iters(0), maxits(1500),prec(prec),info(&info),verbosity(0)
{
init = false;
}
......@@ -1377,6 +1381,7 @@ namespace INMOST
last_it = b->last_it;
resid = b->resid;
Alink = b->Alink;
verbosity = b->verbosity;
if (b->prec != NULL)
{
if (prec == NULL) prec = b->prec->Duplicate();
......@@ -1444,16 +1449,17 @@ namespace INMOST
for(INMOST_DATA_INTEGER_TYPE k = ivbeg; k < ivend; ++k) // r_tilde = r[0] / dot(r[0],r[0])
r0[k] /= resid;
last_it = 0;
#if defined(REPORT_RESIDUAL)
if( info->GetRank() == 0 )
if( verbosity )
{
//std::cout << "iter " << last_it << " residual " << resid << std::endl;
//std::cout << "iter " << last_it << " resid " << resid << "\r";
//printf("iter %3d resid %12g | %12g relative %12g | %12g\r",last_it,resid,atol,resid/resid0,rtol);
printf("iter %3d resid %12g | %g\r", last_it, resid, atol);
fflush(stdout);
if( info->GetRank() == 0 )
{
//std::cout << "iter " << last_it << " residual " << resid << std::endl;
//std::cout << "iter " << last_it << " resid " << resid << "\r";
//printf("iter %3d resid %12g | %12g relative %12g | %12g\r",last_it,resid,atol,resid/resid0,rtol);
printf("iter %3d resid %12g | %g\r", last_it, resid, atol);
fflush(stdout);
}
}
#endif
bool halt = false;
if( last_resid < atol || last_resid < rtol*resid0 )
......@@ -1775,22 +1781,23 @@ namespace INMOST
{
resid = sqrt(resid);
}
#if defined(REPORT_RESIDUAL)
if( info->GetRank() == 0 )
if( verbosity )
{
//std::cout << "iter " << last_it << " residual " << resid << " time " << tt << " matvec " << ts*0.5 << " precond " << tp*0.5 << std::endl;
//std::cout << "iter " << last_it << " resid " << resid << "\r";
//printf("iter %3d resid %12g | %12g relative %12g | %12g\r", last_it, resid, atol, resid / resid0, rtol);
if( info->GetRank() == 0 )
{
//std::cout << "iter " << last_it << " residual " << resid << " time " << tt << " matvec " << ts*0.5 << " precond " << tp*0.5 << std::endl;
//std::cout << "iter " << last_it << " resid " << resid << "\r";
//printf("iter %3d resid %12g | %12g relative %12g | %12g\r", last_it, resid, atol, resid / resid0, rtol);
#if defined(USE_OMP)
#pragma omp single
#endif
{
printf("iter %3d resid %12g | %g\r", last_it, resid, atol);
//printf("iter %3d resid %12g | %g rho %e beta %e alpha %e omega %e\n", last_it, resid, atol,rho,beta,alpha,omega);
fflush(stdout);
{
printf("iter %3d resid %12g | %g\r", last_it, resid, atol);
//printf("iter %3d resid %12g | %g rho %e beta %e alpha %e omega %e\n", last_it, resid, atol,rho,beta,alpha,omega);
fflush(stdout);
}
}
}
#endif
#if defined(USE_OMP)
#pragma omp single
#endif
......
......@@ -13,6 +13,7 @@ namespace INMOST {
drop_tolerance = 0.005;
reuse_tolerance = 0.00005;
fill_level = 3;
verbosity = 0;
}
SolverILU2::SolverILU2(const SolverInterface *other) {
......@@ -33,6 +34,8 @@ namespace INMOST {
solver->RealParameter(":tau2") = reuse_tolerance;
solver->EnumParameter(":scale_iters") = rescale_iterations;
solver->EnumParameter(":fill") = fill_level;
solver->EnumParameter(":verbosity") = verbosity;
solver->EnumParameter("verbosity") = verbosity;
if (sizeof(KSOLVER) == sizeof(BCGSL_solver)) {
solver->EnumParameter("levels") = gmres_substeps;
......@@ -55,6 +58,7 @@ namespace INMOST {
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 == "verbosity") verbosity = static_cast<INMOST_DATA_ENUM_TYPE>(atoi(val));
else SolverInner::SetParameter(name, value);
}
......
......@@ -8,7 +8,7 @@
namespace INMOST {
class SolverILU2 : public SolverInner {
INMOST_DATA_ENUM_TYPE rescale_iterations, schwartz_overlap, gmres_substeps, reorder_nnz, fill_level;
INMOST_DATA_ENUM_TYPE rescale_iterations, schwartz_overlap, gmres_substeps, reorder_nnz, fill_level, verbosity;
INMOST_DATA_REAL_TYPE drop_tolerance, reuse_tolerance;
public:
SolverILU2();
......@@ -26,4 +26,4 @@ namespace INMOST {
}
#endif
\ No newline at end of file
#endif
......@@ -6,8 +6,6 @@
#include "inmost_solver.h"
#include "../solver_prototypes.hpp"
//#define REPORT_ILU
//#define REPORT_ILU_PROGRESS
using namespace INMOST;
#define DEFAULT_TAU 0.005
......@@ -27,11 +25,12 @@ private:
std::vector<INMOST_DATA_REAL_TYPE> luv;
std::vector<INMOST_DATA_ENUM_TYPE> lui;
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> ilu, iu;
INMOST_DATA_ENUM_TYPE Lfill;
INMOST_DATA_ENUM_TYPE Lfill,verbosity;
INMOST_DATA_REAL_TYPE tau, tau2;
Sparse::Vector DL, DR;
INMOST_DATA_ENUM_TYPE nnz, sciters;
bool init;
public:
INMOST_DATA_REAL_TYPE &RealParameter(std::string name)
{
......@@ -44,6 +43,7 @@ public:
{
if (name == "fill") return Lfill;
else if (name == "scale_iters") return sciters;
else if (name == "verbosity" ) return verbosity;
throw -1;
}
......@@ -54,6 +54,7 @@ public:
init = false;
sciters = 12;
Lfill = 1;
verbosity = 0;
}
bool Initialize()
......@@ -99,11 +100,13 @@ public:
iu.set_interval_end(moend);
ilu[mobeg] = 0;
ir[mobeg] = 0;
#if defined(REPORT_ILU)
std::cout << "Matrix overlap " << mobeg << ".." << moend << std::endl;
std::cout << "Local vector part " << vlocbeg << ".." << vlocend << std::endl;
std::cout << "Entire vector " << vbeg << ".." << vend << std::endl;
#endif
if( verbosity > 1 )
{
std::cout << "Matrix overlap " << mobeg << ".." << moend << std::endl;
std::cout << "Local vector part " << vlocbeg << ".." << vlocend << std::endl;
std::cout << "Entire vector " << vbeg << ".." << vend << std::endl;
}
//Rescale Matrix
DL.SetInterval(mobeg, moend);
info->PrepareVector(DR);
......@@ -143,15 +146,19 @@ public:
//for(k = mobeg; k != moend; k++) nza += A[k].Size();
for (k = mobeg; k != moend; k++)
{
#if defined(REPORT_ILU_PROGRESS)
if (k % 1000 == 0)
{
//std::cout << "precond: " << (double)(k-mobeg)/(double)(moend-mobeg)*100 << "\r";
//printf("%6.2f nza %12d nzl %12d nzu %12d nzu2 %12d\r", (double)(k-mobeg)/(double)(moend-mobeg)*100,nza,nzl,nzu,nzu2);
printf("precond: %6.2f\r", (double)(k - mobeg) / (double)(moend - mobeg) * 100);
fflush(stdout);
}
#endif
if( verbosity)
{
if (k % 100 == 0)
{
//std::cout << "precond: " << (double)(k-mobeg)/(double)(moend-mobeg)*100 << "\r";
if( verbosity )//== 1 )
printf("precond: %6.2f\r", (double)(k - mobeg) / (double)(moend - mobeg) * 100);
//else //statistics computation removed!
// printf("%6.2f nza %12d nzl %12d nzu %12d nzu2 %12d\r", (double)(k-mobeg)/(double)(moend-mobeg)*100,nza,nzl,nzu,nzu2);
fflush(stdout);
}
}
//Uncompress row
//row_uncompr
Sparse::Row &Ak = (*Alink)[k];
......@@ -425,24 +432,25 @@ public:
rit->second = rit->second / DL[k] / DR[rit->first];
//std::cout << "matisc: " << Timer() - timer << std::endl;
#if defined(REPORT_ILU)
INMOST_DATA_ENUM_TYPE nzu,nzl, nza;
nzu = 0;
nzl = 0;
nza = 0;
for(INMOST_DATA_ENUM_TYPE k = mobeg; k < moend; k++)
{
nzl += iu[k] - ilu[k];
nzu += ilu[k+1] - iu[k] - 1;
nza += (*Alink)[k].Size();
}
std::cout << " nonzeros in A = " << nza << std::endl;
std::cout << " nonzeros in L = " << nzl - (moend-mobeg) << std::endl;
std::cout << " nonzeros in U = " << nzu << std::endl;
std::cout << " nonzeros in LU = " << ilu[moend] - 1 << std::endl;
std::cout << " nonzeros in U2 = " << ir[moend] - 1 << std::endl;
if( verbosity > 1 )
{
INMOST_DATA_ENUM_TYPE nzu,nzl, nza;
nzu = 0;
nzl = 0;
nza = 0;
for(INMOST_DATA_ENUM_TYPE k = mobeg; k < moend; k++)
{
nzl += iu[k] - ilu[k];
nzu += ilu[k+1] - iu[k] - 1;
nza += (*Alink)[k].Size();
}
std::cout << " nonzeros in A = " << nza << std::endl;
std::cout << " nonzeros in L = " << nzl - (moend-mobeg) << std::endl;
std::cout << " nonzeros in U = " << nzu << std::endl;
std::cout << " nonzeros in LU = " << ilu[moend] - 1 << std::endl;
std::cout << " nonzeros in U2 = " << ir[moend] - 1 << std::endl;
}
//std::cout << __FUNCTION__ << " done" << std::endl;
#endif
/*
info->PrepareVector(div);
......@@ -488,6 +496,7 @@ public:
lui = b->lui;
iu = b->iu;
ilu = b->ilu;
verbosity = b->verbosity;
}
ILU2_preconditioner(const ILU2_preconditioner &other)
......
......@@ -20,6 +20,7 @@ namespace INMOST
pivot_condition = 1.0e+6;
pivot_diag = 1.0e+6;
fill_level = 3;
verbosity = 0;
}
SolverMLMPTILUC::SolverMLMPTILUC(const SolverInterface *other)
......@@ -41,6 +42,8 @@ namespace INMOST
solver->RealParameter(":tau2") = reuse_tolerance;
solver->EnumParameter(":scale_iters") = rescale_iterations;
solver->EnumParameter(":estimator") = condition_estimation;
solver->EnumParameter(":verbosity") = verbosity;
solver->EnumParameter("verbosity") = verbosity;
solver->RealParameter(":pivot_cond") = pivot_condition;
solver->RealParameter(":pivot_diag") = pivot_diag;
......@@ -67,6 +70,7 @@ namespace INMOST
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 if (name == "verbosity") verbosity = static_cast<INMOST_DATA_ENUM_TYPE>(atoi(val));
else SolverInner::SetParameter(name, value);
}
......
......@@ -10,7 +10,7 @@ namespace INMOST
class SolverMLMPTILUC : public SolverInner
{
INMOST_DATA_ENUM_TYPE rescale_iterations, condition_estimation, schwartz_overlap, gmres_substeps, reorder_nnz, fill_level;
INMOST_DATA_ENUM_TYPE rescale_iterations, condition_estimation, schwartz_overlap, gmres_substeps, reorder_nnz, fill_level, verbosity;
INMOST_DATA_REAL_TYPE drop_tolerance, reuse_tolerance, pivot_condition, pivot_diag;
public:
SolverMLMPTILUC();
......
......@@ -30,7 +30,7 @@ class MLMTILUC_preconditioner : public Method
INMOST_DATA_REAL_TYPE pivot_cond;
INMOST_DATA_REAL_TYPE pivot_diag;
INMOST_DATA_REAL_TYPE tau, eps;
INMOST_DATA_ENUM_TYPE sciters;
INMOST_DATA_ENUM_TYPE sciters, verbosity;
Sparse::Matrix * Alink;
Solver::OrderInfo * info;
bool init;
......
......@@ -10,7 +10,7 @@ namespace INMOST {
schwartz_overlap = 1;
gmres_substeps = 2;
reorder_nnz = 1;
verbosity = 0;
drop_tolerance = 0.005;
reuse_tolerance = 0.00005;
fill_level = 3;
......@@ -33,6 +33,8 @@ namespace INMOST {
solver->RealParameter(":tau") = drop_tolerance;
solver->RealParameter(":tau2") = reuse_tolerance;
solver->EnumParameter(":scale_iters") = rescale_iterations;
solver->EnumParameter(":verbosity") = verbosity;
solver->EnumParameter("verbosity") = verbosity;
if (sizeof(KSOLVER) == sizeof(BCGSL_solver)) {
solver->EnumParameter("levels") = gmres_substeps;
......@@ -58,6 +60,7 @@ namespace INMOST {
else if (name == "drop_tolerance") drop_tolerance = atof(val);
else if (name == "reuse_tolerance") reuse_tolerance = atof(val);
else if (name == "fill_level") fill_level = static_cast<INMOST_DATA_ENUM_TYPE>(atoi(val));
else if (name == "verbosity") verbosity = static_cast<INMOST_DATA_ENUM_TYPE>(atoi(val));
else SolverInner::SetParameter(name, value);
}
......
......@@ -9,7 +9,7 @@
namespace INMOST {
class SolverMPTILU2 : public SolverInner {
INMOST_DATA_ENUM_TYPE rescale_iterations, schwartz_overlap, gmres_substeps, reorder_nnz, fill_level;
INMOST_DATA_ENUM_TYPE rescale_iterations, schwartz_overlap, gmres_substeps, reorder_nnz, fill_level, verbosity;
INMOST_DATA_REAL_TYPE drop_tolerance, reuse_tolerance;
public:
SolverMPTILU2();
......
......@@ -111,6 +111,7 @@ void MTILU2_preconditioner::DumpMatrix(interval<INMOST_DATA_ENUM_TYPE, INMOST_DA
{
if (name == "fill") return Lfill;
else if (name == "scale_iters") return sciters;
else if (name == "verbosity" ) return verbosity;
throw -1;
}
MTILU2_preconditioner::MTILU2_preconditioner(Solver::OrderInfo & info)
......@@ -503,11 +504,13 @@ void MTILU2_preconditioner::DumpMatrix(interval<INMOST_DATA_ENUM_TYPE, INMOST_DA
iu.set_interval_end(moend);
ilu[mobeg] = 0;
ir[mobeg] = 0;
#if defined(REPORT_ILU)
std::cout << "Matrix overlap " << mobeg << ".." << moend << std::endl;
std::cout << "Local vector part " << vlocbeg << ".." << vlocend << std::endl;
std::cout << "Entire vector " << vbeg << ".." << vend << std::endl;
#endif
if( verbosity )
{
std::cout << "Matrix overlap " << mobeg << ".." << moend << std::endl;
std::cout << "Local vector part " << vlocbeg << ".." << vlocend << std::endl;
std::cout << "Entire vector " << vbeg << ".." << vend << std::endl;
}
#if defined(RESCALE_EQUALIZE_1NORM)
{
......@@ -590,15 +593,19 @@ void MTILU2_preconditioner::DumpMatrix(interval<INMOST_DATA_ENUM_TYPE, INMOST_DA
//for(k = mobeg; k != moend; k++) nza += A[k].Size();
for (k = mobeg; k != moend; k++)
{
#if defined(REPORT_ILU_PROGRESS)
if (k % 1000 == 0)
if( verbosity )
{
//std::cout << "precond: " << (double)(k-mobeg)/(double)(moend-mobeg)*100 << "\r";
//printf("%6.2f nza %12d nzl %12d nzu %12d nzu2 %12d\r", (double)(k-mobeg)/(double)(moend-mobeg)*100,nza,nzl,nzu,nzu2);
printf("precond: %6.2f\r", (double)(k - mobeg) / (double)(moend - mobeg) * 100);
fflush(stdout);
if (k % 100 == 0)
{
//std::cout << "precond: " << (double)(k-mobeg)/(double)(moend-mobeg)*100 << "\r";
if( verbosity )//== 1 )
printf("precond: %6.2f\r", (double)(k - mobeg) / (double)(moend - mobeg) * 100);
//else no statistics
// printf("%6.2f nza %12d nzl %12d nzu %12d nzu2 %12d\r", (double)(k-mobeg)/(double)(moend-mobeg)*100,nza,nzl,nzu,nzu2);
fflush(stdout);
}
}
#endif
//Uncompress row
//row_uncompr
sort_indeces.clear();
......@@ -855,30 +862,32 @@ void MTILU2_preconditioner::DumpMatrix(interval<INMOST_DATA_ENUM_TYPE, INMOST_DA
tfactor = Timer() - tfactor;
ttotal = Timer() - ttotal;
#if defined(REPORT_ILU_SUMMARY)
INMOST_DATA_ENUM_TYPE nzu,nzl;
nzu = 0;
nzl = 0;
for(INMOST_DATA_ENUM_TYPE k = mobeg; k < moend; k++)
if( verbosity > 1 )
{
nzl += iu[k] - ilu[k];
nzu += ilu[k+1] - iu[k] - 1;
}
std::cout << " nonzeros in A = " << nnz << std::endl;
std::cout << " nonzeros in L = " << nzl - (moend-mobeg) << std::endl;
std::cout << " nonzeros in U = " << nzu << std::endl;
std::cout << " nonzeros in LU = " << ilu[moend] - 1 << std::endl;
std::cout << " nonzeros in U2 = " << ir[moend] - 1 << std::endl;
//std::cout << __FUNCTION__ << " done" << std::endl;
printf("total %f\n",ttotal);
printf("reorder %f (%6.2f%%)\n", treorder, 100.0*treorder/ttotal);
INMOST_DATA_ENUM_TYPE nzu,nzl;
nzu = 0;
nzl = 0;
for(INMOST_DATA_ENUM_TYPE k = mobeg; k < moend; k++)
{
nzl += iu[k] - ilu[k];
nzu += ilu[k+1] - iu[k] - 1;
}
std::cout << " nonzeros in A = " << nnz << std::endl;
std::cout << " nonzeros in L = " << nzl - (moend-mobeg) << std::endl;
std::cout << " nonzeros in U = " << nzu << std::endl;
std::cout << " nonzeros in LU = " << ilu[moend] - 1 << std::endl;
std::cout << " nonzeros in U2 = " << ir[moend] - 1 << std::endl;
//std::cout << __FUNCTION__ << " done" << std::endl;
printf("total %f\n",ttotal);
printf("reorder %f (%6.2f%%)\n", treorder, 100.0*treorder/ttotal);
#if defined(REORDER_METIS_ND)
printf("metis graph %f nd %f\n", tmetisgraph, tmetisnd);
printf("metis graph %f nd %f\n", tmetisgraph, tmetisnd);
#endif
printf("rescale %f (%6.2f%%)\n", trescale, 100.0*trescale / ttotal);
printf("factor %f (%6.2f%%)\n", tfactor, 100.0*tfactor / ttotal);
printf("rescale %f (%6.2f%%)\n", trescale, 100.0*trescale / ttotal);
printf("factor %f (%6.2f%%)\n", tfactor, 100.0*tfactor / ttotal);
#endif
}
/*
//partition of unity for unrestricted additive schwartz
......@@ -945,14 +954,14 @@ void MTILU2_preconditioner::DumpMatrix(interval<INMOST_DATA_ENUM_TYPE, INMOST_DA
INMOST_DATA_ENUM_TYPE mobeg, moend, r, k, vbeg,vend; //, end;
info->GetOverlapRegion(info->GetRank(),mobeg,moend);
info->GetVectorRegion(vbeg,vend);
for(k = vbeg; k < mobeg; k++) temporary[k] = 0;
for(k = vbeg; k < mobeg; k++) temporary[k] = 0;
for(k = mobeg; k < moend; k++) temporary[k] = input[k];
for(k = moend; k < vend; k++) temporary[k] = 0;
for(k = moend; k < vend; k++) temporary[k] = 0;
//for(k = vbeg; k < vend; k++) temporary[k] = input[k];
//for(k = vbeg; k < vend; k++) temporary[k] = input[Perm[k]];
//for(k = vbeg; k < vend; k++) temporary[Perm[k]] = input[k];
//for(k = vbeg; k < vend; k++) temporary[k] = input[k];
//for(k = vbeg; k < vend; k++) temporary[k] = input[Perm[k]];
//for(k = vbeg; k < vend; k++) temporary[Perm[k]] = input[k];
for(k = mobeg; k < moend; k++) //iterate over L part
{
for(r = iu[k]-1; r > ilu[k]; r--)
......@@ -966,10 +975,10 @@ void MTILU2_preconditioner::DumpMatrix(interval<INMOST_DATA_ENUM_TYPE, INMOST_DA
temporary[k-1] *= luv[iu[k-1]];
}