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

Fix compilation issues with intel compiler and msvc compiler

parent b78dc3b8
......@@ -26,6 +26,19 @@ typedef Storage::real_array real_array;
typedef Storage::var_array var_array;
bool print_niter = false; //save file on nonlinear iterations
bool output_matrix = true;
bool norne_wells = true;
//data for wells
Cell cc[3] = {InvalidCell(),InvalidCell(),InvalidCell()};
real WI[3] = {50000,50000,50000};
real pbhp[3] = {265,105,110};
real ccnt[3][3] =
{
{4.567151e+05, 7.321079e+06, 2.767665e+03},
{4.609346e+05, 7.323503e+06, 2.597767e+03},
{4.595400e+05, 7.326078e+06, 2.803586e+03}
};
//#define OPTIMIZATION
......@@ -101,8 +114,8 @@ int main(int argc,char ** argv)
Tag tag_P; // Pressure
Tag tag_K; // Diffusion tensor
Tag tag_F; // Forcing term
Tag tag_BC; // Boundary conditions
Tag tag_W; // Gradient matrix acting on harmonic points on faces and returning gradient on faces
TagRealArray tag_BC; // Boundary conditions
......@@ -139,6 +152,18 @@ int main(int argc,char ** argv)
m->ExchangeData(tag_K,CELL,0); //Exchange diffusion tensor
}
if( norne_wells )
{
for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it) if( it->GetStatus() != Element::Ghost )
{
for(int q = 0; q < 3; ++q) if( it->Inside(ccnt[q]) )
{
cc[q] = it->self();
std::cout << "proc " << m->GetProcessorRank() << " found c" << q << " " << cc[q].LocalID() << std::endl;
}
}
}
if( m->HaveTag("PRESSURE") ) //Is there a pressure on the mesh?
tag_P = m->GetTag("PRESSURE"); //Get the pressure
......@@ -219,9 +244,9 @@ int main(int argc,char ** argv)
//~ W+=(rMatrix::Unit(NF) - R*(R.Transpose()*R).CholeskyInvert()*R.Transpose())*
//~ (2.0/(static_cast<real>(NF)*volume)*(NK*K.CholeskyInvert()*NK.Transpose()).Trace());
double ca = 0, ba = ca;
W = (NK+ca*L*R)*((NK+ca*L*R).Transpose()*R).Invert()*(NK+ca*L*R).Transpose();
W+= L - (1+ba)*(L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose();
W = (NK)*((NK).Transpose()*R).Invert()*(NK).Transpose();
W+= L - (L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose();
//~ W = Areas*W*Areas;
//access data structure for gradient matrix in mesh
......@@ -232,7 +257,7 @@ int main(int argc,char ** argv)
std::copy(W.data(),W.data()+NF*NF,store_W.data());
tag_WG[cell].resize(3*NF);
tag_WG(cell,3,NF) = (NK.Transpose()*R+ca*L*R).PseudoInvert(1.0e-12)*(NK+ca*L*R).Transpose();
tag_WG(cell,3,NF) = (NK.Transpose()*R).PseudoInvert(1.0e-12)*(NK).Transpose();
} //end of loop over cells
}
std::cout << "Construct W matrix: " << Timer() - ttt << std::endl;
......@@ -254,7 +279,32 @@ int main(int argc,char ** argv)
{ //Main loop for problem solution
Automatizator aut; // declare class to help manage unknowns
Automatizator::MakeCurrent(&aut);
dynamic_variable P(aut,aut.RegisterTag(tag_P,CELL|FACE)); //register pressure as primary unknown
MarkerType unk = m->CreateMarker();
for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
m->CellByLocalID(q).SetMarker(unk);
for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
{
Face face = m->FaceByLocalID(q);
face.SetMarker(unk);
if( face.Boundary() )
{
double alpha = 0, beta = 1, gamma = 0;
if( tag_BC.isValid() && face.HaveData(tag_BC) )
{
alpha = tag_BC[face][0];
beta = tag_BC[face][1];
gamma = tag_BC[face][2];
}
if( alpha != 0 && beta == 0 )
{
face.RemMarker(unk);
face.Real(tag_P) = gamma/alpha;
}
}
}
dynamic_variable P(aut,aut.RegisterTag(tag_P,CELL|FACE,unk)); //register pressure as primary unknown
aut.EnumerateEntries(); //enumerate all primary variables
std::cout << "Enumeration done, size " << aut.GetLastIndex() - aut.GetFirstIndex() << std::endl;
......@@ -272,10 +322,21 @@ int main(int argc,char ** argv)
for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
{
Face face = m->FaceByLocalID(q);
if( face.GetStatus() != Element::Ghost )
if( face.GetStatus() != Element::Ghost && face.GetMarker(unk) )
{
if( tag_BC.isValid() && face.HaveData(tag_BC) )
Text.SetAnnotation(P.Index(face),"Pressure guided by boundary condition");
if( face.Boundary() )
{
double alpha = 0, beta = 1, gamma = 0;
if( tag_BC.isValid() && face.HaveData(tag_BC) )
{
alpha = tag_BC[face][0];
beta = tag_BC[face][1];
gamma = tag_BC[face][2];
}
std::stringstream str;
str << "Pressure guided by boundary condition " << alpha << " " << beta << " " << gamma;
Text.SetAnnotation(P.Index(face),str.str());
}
else
Text.SetAnnotation(P.Index(face),"Interface pressure");
}
......@@ -311,7 +372,7 @@ int main(int argc,char ** argv)
FLUX.Resize(NF,1);
for(int k = 0; k < NF; ++k)
pF(k,0) = (P(faces[k]) - P(cell));
pF(k,0) = (P[faces[k]] - P[cell]);
FLUX = W*pF; //fluxes on faces
tag_PG(cell,3,1) = tag_WG(cell,3,NF)*pF;
if( cell.GetStatus() != Element::Ghost )
......@@ -322,6 +383,7 @@ int main(int argc,char ** argv)
for(int k = 0; k < NF; ++k) //loop over faces of current cell
{
if( faces[k].GetStatus() == Element::Ghost ) continue;
if( !faces[k].GetMarker(unk) ) continue;
int index = P.Index(faces[k]);
Locks.Lock(index);
......@@ -363,6 +425,12 @@ int main(int argc,char ** argv)
}
}
if( norne_wells )
{
for(int q = 0; q < 3; ++q) if( cc[q].isValid() )
R[P.Index(cc[q])] += WI[q]*(pbhp[q] - P[cc[q]])*cc[q].Volume();
}
//R[P.Index(m->BeginCell()->self())] = P[m->BeginCell()->self()];
}
std::cout << "assembled in " << Timer() - tttt << "\t\t\t" << std::endl;
......@@ -371,20 +439,53 @@ int main(int argc,char ** argv)
if( R.Norm() < 1.0e-4 ) break;
//~ R.GetJacobian().Save("A.mtx");
//~ R.GetResidual().Save("b.mtx");
//Solver S(Solver::INNER_ILU2);
//~ Solver S(Solver::INNER_MPTILUC);
Solver S(Solver::K3BIILU2);
Solver S(Solver::INNER_MPTILUC);
//~ Solver S(Solver::INNER_MLMPTILUC);
//~ Solver S(Solver::K3BIILU2);
//Solver S("superlu");
S.SetParameter("verbosity","1");
S.SetParameter("verbosity","3");
S.SetParameter("rescale_iterations", "10");
S.SetParameter("relative_tolerance", "1.0e-14");
S.SetParameter("absolute_tolerance", "1.0e-12");
S.SetParameter("drop_tolerance", "1.0e-3");
//~ S.SetParameter("drop_tolerance", "5.0e-2");
//~ S.SetParameter("reuse_tolerance", "2.5e-3");
S.SetParameter("drop_tolerance", "1.0e-2");
S.SetParameter("reuse_tolerance", "1.0e-4");
S.SetParameter("pivot_condition", "20");
double tset = Timer(), titr;
S.SetMatrix(R.GetJacobian());
tset = Timer() - tset;
if( output_matrix )
{
std::cout << "write A.mtx" << std::endl;
R.GetJacobian().Save("A.mtx");//,&Text);
std::cout << "write b.mtx" << std::endl;
R.GetResidual().Save("b.mtx");
std::cout << "write done, solve" << std::endl;
}
titr = Timer();
if( S.Solve(R.GetResidual(),Update) )
bool success = S.Solve(R.GetResidual(),Update);
titr = Timer() - titr;
if( success )
{
std::cout << "Solved system in " << S.Iterations() << " iterations, time for setup " << tset << "s iterations " << titr << "s" << std::endl;
if( output_matrix )
{
std::cout << "write x.mtx" << std::endl;
Update.Save("x.mtx");
//std::cout << "exit" << std::endl;
//exit(-1);
}
#if defined(USE_OMP)
#pragma omp parallel for
#endif
......@@ -392,6 +493,7 @@ int main(int argc,char ** argv)
{
Cell cell = m->CellByLocalID(q);
if( cell->GetStatus() == Element::Ghost ) continue;
if( !cell->GetMarker(unk) ) continue;
cell->Real(tag_P) -= Update[P.Index(cell)];
}
#if defined(USE_OMP)
......@@ -401,6 +503,7 @@ int main(int argc,char ** argv)
{
Face face = m->FaceByLocalID(q);
if( face->GetStatus() == Element::Ghost ) continue;
if( !face->GetMarker(unk) ) continue;
face->Real(tag_P) -= Update[P.Index(face)];
}
m->ExchangeData(tag_P, CELL|FACE, 0);
......
......@@ -582,12 +582,14 @@ namespace INMOST
template<typename typeB>
Matrix<typename Promote<Var,typeB>::type, pool_array_t<typename Promote<Var,typeB>::type> >
operator*(typeB coef) const;
#if defined(USE_AUTODIFF)
/// Multiply the matrix by a coefficient.
/// @param coef Coefficient.
/// @return Matrix multiplied by the coefficient.
template<class A>
Matrix<typename Promote<Var,variable>::type, pool_array_t<typename Promote<Var,variable>::type> >
operator*(shell_expression<A> const & coef) const {return operator*(variable(coef));}
#endif //USE_AUTODIFF
/// Multiply the matrix by the coefficient of the same type and store the result.
/// @param coef Coefficient.
/// @return Reference to the current matrix.
......@@ -599,12 +601,14 @@ namespace INMOST
template<typename typeB>
Matrix<typename Promote<Var,typeB>::type, pool_array_t<typename Promote<Var,typeB>::type> >
operator/(typeB coef) const;
#if defined(USE_AUTODIFF)
/// Divide the matrix by a coefficient of a different type.
/// @param coef Coefficient.
/// @return Matrix divided by the coefficient.
template<class A>
Matrix<typename Promote<Var,variable>::type, pool_array_t<typename Promote<Var,variable>::type> >
operator/(shell_expression<A> const & coef) const {return operator/(variable(coef));}
#endif //USE_AUTODIFF
/// Divide the matrix by the coefficient of the same type and store the result.
/// @param coef Coefficient.
/// @return Reference to the current matrix.
......@@ -2342,15 +2346,16 @@ namespace INMOST
assign((*this)(i,j),(*this)(i,j)+other(i,j));
return *this;
}
#if defined(USE_AUTODIFF)
template<>
template<>
__INLINE Matrix<typename Promote<INMOST_DATA_REAL_TYPE,variable>::type, pool_array_t<typename Promote<INMOST_DATA_REAL_TYPE,variable>::type> >
__INLINE Matrix<Promote<INMOST_DATA_REAL_TYPE,variable>::type, pool_array_t<Promote<INMOST_DATA_REAL_TYPE,variable>::type> >
AbstractMatrix<INMOST_DATA_REAL_TYPE>::operator*<variable>(const AbstractMatrix<variable> & other) const
{
//~ std::cout << __FUNCTION__ << std::endl;
assert(Cols() == other.Rows());
Matrix<typename Promote<INMOST_DATA_REAL_TYPE,variable>::type, pool_array_t<typename Promote<INMOST_DATA_REAL_TYPE,variable>::type> > ret(Rows(),other.Cols()); //check RVO
Matrix<Promote<INMOST_DATA_REAL_TYPE,variable>::type, pool_array_t<Promote<INMOST_DATA_REAL_TYPE,variable>::type> > ret(Rows(),other.Cols()); //check RVO
for(enumerator i = 0; i < Rows(); ++i) //loop rows
{
for(enumerator j = 0; j < other.Cols(); ++j) //loop columns
......@@ -2384,12 +2389,12 @@ namespace INMOST
template<>
template<>
__INLINE Matrix<typename Promote<variable,INMOST_DATA_REAL_TYPE>::type, pool_array_t<typename Promote<variable,INMOST_DATA_REAL_TYPE>::type> >
__INLINE Matrix<Promote<variable,INMOST_DATA_REAL_TYPE>::type, pool_array_t<Promote<variable,INMOST_DATA_REAL_TYPE>::type> >
AbstractMatrix<variable>::operator*<INMOST_DATA_REAL_TYPE>(const AbstractMatrix<INMOST_DATA_REAL_TYPE> & other) const
{
//~ std::cout << __FUNCTION__ << std::endl;
assert(Cols() == other.Rows());
Matrix<typename Promote<variable,INMOST_DATA_REAL_TYPE>::type, pool_array_t<typename Promote<variable,INMOST_DATA_REAL_TYPE>::type> > ret(Rows(),other.Cols()); //check RVO
Matrix<Promote<variable,INMOST_DATA_REAL_TYPE>::type, pool_array_t<Promote<variable,INMOST_DATA_REAL_TYPE>::type> > ret(Rows(),other.Cols()); //check RVO
for(enumerator i = 0; i < Rows(); ++i) //loop rows
{
for(enumerator j = 0; j < other.Cols(); ++j) //loop columns
......@@ -2422,12 +2427,12 @@ namespace INMOST
template<>
template<>
__INLINE Matrix<typename Promote<variable,variable>::type, pool_array_t<typename Promote<variable,variable>::type> >
__INLINE Matrix<Promote<variable,variable>::type, pool_array_t<Promote<variable,variable>::type> >
AbstractMatrix<variable>::operator*<variable>(const AbstractMatrix<variable> & other) const
{
//~ std::cout << __FUNCTION__ << std::endl;
assert(Cols() == other.Rows());
Matrix<typename Promote<variable,variable>::type, pool_array_t<typename Promote<variable,variable>::type> > ret(Rows(),other.Cols()); //check RVO
Matrix<Promote<variable,variable>::type, pool_array_t<Promote<variable,variable>::type> > ret(Rows(),other.Cols()); //check RVO
for(enumerator i = 0; i < Rows(); ++i) //loop rows
{
for(enumerator j = 0; j < other.Cols(); ++j) //loop columns
......@@ -2458,7 +2463,7 @@ namespace INMOST
}
return ret;
}
#endif //USE_AUTODIFF
template<typename Var>
template<typename typeB>
Matrix<typename Promote<Var,typeB>::type, pool_array_t<typename Promote<Var,typeB>::type> >
......@@ -2480,14 +2485,15 @@ namespace INMOST
return ret;
}
#if defined(USE_AUTODIFF)
template<>
template<>
__INLINE typename Promote<INMOST_DATA_REAL_TYPE,variable>::type
__INLINE Promote<INMOST_DATA_REAL_TYPE,variable>::type
AbstractMatrix<INMOST_DATA_REAL_TYPE>::DotProduct<variable>(const AbstractMatrix<variable> & other) const
{
assert(Cols() == other.Cols());
assert(Rows() == other.Rows());
typename Promote<INMOST_DATA_REAL_TYPE,variable>::type ret = 0.0;
Promote<INMOST_DATA_REAL_TYPE,variable>::type ret = 0.0;
if( CheckCurrentAutomatizator() )
{
Sparse::RowMerger & merger = GetCurrentMerger();
......@@ -2513,12 +2519,12 @@ namespace INMOST
template<>
template<>
__INLINE typename Promote<variable,INMOST_DATA_REAL_TYPE>::type
__INLINE Promote<variable,INMOST_DATA_REAL_TYPE>::type
AbstractMatrix<variable>::DotProduct<INMOST_DATA_REAL_TYPE>(const AbstractMatrix<INMOST_DATA_REAL_TYPE> & other) const
{
assert(Cols() == other.Cols());
assert(Rows() == other.Rows());
typename Promote<variable,INMOST_DATA_REAL_TYPE>::type ret = 0.0;
Promote<variable,INMOST_DATA_REAL_TYPE>::type ret = 0.0;
if( CheckCurrentAutomatizator() )
{
Sparse::RowMerger & merger = GetCurrentMerger();
......@@ -2544,12 +2550,12 @@ namespace INMOST
template<>
template<>
__INLINE typename Promote<variable,variable>::type
__INLINE Promote<variable,variable>::type
AbstractMatrix<variable>::DotProduct<variable>(const AbstractMatrix<variable> & other) const
{
assert(Cols() == other.Cols());
assert(Rows() == other.Rows());
typename Promote<variable,variable>::type ret = 0.0;
Promote<variable,variable>::type ret = 0.0;
if( CheckCurrentAutomatizator() )
{
Sparse::RowMerger & merger = GetCurrentMerger();
......@@ -2573,7 +2579,7 @@ namespace INMOST
}
return ret;
}
#endif //USE_AUTODIFF
......@@ -2814,10 +2820,10 @@ namespace INMOST
if( ierr ) *ierr = 0;
return ret;
}
#if defined(USE_AUTODIFF)
template<>
template<>
__INLINE Matrix<typename Promote<variable,variable>::type, pool_array_t<typename Promote<variable,variable>::type> >
__INLINE Matrix<Promote<variable,variable>::type, pool_array_t<Promote<variable,variable>::type> >
AbstractMatrix<variable>::CholeskySolve(const AbstractMatrix<variable> & B, int * ierr) const
{
const AbstractMatrix<variable> & A = *this;
......@@ -2825,7 +2831,7 @@ namespace INMOST
assert(A.Rows() == B.Rows());
enumerator n = A.Rows();
enumerator l = B.Cols();
Matrix<typename Promote<variable,variable>::type, pool_array_t<typename Promote<variable,variable>::type> > ret(B);
Matrix<Promote<variable,variable>::type, pool_array_t<Promote<variable,variable>::type> > ret(B);
SymmetricMatrix<variable, pool_array_t<variable> > L(A);
//Outer product
......@@ -2865,7 +2871,7 @@ namespace INMOST
}
}
// LY=B
Matrix<typename Promote<variable,variable>::type, pool_array_t<typename Promote<variable,variable>::type> > & Y = ret;
Matrix<Promote<variable,variable>::type, pool_array_t<Promote<variable,variable>::type> > & Y = ret;
for(enumerator i = 0; i < n; ++i)
{
for(enumerator k = 0; k < l; ++k)
......@@ -2894,7 +2900,7 @@ namespace INMOST
}
}
// L^TX = Y
Matrix<typename Promote<variable,variable>::type, pool_array_t<typename Promote<variable,variable>::type> > & X = ret;
Matrix<Promote<variable,variable>::type, pool_array_t<Promote<variable,variable>::type> > & X = ret;
for(enumerator it = n; it > 0; --it)
{
enumerator i = it-1;
......@@ -2934,7 +2940,7 @@ namespace INMOST
template<>
template<>
__INLINE Matrix<typename Promote<INMOST_DATA_REAL_TYPE,variable>::type, pool_array_t<typename Promote<INMOST_DATA_REAL_TYPE,variable>::type> >
__INLINE Matrix<Promote<INMOST_DATA_REAL_TYPE,variable>::type, pool_array_t<Promote<INMOST_DATA_REAL_TYPE,variable>::type> >
AbstractMatrix<INMOST_DATA_REAL_TYPE>::CholeskySolve(const AbstractMatrix<variable> & B, int * ierr) const
{
const AbstractMatrix<INMOST_DATA_REAL_TYPE> & A = *this;
......@@ -2942,7 +2948,7 @@ namespace INMOST
assert(A.Rows() == B.Rows());
enumerator n = A.Rows();
enumerator l = B.Cols();
Matrix<typename Promote<INMOST_DATA_REAL_TYPE,variable>::type, pool_array_t<typename Promote<INMOST_DATA_REAL_TYPE,variable>::type> > ret(B);
Matrix<Promote<INMOST_DATA_REAL_TYPE,variable>::type, pool_array_t<Promote<INMOST_DATA_REAL_TYPE,variable>::type> > ret(B);
SymmetricMatrix<INMOST_DATA_REAL_TYPE, pool_array_t<INMOST_DATA_REAL_TYPE> > L(A);
//Outer product
......@@ -2982,7 +2988,7 @@ namespace INMOST
}
}
// LY=B
Matrix<typename Promote<variable,variable>::type, pool_array_t<typename Promote<variable,variable>::type> > & Y = ret;
Matrix<Promote<variable,variable>::type, pool_array_t<Promote<variable,variable>::type> > & Y = ret;
for(enumerator i = 0; i < n; ++i)
{
for(enumerator k = 0; k < l; ++k)
......@@ -3010,7 +3016,7 @@ namespace INMOST
}
}
// L^TX = Y
Matrix<typename Promote<variable,variable>::type, pool_array_t<typename Promote<variable,variable>::type> > & X = ret;
Matrix<Promote<variable,variable>::type, pool_array_t<Promote<variable,variable>::type> > & X = ret;
for(enumerator it = n; it > 0; --it)
{
enumerator i = it-1;
......@@ -3048,7 +3054,7 @@ namespace INMOST
template<>
template<>
__INLINE Matrix<typename Promote<variable,INMOST_DATA_REAL_TYPE>::type, pool_array_t<typename Promote<variable,INMOST_DATA_REAL_TYPE>::type> >
__INLINE Matrix<Promote<variable,INMOST_DATA_REAL_TYPE>::type, pool_array_t<Promote<variable,INMOST_DATA_REAL_TYPE>::type> >
AbstractMatrix<variable>::CholeskySolve(const AbstractMatrix<INMOST_DATA_REAL_TYPE> & B, int * ierr) const
{
const AbstractMatrix<variable> & A = *this;
......@@ -3056,7 +3062,7 @@ namespace INMOST
assert(A.Rows() == B.Rows());
enumerator n = A.Rows();
enumerator l = B.Cols();
Matrix<typename Promote<variable,INMOST_DATA_REAL_TYPE>::type, pool_array_t<typename Promote<variable,INMOST_DATA_REAL_TYPE>::type> > ret(B);
Matrix<Promote<variable,INMOST_DATA_REAL_TYPE>::type, pool_array_t<Promote<variable,INMOST_DATA_REAL_TYPE>::type> > ret(B);
SymmetricMatrix<variable, pool_array_t<variable> > L(A);
//Outer product
......@@ -3096,7 +3102,7 @@ namespace INMOST
}
}
// LY=B
Matrix<typename Promote<variable,INMOST_DATA_REAL_TYPE>::type, pool_array_t<typename Promote<variable,INMOST_DATA_REAL_TYPE>::type> > & Y = ret;
Matrix<Promote<variable,INMOST_DATA_REAL_TYPE>::type, pool_array_t<Promote<variable,INMOST_DATA_REAL_TYPE>::type> > & Y = ret;
for(enumerator i = 0; i < n; ++i)
{
for(enumerator k = 0; k < l; ++k)
......@@ -3125,7 +3131,7 @@ namespace INMOST
}
}
// L^TX = Y
Matrix<typename Promote<variable,INMOST_DATA_REAL_TYPE>::type, pool_array_t<typename Promote<variable,INMOST_DATA_REAL_TYPE>::type> > & X = ret;
Matrix<Promote<variable,INMOST_DATA_REAL_TYPE>::type, pool_array_t<Promote<variable,INMOST_DATA_REAL_TYPE>::type> > & X = ret;
for(enumerator it = n; it > 0; --it)
{
enumerator i = it-1;
......@@ -3161,7 +3167,7 @@ namespace INMOST
if( ierr ) *ierr = 0;
return ret;
}
#endif //USE_AUTODIFF
template<typename Var>
template<typename typeB>
Matrix<typename Promote<Var,typeB>::type, pool_array_t<typename Promote<Var,typeB>::type> >
......
......@@ -97,6 +97,11 @@ namespace INMOST
heap[size+1] = ENUMUNDEF;
return min;
}
INMOST_DATA_ENUM_TYPE PeekHeap()
{
if( size == 0 ) return ENUMUNDEF;
return heap[1];
}
void DecreaseKey(INMOST_DATA_ENUM_TYPE i, INMOST_DATA_REAL_TYPE key)
{
keys[i] = key;
......@@ -163,7 +168,7 @@ namespace INMOST
for(INMOST_DATA_ENUM_TYPE i = 0; i <= size_max; i++)
index[i] = ENUMUNDEF;
}
void PushHeap(INMOST_DATA_ENUM_TYPE i, INMOST_DATA_REAL_TYPE key)
void PushHeap(INMOST_DATA_ENUM_TYPE i, KeyType key)
{
size++;
index[i] = size;
......@@ -186,17 +191,17 @@ namespace INMOST
if( size == 0 ) return ENUMUNDEF;
return heap[1];
}
void DecreaseKey(INMOST_DATA_ENUM_TYPE i, INMOST_DATA_REAL_TYPE key)
void DecreaseKey(INMOST_DATA_ENUM_TYPE i, KeyType key)
{
keys[i] = key;
bubbleUp(index[i]);
}
void IncreaseKey(INMOST_DATA_ENUM_TYPE i, INMOST_DATA_REAL_TYPE key)
void IncreaseKey(INMOST_DATA_ENUM_TYPE i, KeyType key)
{
keys[i] = key;
bubbleDown(index[i]);
}
void ChangeKey(INMOST_DATA_ENUM_TYPE i, INMOST_DATA_REAL_TYPE key)
void ChangeKey(INMOST_DATA_ENUM_TYPE i, KeyType key)
{
keys[i] = key;
if( !Comparator()(keys[i],key) )
......
......@@ -1284,7 +1284,7 @@ namespace INMOST
std::vector< int > points_cluster(total_points,-1);
std::vector< double > cluster_center(K*3,0.0);
std::vector< int > cluster_npoints(K,0);
std::vector< double > cluster_weight(K,1/(double)K);
//std::vector< double > cluster_weight(K,1/(double)K);
std::vector< double > cluster_shift(K*3,0);
#if defined(USE_MPI)
std::vector< double > cluster_center_tmp(K*3);
......@@ -1503,7 +1503,7 @@ namespace INMOST
v[2] = (points_center[i*3+2] - cluster_center[j*3+2]);
//the bigger is the cluster weight the further is the distance
double l = (v[0]*v[0] + v[1]*v[1] + v[2]*v[2]) + cluster_weight[j]*mesh_dist;
double l = (v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);// + cluster_weight[j]*mesh_dist;
if( l < lmin )
{
lmin = l;
......@@ -1579,7 +1579,7 @@ namespace INMOST
cluster_center[i*3+1] /= (double) cluster_npoints[i];
cluster_center[i*3+2] /= (double) cluster_npoints[i];
//recompute cluster weights based on the number of points each cluster possess
cluster_weight[i] = (cluster_weight[i]*0.25+cluster_npoints[i]/(double)total_global_points*0.75);
//cluster_weight[i] = (cluster_weight[i]*0.25+cluster_npoints[i]/(double)total_global_points*0.75);
//if( m->GetProcessorRank() == 0 ) std::cout << cluster_weight[i]*K << " (" <<cluster_npoints[i]<< ") ";
}
//if( m->GetProcessorRank() == 0 ) std::cout << std::endl;
......
......@@ -24,7 +24,7 @@ using namespace INMOST;
//~ #define REORDER_RCM
#define REORDER_WRCM
// #define REORDER_BRCM
//~ #define REORDER_BRCM
//#define REORDER_NNZ
#if defined(USE_SOLVER_METIS)
//~ #define REORDER_METIS_ND
......@@ -42,29 +42,29 @@ const INMOST_DATA_ENUM_TYPE UNDEF = ENUMUNDEF, EOL = ENUMUNDEF - 1;
//~ #define PREMATURE_DROPPING
//#define EQUALIZE_1NORM
//~ #define EQUALIZE_1NORM
//~ #define EQUALIZE_2NORM
#define EQUALIZE_IDOMINANCE
#define PIVOT_THRESHOLD
#define PIVOT_THRESHOLD_VALUE 1.0e-12
#define PIVOT_THRESHOLD_VALUE 1.0e-9
//#define DIAGONAL_PERTURBATION
#define DIAGONAL_PERTURBATION_REL 1.0e-7
#define DIAGONAL_PERTURBATION_ABS 1.0e-10
#define ILUC2
#define ILUC2_SCHUR
//~ #define ILUC2_SCHUR
//#define PIVOT_COND_DEFAULT 0.1/tau
#define PIVOT_COND_DEFAULT 1.0e+2
#define PIVOT_DIAG_DEFAULT 1.0e+5
//~ #define SCHUR_DROPPING_LF
//~ #define SCHUR_DROPPING_EU
//~ #define SCHUR_DROPPING_S
#define SCHUR_DROPPING_S
#define SCHUR_DROPPING_LF_PREMATURE