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

New features and fixes.

Changed how matrix annotation is handled to decrease memory required to
store data of type DATA_VARIABLE.

Added multivar_expression_reference that can reference individual matrix
row.

Added Residual class.

Fixed that Sparse::Matrix::Load and Sparse::Vector::Load may required
MPI initialization.
parent f53db33d
This diff is collapsed.
......@@ -19,7 +19,7 @@ namespace INMOST
void SetAutodiffPrint(bool set) {print_ad_ctor = set;}
bool CheckCurrentAutomatizator() {return Automatizator::HaveCurrent();}
void multivar_expression::FromBasicExpression(const basic_expression & expr)
void FromBasicExpression(Sparse::Row & entries, const basic_expression & expr)
{
Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
expr.GetDerivative(1.0,merger);
......@@ -27,7 +27,7 @@ namespace INMOST
merger.Clear();
}
void multivar_expression::AddBasicExpression(INMOST_DATA_REAL_TYPE multme, INMOST_DATA_REAL_TYPE multit, const basic_expression & expr)
void AddBasicExpression(Sparse::Row & entries, INMOST_DATA_REAL_TYPE multme, INMOST_DATA_REAL_TYPE multit, const basic_expression & expr)
{
Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
merger.PushRow(multme,entries);
......@@ -36,10 +36,10 @@ namespace INMOST
merger.Clear();
}
void multivar_expression::FromGetDerivative(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
void FromGetDerivative(const basic_expression & expr, INMOST_DATA_REAL_TYPE mult, Sparse::Row & r)
{
Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
GetDerivative(mult,merger);
expr.GetDerivative(mult,merger);
merger.AddRow(1.0,r);
merger.RetriveRow(r);
merger.Clear();
......
......@@ -660,7 +660,7 @@ namespace INMOST
for(enumerator l = 0; l < m; ++l)
{
if( fabs(get_value((*this)(k,l))) > threshold )
std::cout << std::setw(10) << (*this)(k,l);
std::cout << std::setw(10) << get_value((*this)(k,l));
else
std::cout << std::setw(10) << 0;
std::cout << " ";
......
......@@ -24,9 +24,7 @@
#if defined(USE_AUTODIFF)
namespace INMOST
{
bool CheckCurrentAutomatizator();
//bool GetAutodiffPrint();
//void SetAutodiffPrint(bool set);
class basic_expression
{
......@@ -39,6 +37,13 @@ namespace INMOST
virtual ~basic_expression() {}//if( GetAutodiffPrint() ) std::cout << this << " Destroied" << std::endl;}
};
bool CheckCurrentAutomatizator();
void FromBasicExpression(Sparse::Row & entries, const basic_expression & expr);
void AddBasicExpression(Sparse::Row & entries, INMOST_DATA_REAL_TYPE multme, INMOST_DATA_REAL_TYPE multit, const basic_expression & expr);
void FromGetDerivative(const basic_expression & expr, INMOST_DATA_REAL_TYPE mult, Sparse::Row & r);
//bool GetAutodiffPrint();
//void SetAutodiffPrint(bool set);
template<class Derived>
class shell_expression : virtual public basic_expression
{
......@@ -95,13 +100,11 @@ namespace INMOST
}
};
class multivar_expression : public shell_expression<multivar_expression>
{
INMOST_DATA_REAL_TYPE value;
Sparse::Row entries;
void FromBasicExpression(const basic_expression & expr);
void AddBasicExpression(INMOST_DATA_REAL_TYPE multme, INMOST_DATA_REAL_TYPE multit, const basic_expression & expr);
void FromGetDerivative(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const;
public:
multivar_expression() :value(0) {}
multivar_expression(INMOST_DATA_REAL_TYPE pvalue) : value(pvalue) {}
......@@ -117,7 +120,7 @@ namespace INMOST
{
value = expr.GetValue();
if( CheckCurrentAutomatizator() )
FromBasicExpression(expr); //Optimized version
FromBasicExpression(entries,expr); //Optimized version
else expr.GetDerivative(1.0,entries);
}
__INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
......@@ -130,7 +133,7 @@ namespace INMOST
__INLINE void GetDerivative(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
{
if( CheckCurrentAutomatizator() )
FromGetDerivative(mult,r);
FromGetDerivative(*this,mult,r);
else
{
for(Sparse::Row::const_iterator it = entries.Begin(); it != entries.End(); ++it)
......@@ -147,7 +150,7 @@ namespace INMOST
{
value = expr.GetValue();
if( CheckCurrentAutomatizator() )
FromBasicExpression(expr);
FromBasicExpression(entries,expr);
else
{
Sparse::Row tmp;
......@@ -168,7 +171,7 @@ namespace INMOST
{
value += expr.GetValue();
if( CheckCurrentAutomatizator() )
AddBasicExpression(1.0,1.0,expr);
AddBasicExpression(entries,1.0,1.0,expr);
else
{
Sparse::Row tmp(entries);
......@@ -181,7 +184,7 @@ namespace INMOST
{
value -= expr.GetValue();
if( CheckCurrentAutomatizator() )
AddBasicExpression(1.0,-1.0,expr);
AddBasicExpression(entries,1.0,-1.0,expr);
else
{
Sparse::Row tmp(entries);
......@@ -194,7 +197,7 @@ namespace INMOST
{
INMOST_DATA_REAL_TYPE lval = value, rval = expr.GetValue();
if( CheckCurrentAutomatizator() )
AddBasicExpression(rval,lval,expr);
AddBasicExpression(entries,rval,lval,expr);
else
{
Sparse::Row tmp(entries);
......@@ -211,7 +214,7 @@ namespace INMOST
INMOST_DATA_REAL_TYPE reciprocial_rval = 1.0/rval;
value *= reciprocial_rval;
if( CheckCurrentAutomatizator() )
AddBasicExpression(reciprocial_rval,-value*reciprocial_rval,expr);
AddBasicExpression(entries,reciprocial_rval,-value*reciprocial_rval,expr);
else
{
Sparse::Row tmp(entries);
......@@ -299,6 +302,161 @@ namespace INMOST
}
};
class multivar_expression_reference : public shell_expression<multivar_expression_reference>
{
INMOST_DATA_REAL_TYPE & value;
Sparse::Row & entries;
public:
/// Constructor, set links to the provided value and 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)
: value(other.value), entries(other.entries) {}
/// Retrive value
__INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
/// Set value without changing derivatives
__INLINE void SetValue(INMOST_DATA_REAL_TYPE val) { value = val; }
/// Retrive derivatives with multiplier into Sparse::RowMerger structure.
__INLINE void GetDerivative(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
{
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.
__INLINE void GetDerivative(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
{
if( CheckCurrentAutomatizator() )
FromGetDerivative(*this,mult,r);
else
{
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();
return *this;
}
__INLINE multivar_expression_reference & operator = (basic_expression const & expr)
{
value = expr.GetValue();
if( CheckCurrentAutomatizator() )
FromBasicExpression(entries,expr);
else
{
Sparse::Row tmp;
expr.GetDerivative(1.0,tmp);
entries.Swap(tmp);
}
return *this;
}
__INLINE multivar_expression_reference & operator = (multivar_expression_reference const & other)
{
value = other.GetValue();
entries = other.GetRow();
return *this;
}
__INLINE multivar_expression_reference & operator = (multivar_expression const & other)
{
value = other.GetValue();
entries = other.GetRow();
return *this;
}
__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);
else
{
Sparse::Row tmp(entries);
expr.GetDerivative(1.0,tmp);
entries.Swap(tmp);
}
return *this;
}
__INLINE multivar_expression_reference & operator -=(basic_expression const & expr)
{
value -= expr.GetValue();
if( CheckCurrentAutomatizator() )
AddBasicExpression(entries,1.0,-1.0,expr);
else
{
Sparse::Row tmp(entries);
expr.GetDerivative(-1.0,tmp);
entries.Swap(tmp);
}
return *this;
}
__INLINE multivar_expression_reference & operator *=(basic_expression const & expr)
{
INMOST_DATA_REAL_TYPE lval = value, rval = expr.GetValue();
if( CheckCurrentAutomatizator() )
AddBasicExpression(entries,rval,lval,expr);
else
{
Sparse::Row tmp(entries);
for(Sparse::Row::iterator it = tmp.Begin(); it != tmp.End(); ++it) it->second *= rval;
expr.GetDerivative(lval,tmp);
entries.Swap(tmp);
}
value *= rval;
return *this;
}
__INLINE multivar_expression_reference & operator /=(basic_expression const & expr)
{
INMOST_DATA_REAL_TYPE rval = expr.GetValue();
INMOST_DATA_REAL_TYPE reciprocial_rval = 1.0/rval;
value *= reciprocial_rval;
if( CheckCurrentAutomatizator() )
AddBasicExpression(entries,reciprocial_rval,-value*reciprocial_rval,expr);
else
{
Sparse::Row tmp(entries);
for(Sparse::Row::iterator it = tmp.Begin(); it != tmp.End(); ++it) it->second *= reciprocial_rval;
expr.GetDerivative(-value*reciprocial_rval,tmp);
entries.Swap(tmp);
}
return *this;
}
__INLINE multivar_expression_reference & operator +=(INMOST_DATA_REAL_TYPE right)
{
value += right;
return *this;
}
__INLINE multivar_expression_reference & operator -=(INMOST_DATA_REAL_TYPE right)
{
value -= 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;
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;
return *this;
}
bool check_nans() const
{
if( value != value ) return true;
for(Sparse::Row::iterator it = entries.Begin(); it != entries.End(); ++it)
if( it->second != it->second ) return true;
return false;
}
};
template<class A>
class const_multiplication_expression : public shell_expression<const_multiplication_expression<A> >
{
......@@ -1103,6 +1261,65 @@ 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()
{
ClearResidual();
ClearJacobian();
}
INMOST_DATA_REAL_TYPE Norm()
{
INMOST_DATA_REAL_TYPE ret = 0;
for(Sparse::Vector::iterator it = residual.Begin(); it != residual.End(); ++it) ret += (*it)*(*it);
return sqrt(ret);
}
};
}
......
......@@ -7,8 +7,6 @@
#if defined(USE_SOLVER)
#define MTX_ALLOW_ANNOTATION
#define MTX_ANNOTATION_SIZE 1024
namespace INMOST
{
......@@ -59,9 +57,12 @@ namespace INMOST
void SetInterval(INMOST_DATA_ENUM_TYPE start, INMOST_DATA_ENUM_TYPE end) {assert(start<=end); data.set_interval_beg(start); data.set_interval_end(end);}
/// Get the start and the end of the distributed vector interval.
void GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const {start = data.get_interval_beg(); end = data.get_interval_end();}
/// Move starting position of local indexes
void ShiftInterval(INMOST_DATA_ENUM_TYPE shift) {data.shift_interval(shift);}
/// Get the first index of the distributed vector interval.
INMOST_DATA_ENUM_TYPE GetFirstIndex() const {return data.get_interval_beg();}
/// Get the last index of the distributed vector interval.
INMOST_DATA_ENUM_TYPE GetLastIndex() const {return data.get_interval_end();}
/// Get the communicator which the vector is associated with.
INMOST_MPI_Comm GetCommunicator() const {return comm;}
......@@ -90,10 +91,7 @@ namespace INMOST
class Row
{
public:
#if defined(MTX_ALLOW_ANNOTATION)
char annotation[MTX_ANNOTATION_SIZE];
#endif
/// Entry of the sparse matrix row.
/// Entry of the sparse matrix row.
typedef struct entry_s
{
INMOST_DATA_ENUM_TYPE first; ///< the column number of the row element.
......@@ -133,10 +131,7 @@ namespace INMOST
bool marker;
Entries data;
public:
std::string GetAnnotation();
void SetAnnotation(std::string input);
void Report() {data.report_addr();}
void Report() {data.report_addr();}
void SetMarker() { marker = true; }
void RemMarker() { marker = false; }
bool GetMarker() { return marker; }
......@@ -244,12 +239,14 @@ namespace INMOST
class Matrix
{
public:
typedef interval<INMOST_DATA_ENUM_TYPE,std::string> Text;
typedef interval<INMOST_DATA_ENUM_TYPE,Row> Rows;
typedef Rows::iterator iterator;
typedef Rows::const_iterator const_iterator;
private:
INMOST_MPI_Comm comm;
Rows data;
Text text;
std::string name;
bool is_parallel;
public:
......@@ -274,18 +271,21 @@ namespace INMOST
const_iterator Begin() const {return data.begin();}
const_iterator End() const {return data.end();}
/// Set the start and the end row numbers of the distributed matrix interval.
void SetInterval(INMOST_DATA_ENUM_TYPE start, INMOST_DATA_ENUM_TYPE end) {data.set_interval_beg(start); data.set_interval_end(end);}
void SetInterval(INMOST_DATA_ENUM_TYPE start, INMOST_DATA_ENUM_TYPE end);
/// Get the start and the end row numbers of the distributed matrix interval.
void GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const {start = data.get_interval_beg(); end = data.get_interval_end();}
void ShiftInterval(INMOST_DATA_ENUM_TYPE shift) {data.shift_interval(shift);}
void ShiftInterval(INMOST_DATA_ENUM_TYPE shift);
/// Get the first row index of the distributed matrix interval.
INMOST_DATA_ENUM_TYPE GetFirstIndex() const {return data.get_interval_beg();}
/// Get the last row index of the distributed matrix interval.
INMOST_DATA_ENUM_TYPE GetLastIndex() const {return data.get_interval_end();}
/// Get the communicator which the matrix is associated with.
INMOST_MPI_Comm GetCommunicator() const {return comm;}
void MoveRows(INMOST_DATA_ENUM_TYPE from, INMOST_DATA_ENUM_TYPE to, INMOST_DATA_ENUM_TYPE size); //for parallel
void Swap(Matrix & other)
{
data.swap(other.data);
text.swap(other.text);
name.swap(other.name);
INMOST_MPI_Comm ctmp = comm;
comm = other.comm;
......@@ -311,8 +311,12 @@ namespace INMOST
void Save(std::string file);
/// Check that matrix is in parallel state.
bool & isParallel() { return is_parallel; }
const bool & isParallel() const { return is_parallel; }
/// Get the matrix name specified in the main constructor.
std::string GetName() {return name;}
std::string GetName() const {return name;}
std::string & Annotation(INMOST_DATA_ENUM_TYPE row);
const std::string & Annotation(INMOST_DATA_ENUM_TYPE row) const;
};
/// This class may be used to sum multiple sparse rows.
......
......@@ -7,6 +7,7 @@ namespace INMOST
{
namespace Sparse
{
const std::string stubstring = "";
bool _hasRowEntryType = false;
INMOST_MPI_Type RowEntryType = INMOST_MPI_DATATYPE_NULL;
......@@ -347,6 +348,13 @@ namespace INMOST
INMOST_DATA_ENUM_TYPE i = to + size, j = from + size;
if( size > 0 && to != from )
while( j != from ) data[--j].MoveRow(data[--i]);
if( !text.empty() )
{
i = to + size;
j = from + size;
if( size > 0 && to != from )
while( j != from ) text[--j] = text[--i];
}
}
......@@ -360,7 +368,9 @@ namespace INMOST
INMOST_DATA_REAL_TYPE val;
int size = 1, rank = 0;
#if defined(USE_MPI)
if( mend == ENUMUNDEF && mbeg == ENUMUNDEF )
int flag = 0;
MPI_Initialized(&flag);
if( flag && mend == ENUMUNDEF && mbeg == ENUMUNDEF )
{
MPI_Comm_rank(GetCommunicator(),&rank);
MPI_Comm_size(GetCommunicator(),&size);
......@@ -410,7 +420,9 @@ namespace INMOST
INMOST_DATA_REAL_TYPE val;
int size = 1, rank = 0;
#if defined(USE_MPI)
if( mend == ENUMUNDEF && mbeg == ENUMUNDEF )
int flag = 0;
MPI_Initialized(&flag);
if( flag && mend == ENUMUNDEF && mbeg == ENUMUNDEF )
{
MPI_Comm_rank(GetCommunicator(),&rank);
MPI_Comm_size(GetCommunicator(),&size);
......@@ -545,9 +557,7 @@ namespace INMOST
mtx.precision(15);
for(iterator it = Begin(); it != End(); ++it)
{
#if defined(MTX_ALLOW_ANNOTATION)
if( it->GetAnnotation() != "" ) mtx << "% " << it->GetAnnotation() << "\n";
#endif
if( !text.empty() ) mtx << "% " << Annotation(it-Begin()).c_str() << "\n";
for(Row::iterator jt = it->Begin(); jt != it->End(); ++jt)
{
mtx << row << " " << jt->first+1 << " " << jt->second << "\n";
......@@ -621,25 +631,22 @@ namespace INMOST
output << mtx.rdbuf();
#endif
}
std::string Row::GetAnnotation()
std::string & Matrix::Annotation(INMOST_DATA_ENUM_TYPE row)
{
#if defined(MTX_ALLOW_ANNOTATION)
return std::string(annotation,strnlen(annotation,MTX_ANNOTATION_SIZE));
#else
return "";
#endif
if( text.empty() )
{
text.set_interval_beg(GetFirstIndex());
text.set_interval_end(GetLastIndex());
}
return text[row];
}
void Row::SetAnnotation(std::string input)
const std::string & Matrix::Annotation(INMOST_DATA_ENUM_TYPE row) const
{
#if defined(MTX_ALLOW_ANNOTATION)
strncpy(annotation,input.c_str(),MTX_ANNOTATION_SIZE);
#endif
if( text.empty() ) return stubstring;
return text[row];
}
Row::Row() :data()
{
#if defined(MTX_ALLOW_ANNOTATION)
annotation[0] = '\0';
#endif
#if defined(USE_OMP)
omp_init_lock(&lock);
#endif
......@@ -647,9 +654,6 @@ namespace INMOST
}
Row::Row(const Row & other) :marker(other.marker),data(other.data)
{
#if defined(MTX_ALLOW_ANNOTATION)
strncpy(annotation,other.annotation,MTX_ANNOTATION_SIZE);
#endif
#if defined(USE_OMP)
omp_init_lock(&lock);
#endif
......@@ -657,9 +661,6 @@ namespace INMOST
}
Row::Row(entry * pbegin, entry * pend) :data(pbegin, pend)
{
#if defined(MTX_ALLOW_ANNOTATION)
annotation[0] = '\0';
#endif
#if defined(USE_OMP)
omp_init_lock(&lock);
#endif
......@@ -669,9 +670,6 @@ namespace INMOST
{
data = other.data;
marker = other.marker;
#if defined(MTX_ALLOW_ANNOTATION)
strncpy(annotation,other.annotation,MTX_ANNOTATION_SIZE);
#endif
return *this;
}
Row::~Row()
......@@ -686,15 +684,26 @@ namespace INMOST
bool tmp = marker;
marker = other.marker;
other.marker = tmp;
#if defined(MTX_ALLOW_ANNOTATION)
char tmpstr[MTX_ANNOTATION_SIZE];
strncpy(tmpstr,annotation,MTX_ANNOTATION_SIZE);
strncpy(annotation,other.annotation,MTX_ANNOTATION_SIZE);
strncpy(other.annotation,tmpstr,MTX_ANNOTATION_SIZE);
#endif
#if defined(USE_OMP)
//swap locks?
#endif
}
void Matrix::SetInterval(INMOST_DATA_ENUM_TYPE start, INMOST_DATA_ENUM_TYPE end)
{
data.set_interval_beg(start);
data.set_interval_end(end);
if( !text.empty() )
{
text.set_interval_beg(start);
text.set_interval_end(end);
}
}