Commit 635cc476 authored by Kirill Terekhov's avatar Kirill Terekhov

Synchronize

Cleanup of Row class and Matrix class.

Significantly reduced memory required for Row class and consequently for
multivar_expression class. Moved out openmp locks and row annotation as
additional service classes.

Introduced classes HessianRow and HessianMatrix for future possibility
of hessian computation in automatic differentiation.
parent e3333185
......@@ -134,6 +134,7 @@ int main(int argc,char ** argv)
Solver S(Solver::INNER_ILU2); // Specify the linear solver to ASM+ILU2+BiCGStab one
S.SetParameterReal("absolute_tolerance",1e-8);
Residual R; // Residual vector
Sparse::LockService Locks;
Sparse::Vector Update; // Declare the solution and the right-hand side vectors
Mesh::GeomParam table;
......@@ -155,6 +156,7 @@ int main(int argc,char ** argv)
// Set the indeces intervals for the matrix and vectors
R.SetInterval(aut.GetFirstIndex(),aut.GetLastIndex());
Locks.SetInterval(aut.GetFirstIndex(),aut.GetLastIndex());
Update.SetInterval(aut.GetFirstIndex(),aut.GetLastIndex());
//~ std::cout << m->GetProcessorRank() << " A,x,b interval " << idmin << ":" << idmax << " size " << idmax-idmin << std::endl;
dynamic_variable Phi(aut,iphi);
......@@ -193,9 +195,9 @@ int main(int argc,char ** argv)
bnd_pnt[2] = r1_cnt[2] + dist * f_nrm[2];
T = r1->Real(tensor_K) * f_area / dist;
//flux = T * (func(bnd_pnt,0) - variable(aut,r1,iphi));
R[i1].Lock();
Locks.Lock(i1);
R[i1] -= T * (func(bnd_pnt,0) - Phi(r1));
R[i1].Unlock();
Locks.UnLock(i1);
}
else
{
......@@ -209,15 +211,15 @@ int main(int argc,char ** argv)
flux = T * (Phi(r2) - Phi(r1));
if( s1 != Element::Ghost )
{
R[i1].Lock();
Locks.Lock(i1);
R[i1] -= flux;
R[i1].Unlock();
Locks.UnLock(i1);
}
if( s2 != Element::Ghost )
{
R[i2].Lock();
Locks.Lock(i2);
R[i2] += flux;
R[i2].Unlock();
Locks.UnLock(i2);
}
}
}
......
......@@ -37,6 +37,11 @@ int main(int argc,char ** argv)
#endif
if( argc > 1 )
{
std::cout << "Row: " << sizeof(Sparse::Row) << std::endl;
std::cout << "HessianRow: " << sizeof(Sparse::HessianRow) << std::endl;
std::cout << "Matrix: " << sizeof(Sparse::Matrix) << std::endl;
std::cout << "HessianMatrix: " << sizeof(Sparse::HessianMatrix) << std::endl;
std::cout << "variable: " << sizeof(multivar_expression) << std::endl;
double ttt; // Variable used to measure timing
bool repartition = false; // Is it required to redistribute the mesh?
Mesh * m = new Mesh(); // Create an empty mesh
......@@ -365,21 +370,22 @@ int main(int argc,char ** argv)
Residual R("",aut.GetFirstIndex(),aut.GetLastIndex());
Sparse::LockService Locks(aut.GetFirstIndex(),aut.GetLastIndex());
Sparse::AnnotationService Text(aut.GetFirstIndex(),aut.GetLastIndex());
Sparse::Vector Update ("",aut.GetFirstIndex(),aut.GetLastIndex()); //vector for update
{//Annotate matrix
for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
{
Cell cell = m->CellByLocalID(q);
R.GetJacobian().Annotation(P.Index(cell)) = "Cell-centered pressure value";
Text.SetAnnotation(P.Index(cell),"Cell-centered pressure value");
}
for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
{
Face face = m->FaceByLocalID(q);
if( tag_BC.isValid() && face.HaveData(tag_BC) )
R.GetJacobian().Annotation(P.Index(face)) = "Pressure guided by boundary condition";
Text.SetAnnotation(P.Index(face),"Pressure guided by boundary condition");
else
R.GetJacobian().Annotation(P.Index(face)) = "Interface pressure";
Text.SetAnnotation(P.Index(face),"Interface pressure");
}
}
......@@ -388,6 +394,7 @@ int main(int argc,char ** argv)
do
{
R.Clear(); //clean up the residual
double tttt = Timer();
#if defined(USE_OMP)
#pragma omp parallel for
#endif
......@@ -407,7 +414,7 @@ int main(int argc,char ** argv)
for(int k = 0; k < NF; ++k) //loop over faces of current cell
{
int index = P.Index(faces[k]);
R[index].Lock();
Locks.Lock(index);
if( tag_BC.isValid() && faces[k].HaveData(tag_BC) )
{
real_array BC = faces[k].RealArray(tag_BC);
......@@ -415,7 +422,7 @@ int main(int argc,char ** argv)
}
else
R[index] -= FLUX(k,0);
R[index].Unlock();
Locks.UnLock(index);
}
} //end of loop over cells
......@@ -432,9 +439,11 @@ int main(int argc,char ** argv)
}
}
std::cout << "assembled in " << Timer() - tttt << "\t\t\t" << std::endl;
R.Rescale();
R.GetJacobian().Save("jacobian.mtx");
R.GetJacobian().Save("jacobian.mtx",&Text);
R.GetResidual().Save("residual.mtx");
......
......@@ -118,6 +118,7 @@ public:
class Reorder_ARMS
{
interval<INMOST_DATA_ENUM_TYPE, bool> markers;
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> reorder;
interval<INMOST_DATA_ENUM_TYPE, std::vector< INMOST_DATA_ENUM_TYPE> > neighbours;
std::vector< std::pair<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> > node_neighbours;
......@@ -129,7 +130,7 @@ class Reorder_ARMS
Sparse::Matrix & A;
void add_block_entry(INMOST_DATA_ENUM_TYPE ind)
{
A[ind].SetMarker();
markers[ind] = true;
block_entries.push_back(ind);
reorder[ind] = visited++;
stack.push_back(ind);
......@@ -139,7 +140,7 @@ class Reorder_ARMS
for (INMOST_DATA_ENUM_TYPE j = 0; j < block_entries.size(); j++)
{
for (INMOST_DATA_ENUM_TYPE k = 0; k < neighbours[block_entries[j]].size(); ++k)
A[neighbours[block_entries[j]][k]].SetMarker();
markers[neighbours[block_entries[j]][k]] = true;
}
sblock += static_cast<int>(block_entries.size());
block_ends.push_back(sblock);
......@@ -172,7 +173,7 @@ class Reorder_ARMS
return static_cast<Storage::enumerator>(neighbours[ind].size());
INMOST_DATA_ENUM_TYPE ret = 0, k;
for (k = 0; k < neighbours[ind].size(); k++)
if (!A[neighbours[ind][k]].GetMarker())
if (!markers[neighbours[ind][k]])
ret++;
return ret;
}
......@@ -204,6 +205,8 @@ public:
reorder.set_interval_end(mend);
neighbours.set_interval_beg(mbeg);
neighbours.set_interval_end(mend);
markers.set_interval_beg(mbeg);
markers.set_interval_end(mend);
for (INMOST_DATA_ENUM_TYPE k = mbeg; k < mend; ++k)
{
......@@ -243,7 +246,7 @@ public:
select = ENUMUNDEF;
for (INMOST_DATA_ENUM_TYPE k = mbeg; k < mend; ++k)
{
if (!A[k].GetMarker())
if (!markers[k])
{
if (select == ENUMUNDEF || GetOrder(k) < GetOrder(select))
select = k;
......@@ -266,7 +269,7 @@ public:
for (INMOST_DATA_ENUM_TYPE k = 0; k != neighbours[select].size(); ++k)
{
INMOST_DATA_ENUM_TYPE j = neighbours[select][k];
if (!A[j].GetMarker())
if (!markers[j])
{
node_neighbours.push_back(std::pair<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE>(GetOrder(j), j));
//A[j].SetMarker();
......@@ -289,7 +292,7 @@ public:
std::cout << "Left untouched: " << mend - mbeg - visited << std::endl;
for (INMOST_DATA_ENUM_TYPE k = mbeg; k < mend; ++k)
{
A[k].RemMarker();
markers[k] = false;
if (reorder[k] == ENUMUNDEF)
reorder[k] = visited++;
}
......
......@@ -130,6 +130,7 @@ int main(int argc,char ** argv)
ttt = Timer();
Solver S(Solver::INNER_ILU2); // Specify the linear solver to ASM+ILU2+BiCGStab one
S.SetParameterReal("absolute_tolerance",1e-8);
Sparse::LockService L;
Sparse::Matrix A; // Declare the matrix of the linear system to be solved
Sparse::Vector x,b; // Declare the solution and the right-hand side vectors
......@@ -155,6 +156,7 @@ int main(int argc,char ** argv)
}
// Set the indeces intervals for the matrix and vectors
L.SetInterval(idmin,idmax);
A.SetInterval(idmin,idmax);
x.SetInterval(idmin,idmax);
b.SetInterval(idmin,idmax);
......@@ -192,10 +194,10 @@ int main(int argc,char ** argv)
bnd_pnt[1] = r1_cnt[1] + dist * f_nrm[1];
bnd_pnt[2] = r1_cnt[2] + dist * f_nrm[2];
Coef = K1 * f_area / dist;
A[id1].Lock();
L.Lock(id1);
A[id1][id1] += -Coef;
b[id1] += -Coef * func(bnd_pnt, 0);
A[id1].Unlock();
L.UnLock(id1);
}
else
{
......@@ -209,17 +211,17 @@ int main(int argc,char ** argv)
if( s1 != Element::Ghost )
{
A[id1].Lock();
L.Lock(id1);
A[id1][id1] += -Coef;
A[id1][id2] += Coef;
A[id1].Unlock();
L.UnLock(id2);
}
if( s2 != Element::Ghost )
{
A[id2].Lock();
L.Lock(id2);
A[id2][id1] += Coef;
A[id2][id2] += -Coef;
A[id2].Unlock();
L.UnLock(id2);
}
}
}
......
......@@ -4,6 +4,7 @@
#include "inmost_common.h"
#include "inmost_mesh.h"
#include "inmost_solver.h"
#include "inmost_variable.h"
#include <sstream> //for debug
//#define NEW_VERSION
......@@ -83,6 +84,94 @@ namespace INMOST
class expr;
class Residual
{
Sparse::Matrix jacobian;
Sparse::Vector residual;
public:
Residual(std::string name = "", INMOST_DATA_ENUM_TYPE start = 0, INMOST_DATA_ENUM_TYPE end = 0, INMOST_MPI_Comm _comm = INMOST_MPI_COMM_WORLD)
: jacobian(name,start,end,_comm),residual(name,start,end,_comm) {}
Residual(const Residual & other)
: jacobian(other.jacobian), residual(other.residual)
{}
Residual & operator =(Residual const & other)
{
jacobian = other.jacobian;
residual = other.residual;
return *this;
}
Sparse::Matrix & GetJacobian() {return jacobian;}
const Sparse::Matrix & GetJacobian() const {return jacobian;}
Sparse::Vector & GetResidual() {return residual;}
const Sparse::Vector & GetResidual() const {return residual;}
INMOST_DATA_ENUM_TYPE GetFirstIndex() const {return residual.GetFirstIndex();}
INMOST_DATA_ENUM_TYPE GetLastIndex() const {return residual.GetLastIndex();}
void GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const
{
start = residual.GetFirstIndex();
end = residual.GetLastIndex();
}
void SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end)
{
jacobian.SetInterval(beg,end);
residual.SetInterval(beg,end);
}
multivar_expression_reference operator [](INMOST_DATA_ENUM_TYPE row)
{
return multivar_expression_reference(residual[row],&jacobian[row]);
}
void ClearResidual()
{
for(Sparse::Vector::iterator it = residual.Begin(); it != residual.End(); ++it) (*it) = 0.0;
}
void ClearJacobian()
{
for(Sparse::Matrix::iterator it = jacobian.Begin(); it != jacobian.End(); ++it)
it->Clear();
}
void Clear()
{
#if defined(USE_OMP)
#pragma omp for
#endif
for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k)
{
residual[k] = 0.0;
jacobian[k].Clear();
}
}
INMOST_DATA_REAL_TYPE Norm()
{
INMOST_DATA_REAL_TYPE ret = 0;
#if defined(USE_OMP)
#pragma omp for
#endif
for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k)
ret += residual[k]*residual[k];
return sqrt(ret);
}
/// Normalize entries in jacobian and right hand side
void Rescale()
{
#if defined(USE_OMP)
#pragma omp for
#endif
for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k)
{
INMOST_DATA_REAL_TYPE norm = 0.0;
for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q)
norm += jacobian[k].GetValue(q)*jacobian[k].GetValue(q);
norm = sqrt(norm);
if( norm )
{
norm = 1.0/norm;
residual[k] *= norm;
for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q)
jacobian[k].GetValue(q) *= norm;
}
}
}
};
class Automatizator
{
......
......@@ -110,6 +110,7 @@ namespace INMOST
{
INMOST_DATA_REAL_TYPE value;
Sparse::Row entries;
Sparse::HessianRow hessian_entries;
public:
multivar_expression() :value(0) {}
multivar_expression(INMOST_DATA_REAL_TYPE pvalue) : value(pvalue) {}
......@@ -311,12 +312,11 @@ namespace INMOST
class multivar_expression_reference : public shell_expression<multivar_expression_reference>
{
INMOST_DATA_REAL_TYPE & value;
Sparse::Row & entries;
Sparse::Row * entries;
Sparse::HessianRow * hessian_entries;
public:
void Lock() {entries.Lock();}
void Unlock() {entries.Unlock();}
/// Constructor, set links to the provided value and entries
multivar_expression_reference(INMOST_DATA_REAL_TYPE & _value, Sparse::Row & _entries)
multivar_expression_reference(INMOST_DATA_REAL_TYPE & _value, Sparse::Row * _entries)
: value(_value), entries(_entries) {}
/// Copy constructor, sets links to the same reference of value and entries
multivar_expression_reference(const multivar_expression_reference & other)
......@@ -328,7 +328,7 @@ namespace INMOST
/// Retrive derivatives with multiplier into Sparse::RowMerger structure.
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
{
for(Sparse::Row::iterator it = entries.Begin(); it != entries.End(); ++it)
for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it)
r[it->first] += it->second*mult;
}
/// Retrive derivatives with multiplier into Sparse::Row structure.
......@@ -338,53 +338,53 @@ namespace INMOST
FromGetJacobian(*this,mult,r);
else
{
for(Sparse::Row::iterator it = entries.Begin(); it != entries.End(); ++it)
for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it)
r[it->first] += it->second*mult;
}
}
__INLINE multivar_expression_reference & operator = (INMOST_DATA_REAL_TYPE pvalue)
{
value = pvalue;
entries.Clear();
entries->Clear();
return *this;
}
__INLINE multivar_expression_reference & operator = (basic_expression const & expr)
{
value = expr.GetValue();
if( CheckCurrentAutomatizator() )
FromBasicExpression(entries,expr);
FromBasicExpression(*entries,expr);
else
{
Sparse::Row tmp;
expr.GetJacobian(1.0,tmp);
entries.Swap(tmp);
entries->Swap(tmp);
}
return *this;
}
__INLINE multivar_expression_reference & operator = (multivar_expression_reference const & other)
{
value = other.GetValue();
entries = other.GetRow();
*entries = other.GetRow();
return *this;
}
__INLINE multivar_expression_reference & operator = (multivar_expression const & other)
{
value = other.GetValue();
entries = other.GetRow();
*entries = other.GetRow();
return *this;
}
__INLINE Sparse::Row & GetRow() {return entries;}
__INLINE const Sparse::Row & GetRow() const {return entries;}
__INLINE Sparse::Row & GetRow() {return *entries;}
__INLINE const Sparse::Row & GetRow() const {return *entries;}
__INLINE multivar_expression_reference & operator +=(basic_expression const & expr)
{
value += expr.GetValue();
if( CheckCurrentAutomatizator() )
AddBasicExpression(entries,1.0,1.0,expr);
AddBasicExpression(*entries,1.0,1.0,expr);
else
{
Sparse::Row tmp(entries);
Sparse::Row tmp(*entries);
expr.GetJacobian(1.0,tmp);
entries.Swap(tmp);
entries->Swap(tmp);
}
return *this;
}
......@@ -392,12 +392,12 @@ namespace INMOST
{
value -= expr.GetValue();
if( CheckCurrentAutomatizator() )
AddBasicExpression(entries,1.0,-1.0,expr);
AddBasicExpression(*entries,1.0,-1.0,expr);
else
{
Sparse::Row tmp(entries);
Sparse::Row tmp(*entries);
expr.GetJacobian(-1.0,tmp);
entries.Swap(tmp);
entries->Swap(tmp);
}
return *this;
}
......@@ -405,13 +405,13 @@ namespace INMOST
{
INMOST_DATA_REAL_TYPE lval = value, rval = expr.GetValue();
if( CheckCurrentAutomatizator() )
AddBasicExpression(entries,rval,lval,expr);
AddBasicExpression(*entries,rval,lval,expr);
else
{
Sparse::Row tmp(entries);
Sparse::Row tmp(*entries);
for(Sparse::Row::iterator it = tmp.Begin(); it != tmp.End(); ++it) it->second *= rval;
expr.GetJacobian(lval,tmp);
entries.Swap(tmp);
entries->Swap(tmp);
}
value *= rval;
return *this;
......@@ -422,13 +422,13 @@ namespace INMOST
INMOST_DATA_REAL_TYPE reciprocial_rval = 1.0/rval;
value *= reciprocial_rval;
if( CheckCurrentAutomatizator() )
AddBasicExpression(entries,reciprocial_rval,-value*reciprocial_rval,expr);
AddBasicExpression(*entries,reciprocial_rval,-value*reciprocial_rval,expr);
else
{
Sparse::Row tmp(entries);
Sparse::Row tmp(*entries);
for(Sparse::Row::iterator it = tmp.Begin(); it != tmp.End(); ++it) it->second *= reciprocial_rval;
expr.GetJacobian(-value*reciprocial_rval,tmp);
entries.Swap(tmp);
entries->Swap(tmp);
}
return *this;
}
......@@ -445,19 +445,19 @@ namespace INMOST
__INLINE multivar_expression_reference & operator *=(INMOST_DATA_REAL_TYPE right)
{
value *= right;
for(Sparse::Row::iterator it = entries.Begin(); it != entries.End(); ++it) it->second *= right;
for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it) it->second *= right;
return *this;
}
__INLINE multivar_expression_reference & operator /=(INMOST_DATA_REAL_TYPE right)
{
value /= right;
for(Sparse::Row::iterator it = entries.Begin(); it != entries.End(); ++it) it->second /= right;
for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it) it->second /= right;
return *this;
}
bool check_nans() const
{
if( value != value ) return true;
for(Sparse::Row::iterator it = entries.Begin(); it != entries.End(); ++it)
for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it)
if( it->second != it->second ) return true;
return false;
}
......@@ -1269,95 +1269,6 @@ namespace INMOST
typedef multivar_expression variable;
typedef var_expression unknown;
class Residual
{
Sparse::Matrix jacobian;
Sparse::Vector residual;
public:
Residual(std::string name = "", INMOST_DATA_ENUM_TYPE start = 0, INMOST_DATA_ENUM_TYPE end = 0, INMOST_MPI_Comm _comm = INMOST_MPI_COMM_WORLD)
: jacobian(name,start,end,_comm),residual(name,start,end,_comm) {}
Residual(const Residual & other)
: jacobian(other.jacobian), residual(other.residual)
{}
Residual & operator =(Residual const & other)
{
jacobian = other.jacobian;
residual = other.residual;
return *this;
}
Sparse::Matrix & GetJacobian() {return jacobian;}
const Sparse::Matrix & GetJacobian() const {return jacobian;}
Sparse::Vector & GetResidual() {return residual;}
const Sparse::Vector & GetResidual() const {return residual;}
INMOST_DATA_ENUM_TYPE GetFirstIndex() const {return residual.GetFirstIndex();}
INMOST_DATA_ENUM_TYPE GetLastIndex() const {return residual.GetLastIndex();}
void GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const
{
start = residual.GetFirstIndex();
end = residual.GetLastIndex();
}
void SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end)
{
jacobian.SetInterval(beg,end);
residual.SetInterval(beg,end);
}
multivar_expression_reference operator [](INMOST_DATA_ENUM_TYPE row)
{
return multivar_expression_reference(residual[row],jacobian[row]);
}
void ClearResidual()
{
for(Sparse::Vector::iterator it = residual.Begin(); it != residual.End(); ++it) (*it) = 0.0;
}
void ClearJacobian()
{
for(Sparse::Matrix::iterator it = jacobian.Begin(); it != jacobian.End(); ++it)
it->Clear();
}
void Clear()
{
#if defined(USE_OMP)
#pragma omp for
#endif
for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k)
{
residual[k] = 0.0;
jacobian[k].Clear();
}
}
INMOST_DATA_REAL_TYPE Norm()
{
INMOST_DATA_REAL_TYPE ret = 0;
#if defined(USE_OMP)
#pragma omp for
#endif
for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k)
ret += residual[k]*residual[k];
return sqrt(ret);
}
/// Normalize entries in jacobian and right hand side
void Rescale()
{
#if defined(USE_OMP)
#pragma omp for
#endif
for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k)
{
INMOST_DATA_REAL_TYPE norm = 0.0;
for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q)
norm += jacobian[k].GetValue(q)*jacobian[k].GetValue(q);
norm = sqrt(norm);
if( norm )
{
norm = 1.0/norm;
residual[k] *= norm;
for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q)
jacobian[k].GetValue(q) *= norm;
}
}
}
};
}
......
......@@ -253,12 +253,12 @@ namespace INMOST
void SetParameterReal(std::string name, INMOST_DATA_REAL_TYPE value);
/// Get the used defined name of the Solver.
std::string GetName() {return name;}
/// Get the package Type.
Type GetPackage() const {return _pack;}
/// Set the matrix and construct the preconditioner.
/// @param A Matrix A in linear problem Ax = b
/// @param ModifiedPattern Indicates whether the structure of the matrix have
/// changed since last call to Solver::SetMatrix.
/// @param OldPreconditioner If this parameter is set to true,
/// then the previous preconditioner will be used,
/// otherwise the new preconditioner will be constructed.
......@@ -269,7 +269,7 @@ namespace INMOST
///
/// Any changes to preconditioner parameters should happen before that point.
/// If you increase gmres_substep after this point, inner methods most likely will fail
void SetMatrix(Sparse::Matrix & A, bool OldPreconditioner = false);
void SetMatrix(Sparse::Matrix & A, bool ModifiedPattern = true, bool OldPreconditioner = false);
/// Solver the linear system: A*x = b.
/// Prior to this call you should call SetMatrix
///
......
This diff is collapsed.
......@@ -1424,16 +1424,14 @@ namespace INMOST
return *this;
}
void Solver::SetMatrix(Sparse::Matrix & A, bool OldPreconditioner)
void Solver::SetMatrix(Sparse::Matrix & A, bool ModifiedPattern, bool OldPreconditioner)
{
(void) OldPreconditioner;
bool ok = false;
#if defined(USE_SOLVER_PETSC)
if( _pack == PETSc )
{
bool modified_pattern = false;
for(Sparse::Matrix::iterator it = A.Begin(); it != A.End() && !modified_pattern; ++it)
modified_pattern |= it->modified_pattern;
bool modified_pattern = ModifiedPattern;
//~ if( A.comm != comm ) throw DifferentCommunicatorInSolver;
if( matrix_data == NULL )
{
......@@ -1531,10 +1529,7 @@ namespace INMOST
#if defined(USE_SOLVER_ANI)
if( _pack == ANI )
{
bool modified_pattern = false;
for(Sparse::Matrix::iterator it = A.Begin(); it != A.End() && !modified_pattern; ++it)
modified_pattern |= it->modified_pattern;
//~ if( A.comm != comm ) throw DifferentCommunicatorInSolver;
bool modified_pattern = ModifiedPattern;
if( matrix_data == NULL )
{
MatrixInitDataAni(&matrix_data,A.GetCommunicator(),A.GetName().c_str());
......@@ -1585,9 +1580,7 @@ namespace INMOST
#if defined(HAVE_SOLVER_FCBIILU2)
if( _pack == FCBIILU2 )
{
bool modified_pattern = false;
for(Sparse::Matrix::iterator it = A.Begin(); it != A.End() && !modified_pattern; ++it)
modified_pattern |= it->modified_pattern;
bool modified_pattern = ModifiedPattern;
//~ if( A.comm != comm ) throw DifferentCommunicatorInSolver;
if( matrix_data == NULL )
{
......@@ -1659,9 +1652,7 @@ namespace INMOST
#if defined(HAVE_SOLVER_K3BIILU2)
if( _pack == K3BIILU2 )
{
bool modified_pattern = false;
for(Sparse::Matrix::iterator it = A.Begin(); it != A.End() && !modified_pattern; ++it)
modified_pattern |= it->modified_pattern;
bool modified_pattern = ModifiedPattern;
//~ if( A.comm != comm ) throw DifferentCommunicatorInSolver;
if( matrix_data == NULL )
{
......@@ -1743,9 +1734,7 @@ namespace INMOST
#endif
//std::cout << Comm.MyPID() << " " << __FUNCTION__ << ":" << __LINE__ << std::endl;
//std::cout << Comm.MyPID() << " " << "Check modified pattern" << std::endl;
bool modified_pattern = false;
for(Sparse::Matrix::iterator it = A.Begin(); it != A.End() && !modified_pattern; ++it)
modified_pattern |= it->modified_pattern;
bool modified_pattern = ModifiedPattern;
INMOST_DATA_ENUM_TYPE mbeg,mend;
A.GetInterval(mbeg,mend);
//std::cout << Comm.MyPID() << " " << "Get interval " << mbeg << ":" << mend << std::endl;
......@@ -1857,7 +1846,6 @@ namespace INMOST
ok = true;
}
for(Sparse::Matrix::iterator it = A.Begin(); it != A.End(); it++) it->modified_pattern = false;
if(!ok) throw NotImplemented;
}
......