Commit 9c33d2c1 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Residual class annotation in inmost_autodiff.h

parent b2fb0655
...@@ -401,6 +401,86 @@ namespace INMOST ...@@ -401,6 +401,86 @@ namespace INMOST
return ret; return ret;
} }
#endif //USE_MESH #endif //USE_MESH
#if defined(USE_SOLVER)
void Residual::GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const
{
start = residual.GetFirstIndex();
end = residual.GetLastIndex();
}
void Residual::SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end)
{
jacobian.SetInterval(beg,end);
residual.SetInterval(beg,end);
}
void Residual::ClearResidual()
{
for(Sparse::Vector::iterator it = residual.Begin(); it != residual.End(); ++it) (*it) = 0.0;
}
void Residual::ClearJacobian()
{
for(Sparse::Matrix::iterator it = jacobian.Begin(); it != jacobian.End(); ++it) it->Clear();
}
void Residual::Clear()
{
#if defined(USE_OMP)
#pragma omp for
#endif //USE_OMP
for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k)
{
residual[k] = 0.0;
jacobian[k].Clear();
}
}
INMOST_DATA_REAL_TYPE Residual::Norm()
{
INMOST_DATA_REAL_TYPE ret = 0;
#if defined(USE_OMP)
#pragma omp parallel for reduction(+:ret)
#endif //USE_OMP
for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k)
ret += residual[k]*residual[k];
#if defined(USE_MPI)
INMOST_DATA_REAL_TYPE tmp = ret;
MPI_Allreduce(&tmp, &ret, 1, INMOST_MPI_DATA_REAL_TYPE, MPI_SUM, jacobian.GetCommunicator());
#endif
return sqrt(ret);
}
Residual & Residual::operator =(Residual const & other)
{
jacobian = other.jacobian;
residual = other.residual;
return *this;
}
void Residual::Rescale()
{
#if defined(USE_OMP)
#pragma omp parallel for
#endif //USE_OMP
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;
}
}
}
Residual::Residual(std::string name, INMOST_DATA_ENUM_TYPE start, INMOST_DATA_ENUM_TYPE end, INMOST_MPI_Comm _comm)
: jacobian(name,start,end,_comm),residual(name,start,end,_comm)
{
}
Residual::Residual(const Residual & other)
: jacobian(other.jacobian), residual(other.residual)
{
}
#endif //USE_SOLVER
}; };
#endif //USE_AUTODIFF #endif //USE_AUTODIFF
...@@ -14,104 +14,78 @@ namespace INMOST ...@@ -14,104 +14,78 @@ namespace INMOST
class Automatizator; //forward declaration class Automatizator; //forward declaration
#if defined(USE_SOLVER) #if defined(USE_SOLVER)
/// The Residual class provides a representation for array of residuals of nonlinear equations.
/// By working with the residual class you automatically assemble right hand side and
/// the jacobian of a nonlinear system of equation.
/// Jacobian matrix has a sparse representation.
/// \todo
/// 1. Extend for hessian calculation.
class Residual class Residual
{ {
Sparse::Matrix jacobian; Sparse::Matrix jacobian; ///< Jacobian matrix.
Sparse::Vector residual; Sparse::Vector residual; ///< Right hand side vector.
Sparse::LockService locks; Sparse::LockService locks; ///< Array of locks for openmp shared access.
public: 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) /// Constructor
: jacobian(name,start,end,_comm),residual(name,start,end,_comm) {} /// @param name Name for the matrix and right hand side. Can be used to set options for linear solver.
Residual(const Residual & other) /// @param start First index of the equation in the local partition. Use Automatizator::GetFirstIndex.
: jacobian(other.jacobian), residual(other.residual) /// @param end Last index of the equation in the local partition. Use Automatizator::GetLastIndex.
{} /// @param _comm MPI Communicator.
Residual & operator =(Residual const & other) Residual(std::string name = "", INMOST_DATA_ENUM_TYPE start = 0, INMOST_DATA_ENUM_TYPE end = 0, INMOST_MPI_Comm _comm = INMOST_MPI_COMM_WORLD);
{ /// Copy constructor.
jacobian = other.jacobian; /// \warning May be expensive if matrices are large.
residual = other.residual; Residual(const Residual & other);
return *this; /// Assignment operator.
} /// \warning May be expensive if matrices are large.
Residual & operator =(Residual const & other);
/// Retrive jacobian matrix. Use in Sparse::Solver::Solve function.
Sparse::Matrix & GetJacobian() {return jacobian;} Sparse::Matrix & GetJacobian() {return jacobian;}
/// Retrive jacobian matrix without right of modificaiton.
const Sparse::Matrix & GetJacobian() const {return jacobian;} const Sparse::Matrix & GetJacobian() const {return jacobian;}
/// Retrive right hand side vector. Use in Sparse::Solver::Solve function.
Sparse::Vector & GetResidual() {return residual;} Sparse::Vector & GetResidual() {return residual;}
/// Retrive right hand side vector without right of modification.
const Sparse::Vector & GetResidual() const {return residual;} const Sparse::Vector & GetResidual() const {return residual;}
/// Retrive the first index of the equations in the local partition.
INMOST_DATA_ENUM_TYPE GetFirstIndex() const {return residual.GetFirstIndex();} INMOST_DATA_ENUM_TYPE GetFirstIndex() const {return residual.GetFirstIndex();}
/// Retrive the last index of the equations in the local partition.
INMOST_DATA_ENUM_TYPE GetLastIndex() const {return residual.GetLastIndex();} INMOST_DATA_ENUM_TYPE GetLastIndex() const {return residual.GetLastIndex();}
void GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const /// Retrive the first and the last indices of the equations in the local partition
{ /// @param start The first index of the equations will be recorded here.
start = residual.GetFirstIndex(); /// @param end The last index of the equations will be recorded here.
end = residual.GetLastIndex(); void GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const;
} /// Assign the new first and last indices of the equations in the local partition.
void SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end) /// @param start The new first index of the equations.
{ /// @param end The new last index of the equations.
jacobian.SetInterval(beg,end); void SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end);
residual.SetInterval(beg,end); /// Retrive a residual value and a jacobian row corresponding to certain equation.
} /// @param row Equation number.
multivar_expression_reference operator [](INMOST_DATA_ENUM_TYPE row) /// @return A structure that can be used in or assigned an automatic differentiation expression.
{ __INLINE multivar_expression_reference operator [](INMOST_DATA_ENUM_TYPE row) {return multivar_expression_reference(residual[row],&jacobian[row]);}
return multivar_expression_reference(residual[row],&jacobian[row]); /// Zero out right hand side vector.
} void ClearResidual();
void ClearResidual() /// Remove all entries in jacobian matrix.
{ void ClearJacobian();
for(Sparse::Vector::iterator it = residual.Begin(); it != residual.End(); ++it) (*it) = 0.0; /// Zero out right hand side vector and remove all entries in jacobian matrix.
} void Clear();
void ClearJacobian() /// Compute the second norm of the right hand side vector over all of the processors.
{ INMOST_DATA_REAL_TYPE Norm();
for(Sparse::Matrix::iterator it = jacobian.Begin(); it != jacobian.End(); ++it) /// Normalize jacobian rows to unit second norms and scale right hand side accordingly.
it->Clear(); void Rescale();
} /// Initialize openmp locks.
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 parallel for reduction(+:ret)
#endif
for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k)
ret += residual[k]*residual[k];
#if defined(USE_MPI)
INMOST_DATA_REAL_TYPE tmp = ret;
MPI_Allreduce(&tmp, &ret, 1, INMOST_MPI_DATA_REAL_TYPE, MPI_SUM, jacobian.GetCommunicator());
#endif
return sqrt(ret);
}
/// Normalize entries in jacobian and right hand side
void Rescale()
{
#if defined(USE_OMP)
#pragma omp parallel 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;
}
}
}
void InitLocks() {locks.SetInterval(GetFirstIndex(),GetLastIndex());} void InitLocks() {locks.SetInterval(GetFirstIndex(),GetLastIndex());}
void Lock(INMOST_DATA_ENUM_TYPE pos) {if(!locks.Empty()) locks.Lock(pos);} /// Lock an equation to avoid simultaneous shared access.
void UnLock(INMOST_DATA_ENUM_TYPE pos) {if(!locks.Empty()) locks.UnLock(pos);} /// @param pos Equation number.
void TestLock(INMOST_DATA_ENUM_TYPE pos) {if(!locks.Empty()) locks.TestLock(pos);} __INLINE void Lock(INMOST_DATA_ENUM_TYPE pos) {if(!locks.Empty()) locks.Lock(pos);}
/// UnLock an equation to allow simultaneous shared access.
/// @param pos Equation number.
__INLINE void UnLock(INMOST_DATA_ENUM_TYPE pos) {if(!locks.Empty()) locks.UnLock(pos);}
/// Try to lock the equation.
/// @param pos Equation number.
/// @return True if equation was locked.
__INLINE bool TestLock(INMOST_DATA_ENUM_TYPE pos) {if(!locks.Empty()) return locks.TestLock(pos); return false;}
}; };
#endif #endif //USE_SOLVER
#if defined(USE_MESH) #if defined(USE_MESH)
/// The Automatizator class helps in defining primary unknowns of the model and /// The Automatizator class helps in defining primary unknowns of the model and
...@@ -123,26 +97,26 @@ namespace INMOST ...@@ -123,26 +97,26 @@ namespace INMOST
/// In addition this class provides automatic differentiation with acceleration structures /// In addition this class provides automatic differentiation with acceleration structures
/// that help evaluate expressions faster. /// that help evaluate expressions faster.
/// \todo /// \todo
/// 1. UnRegisterTag. /// 1. (test) UnRegisterTag.
/// 2. (test) Copy constructor. /// 2. (test) Copy constructor.
class Automatizator class Automatizator
{ {
private: private:
static Automatizator * CurrentAutomatizator; //< Currently used automatizator for automatic differentiation acceleration structures. static Automatizator * CurrentAutomatizator; ///< Currently used automatizator for automatic differentiation acceleration structures.
#if defined(USE_OMP) #if defined(USE_OMP)
std::vector<Sparse::RowMerger> merger; //< Automatic differentiation acceleration structures. std::vector<Sparse::RowMerger> merger; ///< Automatic differentiation acceleration structures.
#else #else
Sparse::RowMerger merger; //< Automatic differentiation acceleration structures. Sparse::RowMerger merger; ///< Automatic differentiation acceleration structures.
#endif #endif //USE_OMP
typedef struct{ Tag t, indices; MarkerType domain_mask; bool active;} tagdata; //< Pair of tag and it's mask marker. typedef struct{ Tag t, indices; MarkerType domain_mask; bool active;} tagdata; ///< Pair of tag and it's mask marker.
typedef std::vector<tagdata> tag_enum; //< A type for an array of registered tags. typedef std::vector<tagdata> tag_enum; ///< A type for an array of registered tags.
typedef std::vector<INMOST_DATA_ENUM_TYPE> del_enum; // A type for an array of deleted positions. typedef std::vector<INMOST_DATA_ENUM_TYPE> del_enum; ///< A type for an array of deleted positions.
private: private:
std::string name; //< Name of the automatizator. std::string name; ///< Name of the automatizator.
del_enum del_tags; //< Array of deleted positions. del_enum del_tags; ///< Array of deleted positions.
tag_enum reg_tags; //< Array of registered tags. tag_enum reg_tags; ///< Array of registered tags.
INMOST_DATA_ENUM_TYPE first_num; //< First index in unknowns of locally owned elements. INMOST_DATA_ENUM_TYPE first_num; ///< First index in unknowns of locally owned elements.
INMOST_DATA_ENUM_TYPE last_num; //< Last index in unknowns of locally owned elements. INMOST_DATA_ENUM_TYPE last_num; ///< Last index in unknowns of locally owned elements.
public: public:
/// Make a copy. /// Make a copy.
/// \warning /// \warning
...@@ -218,13 +192,13 @@ namespace INMOST ...@@ -218,13 +192,13 @@ namespace INMOST
/// Retrive acceleration structure for automatic differentiation. /// Retrive acceleration structure for automatic differentiation.
/// This structure can be used for operations on the matrix. /// This structure can be used for operations on the matrix.
/// @return Acceleration structure. /// @return Acceleration structure.
Sparse::RowMerger & GetMerger() __INLINE Sparse::RowMerger & GetMerger()
{ {
#if defined(USE_OMP) #if defined(USE_OMP)
return merger[omp_get_thread_num()]; return merger[omp_get_thread_num()];
#else #else
return merger; return merger;
#endif #endif //USE_OMP
} }
/// Remove global current automatizator used to set acceleration structures for automatic differentation. /// Remove global current automatizator used to set acceleration structures for automatic differentation.
static void RemoveCurrent() {CurrentAutomatizator = NULL;} static void RemoveCurrent() {CurrentAutomatizator = NULL;}
...@@ -239,7 +213,7 @@ namespace INMOST ...@@ -239,7 +213,7 @@ namespace INMOST
/// @return An array with indices corresponding to all registered tags. /// @return An array with indices corresponding to all registered tags.
std::vector<INMOST_DATA_ENUM_TYPE> ListRegisteredTags() const; std::vector<INMOST_DATA_ENUM_TYPE> ListRegisteredTags() const;
}; };
#endif #endif //USE_MESH
} //namespace INMOST } //namespace INMOST
#endif //USE_AUTODIFF #endif //USE_AUTODIFF
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment