Commit ecaafa49 authored by David Bommes's avatar David Bommes
Browse files

added TinyAD support for finite elements (TinyAD dependency required)

parent ecd800a4
Pipeline #17901 failed with stages
in 6 minutes and 13 seconds
......@@ -344,6 +344,15 @@ else ()
set (COMISO_CPLEX_CONFIG_FILE_SETTINGS "#define COMISO_CPLEX_AVAILABLE 0" )
endif ()
find_package (TINYAD)
if (TINYAD_FOUND)
set (COMISO_TINYAD_CONFIG_FILE_SETTINGS "#define COMISO_TINYAD_AVAILABLE 1" )
list( APPEND COMISO_INCLUDE_DIRECTORIES ${TINYAD_INCLUDE_DIRS} )
else()
set (COMISO_TINYAD_CONFIG_FILE_SETTINGS "#define COMISO_TINYAD_AVAILABLE 0" )
endif ()
set(TMP_CMAKE_FIND_LIBRARY_PREFIXES "${CMAKE_FIND_LIBRARY_PREFIXES}")
set(CMAKE_FIND_LIBRARY_PREFIXES lib "") #Our lapack library is called liblapack.lib. Is there a better way to find it than this?
if (NEED_LAPACK AND NOT SUITESPARSE_FOUND)
......
......@@ -27,3 +27,4 @@
@COMISO_COINUTILS_CONFIG_FILE_SETTINGS@
@COMISO_MOSEK_CONFIG_FILE_SETTINGS@
@COMISO_OSQP_CONFIG_FILE_SETTINGS@
@COMISO_TINYAD_CONFIG_FILE_SETTINGS@
//=============================================================================
//
// CLASS FINITEELEMENTHESSIANPROJECTION
//
//=============================================================================
#ifndef COMISO_FINITEELEMENTHESSIANPROJECTION_HH
#define COMISO_FINITEELEMENTHESSIANPROJECTION_HH
//== INCLUDES =================================================================
#include <CoMISo/Config/CoMISoDefines.hh>
#include<Eigen/Dense>
#include<Eigen/Sparse>
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================
namespace COMISO {
//== CLASS DEFINITION =========================================================
/** \class FINITEELEMENTHESSIANPROJECTION
Implements forward mode algorithmic differentiation (first and second order) for FiniteElements via TinyAD
A more elaborate description follows.
*/
template<class ElementT>
class FiniteElementHessianProjection : public ElementT
{
public:
// import dimensions
const static int NV = ElementT::NV;
const static int NC = ElementT::NC;
using VecI = Eigen::Matrix<size_t,NV,1>;
using VecV = Eigen::Matrix<double,NV,1>;
using VecC = Eigen::Matrix<double,NC,1>;
using Triplet = Eigen::Triplet<double>;
inline double eval_f(const VecV& _x, const VecC& _c) const
{
return ElementT::eval_f(_x,_c);
}
inline void eval_gradient(const VecV& _x, const VecC& _c, VecV& _g) const
{
return ElementT::eval_gradient(_x,_c,_g);
}
inline void eval_hessian (const VecV& _x, const VecC& _c, std::vector<Triplet>& _triplets) const
{
ElementT::eval_hessian(_x,_c,_triplets);
// convert to dense matrix
Eigen::Matrix<double,NV,NV> H;
H.setZero();
for(auto t : _triplets)
H(t.row(), t.col()) += t.value();
// make s.p.d. if necessary
if(!diagonally_dominant(H))
{
project_to_positive_semidefinite(H);
_triplets.clear();
for (unsigned int i = 0; i < NV; ++i)
for (unsigned int j = 0; j < NV; ++j)
_triplets.push_back(Triplet(i, j, H(i, j)));
}
}
inline double max_feasible_step(const VecV& _x, const VecV& _v, const VecC& _c)
{
return ElementT::max_feasible_step(_x, _v, _c);
}
template <class MatT>
inline bool diagonally_dominant(const MatT& _A) const
{
// only works for square matrices
assert(_A.cols() == _A.rows());
for(unsigned int i=0; i<_A.rows(); ++i)
{
double asum(0.0);
for(unsigned int j=0; j<_A.cols(); ++j)
asum += std::abs(_A(i,j));
double d = std::abs(_A(i,i));
if( 2*d < asum)
return false;
}
return true;
}
template <class MatT>
inline void project_to_positive_semidefinite( MatT& _A ) const
{
typename Eigen::SelfAdjointEigenSolver<MatT> es(_A);
typename Eigen::SelfAdjointEigenSolver<MatT>::RealVectorType ev = es.eigenvalues();
for (unsigned int i = 0; i < ev.size(); ++i)
{
ev[i] = std::max(ev[i], 0.0);
}
_A = es.eigenvectors() * ev.asDiagonal() * es.eigenvectors().transpose();
}
};
//=============================================================================
} // namespace COMISO
//=============================================================================
#endif // COMISO_FINITEELEMENTHESSIANPROJECTION_HH defined
//=============================================================================
//=============================================================================
//
// CLASS FINITEELEMENTTINYAD
//
//=============================================================================
#ifndef COMISO_FINITEELEMENTTINYAD_HH
#define COMISO_FINITEELEMENTTINYAD_HH
//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
#if COMISO_TINYAD_AVAILABLE
//== INCLUDES =================================================================
#include <CoMISo/Config/CoMISoDefines.hh>
#include<Eigen/Dense>
#include<Eigen/Sparse>
#include <TinyAD/TinyAD.hh>
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================
namespace COMISO {
//== CLASS DEFINITION =========================================================
/** \class FiniteElementTinyAD
Implements forward mode algorithmic differentiation (first and second order) for FiniteElements via TinyAD
A more elaborate description follows.
*/
template<class ElementT>
class FiniteElementTinyAD : public ElementT
{
public:
// import dimensions
const static int NV = ElementT::NV;
const static int NC = ElementT::NC;
using VecI = Eigen::Matrix<size_t,NV,1>;
using VecV = Eigen::Matrix<double,NV,1>;
using VecC = Eigen::Matrix<double,NC,1>;
using Triplet = Eigen::Triplet<double>;
using TADG = TinyAD::Scalar<NV,double,false>;
using TADH = TinyAD::Scalar<NV,double,true>;
inline double eval_f(const VecV& _x, const VecC& _c) const
{
return ElementT::eval_f(_x,_c);
}
inline void eval_gradient(const VecV& _x, const VecC& _c, VecV& _g) const
{
auto x_ad = TADG::make_active(_x);
auto f = ElementT::eval_f(x_ad,_c);
_g = f.grad;
}
inline void eval_hessian (const VecV& _x, const VecC& _c, std::vector<Triplet>& _triplets) const
{
auto x_ad = TADH::make_active(_x);
auto f = ElementT::eval_f(x_ad,_c);
auto H = f.Hess;
_triplets.clear();
for(unsigned int i=0; i<H.rows(); ++i)
for(unsigned int j=0; j<H.cols(); ++j)
_triplets.push_back(Triplet(i,j,H(i,j)));
}
inline double max_feasible_step(const VecV& _x, const VecV& _v, const VecC& _c)
{
return ElementT::max_feasible_step(_x, _v, _c);
}
};
//=============================================================================
} // namespace COMISO
//=============================================================================
#endif // COMISO_TINYAD_AVAILABLE
//=============================================================================
#endif // COMISO_FINITEELEMENTTINYAD_HH defined
//=============================================================================
......@@ -64,6 +64,12 @@ void NPTiming::store_result ( const double* _x )
print_statistics();
}
double NPTiming::max_feasible_step ( const double* _x, const double* _v)
{
return base_->max_feasible_step(_x, _v);
}
bool NPTiming::constant_gradient() const
{
return base_->constant_gradient();
......
......@@ -12,7 +12,6 @@
//== INCLUDES =================================================================
#include <Base/Utils/StopWatch.hh>
#include <CoMISo/Utils/gmm.hh>
#include "NProblemInterface.hh"
#include <CoMISo/Config/CoMISoDefines.hh>
......@@ -26,7 +25,7 @@ namespace COMISO {
/** \class NProblemGmmInterface NProblemGmmInterface.hh <ACG/.../NProblemGmmInterface.hh>
/** \class NPTiming NPTiming.hh <CoMISo/NSolver/NPTiming.hh>
Brief Description.
......@@ -54,6 +53,7 @@ public:
// advanced properties
virtual bool constant_gradient() const;
virtual bool constant_hessian() const;
virtual double max_feasible_step ( const double* _x, const double* _v);
void start_timing();
......
# - Try to find TINYAD
# Once done this will define
# TINYAD_FOUND - System has TINYAD
# TINYAD_INCLUDE_DIRS - The TINYAD include directories
if (TINYAD_INCLUDE_DIR)
# in cache already
set(TINYAD_FOUND TRUE)
set(TINYAD_INCLUDE_DIRS "${TINYAD_INCLUDE_DIR}" )
else (TINYAD_INCLUDE_DIR)
# Check if the base path is set
if ( NOT CMAKE_WINDOWS_LIBS_DIR )
# This is the base directory for windows library search used in the finders we shipp.
set(CMAKE_WINDOWS_LIBS_DIR "c:/libs" CACHE STRING "Default Library search dir on windows." )
endif()
find_path( TINYAD_INCLUDE_DIR
NAMES TinyAD/TinyAD.hh
PATHS $ENV{TINYAD_DIR}
$ENV{TINYAD_DIR}/include
/usr/include
/usr/include
/usr/local/include
)
set(TINYAD_INCLUDE_DIRS "${TINYAD_INCLUDE_DIR}" )
include(FindPackageHandleStandardArgs)
# handle the QUIETLY and REQUIRED arguments and set TINYAD_FOUND to TRUE
# if all listed variables are TRUE
find_package_handle_standard_args(TINYAD DEFAULT_MSG
TINYAD_INCLUDE_DIR)
mark_as_advanced(TINYAD_INCLUDE_DIR)
endif(TINYAD_INCLUDE_DIR)
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