EigenMatrix::EigenMatrixContractedHessianImp< T1, T2, Dim > Struct Template Reference

#include "src/ftensor/src/MatrixFunctionTemplate.hpp"

Collaboration diagram for EigenMatrix::EigenMatrixContractedHessianImp< T1, T2, Dim >:

Public Types
using	Val = const FTensor::Tensor1< T1, Dim >
using	Vec = const FTensor::Tensor2< T2, Dim, Dim >
using	Fun = boost::function< double(const double)>
using	V = double

Public Member Functions
	EigenMatrixContractedHessianImp (Val &t_val, Vec &t_vec)
template<typename T >
auto	getDiffDiffMat (Fun f, Fun d_f, Fun dd_f, T &t_S)

Static Public Attributes
static constexpr int	sizeSymm = (Dim * (Dim + 1)) / 2

Private Attributes
Val &	tVal
Vec &	tVec

Detailed Description

template<typename T1, typename T2, int Dim>
struct EigenMatrix::EigenMatrixContractedHessianImp< T1, T2, Dim >

Compute S:D2f(A) directly, without materialising sixth-order derivatives.

Definition at line 797 of file MatrixFunctionTemplate.hpp.

Member Typedef Documentation

Fun

template<typename T1 , typename T2 , int Dim>

using EigenMatrix::EigenMatrixContractedHessianImp< T1, T2, Dim >::Fun = boost::function<double(const double)>

Definition at line 801 of file MatrixFunctionTemplate.hpp.

V

template<typename T1 , typename T2 , int Dim>

using EigenMatrix::EigenMatrixContractedHessianImp< T1, T2, Dim >::V = double

Definition at line 802 of file MatrixFunctionTemplate.hpp.

Val

template<typename T1 , typename T2 , int Dim>

using EigenMatrix::EigenMatrixContractedHessianImp< T1, T2, Dim >::Val = const FTensor::Tensor1<T1, Dim>

Definition at line 799 of file MatrixFunctionTemplate.hpp.

Vec

template<typename T1 , typename T2 , int Dim>

using EigenMatrix::EigenMatrixContractedHessianImp< T1, T2, Dim >::Vec = const FTensor::Tensor2<T2, Dim, Dim>

Definition at line 800 of file MatrixFunctionTemplate.hpp.

Constructor & Destructor Documentation

EigenMatrixContractedHessianImp()

template<typename T1 , typename T2 , int Dim>

EigenMatrix::EigenMatrixContractedHessianImp< T1, T2, Dim >::EigenMatrixContractedHessianImp ( Val &  t_val, Vec &  t_vec  )

inline

Definition at line 806 of file MatrixFunctionTemplate.hpp.

      : tVal(t_val), tVec(t_vec) {}

Member Function Documentation

getDiffDiffMat()

template<typename T1 , typename T2 , int Dim>

template<typename T >

auto EigenMatrix::EigenMatrixContractedHessianImp< T1, T2, Dim >::getDiffDiffMat ( Fun  f, Fun  d_f, Fun  dd_f, T &  t_S  )

inline

Definition at line 810 of file MatrixFunctionTemplate.hpp.

                                                               {
    // Coalesce only eigenvalues equal within the same tolerance as getUniqNb.
    V group_ref[Dim];
    V group_val[Dim]{};
    int group_size[Dim]{};
    int group_of[Dim];
    int nb_groups = 0;
 
    auto is_equal = [](const V a, const V b) {
      const V eps = 100 * std::numeric_limits<V>::epsilon();
      const V scale = std::max(V(1), std::max(std::abs(a), std::abs(b)));
      return std::abs(a - b) <= eps * scale;
    };
 
    for (int aa = 0; aa != Dim; ++aa) {
      int gg = 0;
      for (; gg != nb_groups; ++gg)
        if (is_equal(tVal(aa), group_ref[gg]))
          break;
 
      if (gg == nb_groups) {
        group_ref[gg] = tVal(aa);
        ++nb_groups;
      }
      group_of[aa] = gg;
      group_val[gg] += tVal(aa);
      ++group_size[gg];
    }
 
    V f_val[Dim];
    V df_val[Dim];
    V ddf_val[Dim];
    for (int aa = 0; aa != nb_groups; ++aa) {
      group_val[aa] /= group_size[aa];
      f_val[aa] = f(group_val[aa]);
      df_val[aa] = d_f(group_val[aa]);
      ddf_val[aa] = dd_f(group_val[aa]);
    }
 
    auto first_divided_difference = [&](const int aa, const int bb) {
      if (aa == bb)
        return df_val[aa];
      return (f_val[aa] - f_val[bb]) / (group_val[aa] - group_val[bb]);
    };
 
    auto second_divided_difference = [&](const int aa, const int bb,
                                         const int cc) {
      if (aa == bb && bb == cc)
        return ddf_val[aa] / 2;
 
      if (aa == bb) {
        const V d1 = first_divided_difference(aa, cc);
        return (df_val[aa] - d1) / (group_val[aa] - group_val[cc]);
      }
      if (aa == cc) {
        const V d1 = first_divided_difference(aa, bb);
        return (df_val[aa] - d1) / (group_val[aa] - group_val[bb]);
      }
      if (bb == cc) {
        const V d1 = first_divided_difference(aa, bb);
        return (d1 - df_val[bb]) / (group_val[aa] - group_val[bb]);
      }
 
      const V d1_ab = first_divided_difference(aa, bb);
      const V d1_bc = first_divided_difference(bb, cc);
      return (d1_ab - d1_bc) / (group_val[aa] - group_val[cc]);
    };
 
    V ddf[Dim][Dim][Dim];
    for (int aa = 0; aa != Dim; ++aa)
      for (int bb = 0; bb != Dim; ++bb)
        for (int cc = 0; cc != Dim; ++cc)
          ddf[aa][bb][cc] = second_divided_difference(
              group_of[aa], group_of[bb], group_of[cc]);
 
    V s_hat[Dim][Dim]{};
    for (int aa = 0; aa != Dim; ++aa)
      for (int cc = 0; cc != Dim; ++cc)
        for (int ii = 0; ii != Dim; ++ii)
          for (int jj = 0; jj != Dim; ++jj)
            s_hat[aa][cc] += tVec(aa, ii) * t_S(ii, jj) * tVec(cc, jj);
 
    int pair_0[sizeSymm];
    int pair_1[sizeSymm];
    V basis_hat[sizeSymm][Dim][Dim];
    // Off-diagonal directions carry 1/2 to recover Cartesian Ddg components.
    for (int ii = 0; ii != Dim; ++ii)
      for (int jj = ii; jj != Dim; ++jj) {
        const int LL = get_sym_index(ii, jj, Dim);
        pair_0[LL] = ii;
        pair_1[LL] = jj;
        for (int aa = 0; aa != Dim; ++aa)
          for (int bb = 0; bb != Dim; ++bb)
            basis_hat[LL][aa][bb] =
                (tVec(aa, ii) * tVec(bb, jj) + tVec(aa, jj) * tVec(bb, ii)) / 2;
      }
 
    FTensor::Ddg<V, Dim, Dim> t_diff_A;
    for (int LL = 0; LL != sizeSymm; ++LL)
      for (int JJ = 0; JJ <= LL; ++JJ) {
        V v = 0;
        for (int aa = 0; aa != Dim; ++aa)
          for (int bb = 0; bb != Dim; ++bb)
            for (int cc = 0; cc != Dim; ++cc)
              v += s_hat[aa][cc] * ddf[aa][bb][cc] *
                   (basis_hat[LL][aa][bb] * basis_hat[JJ][bb][cc] +
                    basis_hat[JJ][aa][bb] * basis_hat[LL][bb][cc]);
 
        t_diff_A(pair_0[LL], pair_1[LL], pair_0[JJ], pair_1[JJ]) = v;
        if (JJ != LL)
          t_diff_A(pair_0[JJ], pair_1[JJ], pair_0[LL], pair_1[LL]) = v;
      }
 
    return t_diff_A;
  }

Member Data Documentation

sizeSymm

template<typename T1 , typename T2 , int Dim>

constexpr int EigenMatrix::EigenMatrixContractedHessianImp< T1, T2, Dim >::sizeSymm = (Dim * (Dim + 1)) / 2

staticconstexpr

Definition at line 804 of file MatrixFunctionTemplate.hpp.

tVal

template<typename T1 , typename T2 , int Dim>

Val& EigenMatrix::EigenMatrixContractedHessianImp< T1, T2, Dim >::tVal

private

Definition at line 927 of file MatrixFunctionTemplate.hpp.

tVec

template<typename T1 , typename T2 , int Dim>

Vec& EigenMatrix::EigenMatrixContractedHessianImp< T1, T2, Dim >::tVec

private

Definition at line 928 of file MatrixFunctionTemplate.hpp.

---

The documentation for this struct was generated from the following file:

• src/ftensor/src/MatrixFunctionTemplate.hpp