Commit d89bb9ec by Kirill Terekhov

Some updates for hessian calculation

parent 0a8f0b00
......@@ -421,6 +421,10 @@ namespace INMOST
{
for(Sparse::Matrix::iterator it = jacobian.Begin(); it != jacobian.End(); ++it) it->Clear();
}
void Residual::ClearHessian()
{
for(Sparse::HessianMatrix::iterator it = hessian.Begin(); it != hessian.End(); ++it) it->Clear();
}
void Residual::Clear()
{
#if defined(USE_OMP)
......@@ -429,7 +433,8 @@ namespace INMOST
for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k)
{
residual[k] = 0.0;
jacobian[k].Clear();
if( !jacobian.Empty() ) jacobian[k].Clear();
if( !hessian.Empty() ) hessian[k].Clear();
}
}
INMOST_DATA_REAL_TYPE Residual::Norm()
......@@ -448,11 +453,12 @@ namespace INMOST
}
Residual & Residual::operator =(Residual const & other)
{
hessian = other.hessian;
jacobian = other.jacobian;
residual = other.residual;
return *this;
}
void Residual::Rescale()
void Residual::Rescale(INMOST_DATA_ENUM_TYPE p)
{
#if defined(USE_OMP)
#pragma omp parallel for
......@@ -460,24 +466,36 @@ namespace INMOST
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( p == ENUMUNDEF ) //infinite norm
{
for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q)
if( norm < fabs(jacobian[k].GetValue(q)) )
norm = fabs(jacobian[k].GetValue(q));
}
else //p-norm
{
for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q)
norm += pow(fabs(jacobian[k].GetValue(q)),p);
norm = pow(norm,1.0/p);
}
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;
if( !hessian.Empty() )
for(INMOST_DATA_ENUM_TYPE q = 0; q < hessian[k].Size(); ++q)
hessian[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)
: jacobian(name,start,end,_comm),residual(name,start,end,_comm), hessian(name,0,0,_comm)
{
}
Residual::Residual(const Residual & other)
: jacobian(other.jacobian), residual(other.residual)
: jacobian(other.jacobian), residual(other.residual), hessian(other.hessian)
{
}
#endif //USE_SOLVER
......
......@@ -22,6 +22,7 @@ namespace INMOST
/// 1. Extend for hessian calculation.
class Residual
{
Sparse::HessianMatrix hessian; ///< Hessian matrix
Sparse::Matrix jacobian; ///< Jacobian matrix.
Sparse::Vector residual; ///< Right hand side vector.
Sparse::LockService locks; ///< Array of locks for openmp shared access.
......@@ -38,14 +39,6 @@ namespace INMOST
/// 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;}
/// Retrive jacobian matrix without right of modificaiton.
const Sparse::Matrix & GetJacobian() const {return jacobian;}
/// Retrive right hand side vector. Use in Sparse::Solver::Solve function.
Sparse::Vector & GetResidual() {return residual;}
/// Retrive right hand side vector without right of modification.
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();}
/// Retrive the last index of the equations in the local partition.
......@@ -61,17 +54,33 @@ namespace INMOST
/// Retrive a residual value and a jacobian row corresponding to certain equation.
/// @param row Equation number.
/// @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]);}
__INLINE multivar_expression_reference operator [](INMOST_DATA_ENUM_TYPE row)
{return multivar_expression_reference(residual[row],&jacobian[row]);}
/// Retrive hessian matrix. Use in nonlinear solver.
Sparse::HessianMatrix & GetHessian() {return hessian;}
/// Retrive hessian matrix without right of modificaiton.
const Sparse::HessianMatrix & GetHessian() const {return hessian;}
/// Retrive jacobian matrix. Use in Sparse::Solver::Solve function.
Sparse::Matrix & GetJacobian() {return jacobian;}
/// Retrive jacobian matrix without right of modificaiton.
const Sparse::Matrix & GetJacobian() const {return jacobian;}
/// Retrive right hand side vector. Use in Sparse::Solver::Solve function.
Sparse::Vector & GetResidual() {return residual;}
/// Retrive right hand side vector without right of modification.
const Sparse::Vector & GetResidual() const {return residual;}
/// Zero out right hand side vector.
void ClearResidual();
/// Remove all entries in jacobian matrix.
void ClearJacobian();
/// Remove all entries in hessian matrix.
void ClearHessian();
/// Zero out right hand side vector and remove all entries in jacobian matrix.
void Clear();
/// Compute the second norm of the right hand side vector over all of the processors.
INMOST_DATA_REAL_TYPE Norm();
/// Normalize jacobian rows to unit second norms and scale right hand side accordingly.
void Rescale();
/// Normalize jacobian rows to unit p-norms and scale right hand side accordingly.
/// Use ENUMUNDEF as p to scale to infinite-norm.
void Rescale(INMOST_DATA_ENUM_TYPE p = 2);
/// Initialize openmp locks.
void InitLocks() {locks.SetInterval(GetFirstIndex(),GetLastIndex());}
/// Lock an equation to avoid simultaneous shared access.
......
......@@ -15,17 +15,43 @@ namespace INMOST
template<> struct Promote<INMOST_DATA_REAL_TYPE, INMOST_DATA_INTEGER_TYPE> {typedef INMOST_DATA_REAL_TYPE type;};
template<> struct Promote<INMOST_DATA_REAL_TYPE, INMOST_DATA_REAL_TYPE> {typedef INMOST_DATA_REAL_TYPE type;};
#if defined(USE_AUTODIFF)
//For INMOST_DATA_INTEGER_TYPE
template<> struct Promote<INMOST_DATA_INTEGER_TYPE, multivar_expression_reference> {typedef variable type;};
template<> struct Promote<INMOST_DATA_INTEGER_TYPE, variable> {typedef variable type;};
template<> struct Promote<INMOST_DATA_INTEGER_TYPE, hessian_multivar_expression_reference> {typedef hessian_variable type;};
template<> struct Promote<INMOST_DATA_INTEGER_TYPE, hessian_variable> {typedef hessian_variable type;};
//For INMOST_DATA_REAL_TYPE
template<> struct Promote<INMOST_DATA_REAL_TYPE, variable> {typedef variable type;};
template<> struct Promote<INMOST_DATA_REAL_TYPE, multivar_expression_reference> {typedef variable type;};
template<> struct Promote<INMOST_DATA_REAL_TYPE, hessian_multivar_expression_reference> {typedef hessian_variable type;};
template<> struct Promote<INMOST_DATA_REAL_TYPE, hessian_variable> {typedef hessian_variable type;};
//For multivar_expression_reference
template<> struct Promote<multivar_expression_reference, INMOST_DATA_INTEGER_TYPE> {typedef variable type;};
template<> struct Promote<multivar_expression_reference, INMOST_DATA_REAL_TYPE> {typedef variable type;};
template<> struct Promote<multivar_expression_reference, multivar_expression_reference> {typedef variable type;};
template<> struct Promote<multivar_expression_reference, variable> {typedef variable type;};
template<> struct Promote<multivar_expression_reference, hessian_multivar_expression_reference> {typedef variable type;};
template<> struct Promote<multivar_expression_reference, hessian_variable> {typedef variable type;};
//For variable
template<> struct Promote<variable, INMOST_DATA_INTEGER_TYPE> {typedef variable type;};
template<> struct Promote<variable, INMOST_DATA_REAL_TYPE> {typedef variable type;};
template<> struct Promote<variable, multivar_expression_reference> {typedef variable type;};
template<> struct Promote<variable, variable> {typedef variable type;};
template<> struct Promote<INMOST_DATA_INTEGER_TYPE, hessian_variable> {typedef hessian_variable type;};
template<> struct Promote<INMOST_DATA_REAL_TYPE, hessian_variable> {typedef hessian_variable type;};
template<> struct Promote<variable, hessian_multivar_expression_reference> {typedef hessian_variable type;};
template<> struct Promote<variable, hessian_variable> {typedef hessian_variable type;};
//For hessian_multivar_expression_reference
template<> struct Promote<hessian_multivar_expression_reference, INMOST_DATA_INTEGER_TYPE> {typedef hessian_variable type;};
template<> struct Promote<hessian_multivar_expression_reference, INMOST_DATA_REAL_TYPE> {typedef hessian_variable type;};
template<> struct Promote<hessian_multivar_expression_reference, multivar_expression_reference> {typedef hessian_variable type;};
template<> struct Promote<hessian_multivar_expression_reference, variable> {typedef hessian_variable type;};
template<> struct Promote<hessian_multivar_expression_reference, hessian_multivar_expression_reference> {typedef hessian_variable type;};
template<> struct Promote<hessian_multivar_expression_reference, hessian_variable> {typedef hessian_variable type;};
//For hessian_variable
template<> struct Promote<hessian_variable, INMOST_DATA_INTEGER_TYPE> {typedef hessian_variable type;};
template<> struct Promote<hessian_variable, INMOST_DATA_REAL_TYPE> {typedef hessian_variable type;};
template<> struct Promote<hessian_variable, multivar_expression_reference> {typedef hessian_variable type;};
template<> struct Promote<hessian_variable, variable> {typedef hessian_variable type;};
template<> struct Promote<hessian_variable, hessian_multivar_expression_reference> {typedef hessian_variable type;};
template<> struct Promote<hessian_variable, hessian_variable> {typedef hessian_variable type;};
#endif
......
This source diff could not be displayed because it is too large. You can view the blob instead.
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