adolcfunction.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUFEM_FUNCTIONS_ADOLCFUNCTION_HH
#define DUNE_FUFEM_FUNCTIONS_ADOLCFUNCTION_HH
 
#ifdef HAVE_ADOLC
 
#include <vector>
#include <array>
#include <utility>
 
#include <adolc/adolc.h>
#include <adolc/drivers/drivers.h>
#include <adolc/taping.h>
 
 
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
 
#include <dune/istl/matrix.hh>
 
#include <dune/functions/common/signature.hh>
#include <dune/functions/common/differentiablefunctionfromcallables.hh>
 
 
 
namespace Dune::Fufem {
 
 
namespace Impl {
 
template<class T>
struct AdolCConversion;
 
 
 
template<class T>
using AdolCActiveType = typename AdolCConversion<T>::ActiveType;
 
 
 
template<class T>
struct AdolCConstVectorDataHandle;
 
 
 
template<class T>
struct AdolCMutableMatrixDataHandle;
 
 
 
template<class Signature, class F, class D>
auto adolCTapeFunction(short int tapeTag, F&& f, const D& xp)
{
  using Domain = typename Dune::Functions::SignatureTraits<Signature>::Domain;
  using Range = typename Dune::Functions::SignatureTraits<Signature>::Range;
 
  Range yp;
 
  trace_on(tapeTag);
 
  AdolCActiveType<Domain> x;
  AdolCActiveType<Range> y;
 
  AdolCConversion<Domain>::toActiveType(xp, x);
  y = f(x);
  AdolCConversion<Range>::fromActiveType(y, yp);
 
  trace_off();
 
  return yp;
}
 
 
 
template<class DerivativeSignature, class D, class R>
auto adolCJacobian(short int tapeTag, const D& x, const R&y)
{
  using Domain = typename Dune::Functions::SignatureTraits<DerivativeSignature>::Domain;
  using DerivativeRange = typename Dune::Functions::SignatureTraits<DerivativeSignature>::Range;
 
  auto J = DerivativeRange();
 
  auto x_data_handle = AdolCConstVectorDataHandle<const Domain>(x);
  auto J_data_handle = AdolCMutableMatrixDataHandle<DerivativeRange>(J);
  J_data_handle.resizeTo(y, x);
 
  int n = J_data_handle.N();
  int m = J_data_handle.M();
  assert(m == x_data_handle.size());
 
  jacobian(tapeTag, n, m, x_data_handle.data(), J_data_handle.data());
 
  return J;
}
 
 
 
template<class SecondDerivativeSignature, class D>
auto adolCHessian(short int tapeTag, const D& x)
{
  using SecondDerivativeRange = typename Dune::Functions::SignatureTraits<SecondDerivativeSignature>::Range;
 
  auto H = SecondDerivativeRange();
 
  auto x_data_handle = AdolCConstVectorDataHandle<const D>(x);
  auto H_data_handle = AdolCMutableMatrixDataHandle<SecondDerivativeRange>(H);
  H_data_handle.resizeTo(x, x);
 
  int n = H_data_handle.N();
 
  // Hack around a bug in the AdolC hessian interface which requires
  // that this pointer is mutable (although it hopefully does not
  // modify the pointed value).
  auto x_data_pointer = const_cast<double*>(x_data_handle.data());
  hessian(tapeTag, n, x_data_pointer, H_data_handle.data());
 
  // Copy upper half
  for(auto i : range(n))
    for(auto j : range(i+1,n))
      H_data_handle.data()[i][j] = H_data_handle.data()[j][i];
 
  return H;
}
 
 
 
// *****************************************************************************
// Spezializations of AdolCConversion
// *****************************************************************************
 
template<>
struct AdolCConversion<double>
{
  using PassiveType = double;
  using ActiveType = adouble;
 
  static void toActiveType(const PassiveType& px, ActiveType& x) {
    x <<= px;
  }
 
  static void fromActiveType(ActiveType& y, PassiveType& py) {
     y >>= py;
  }
 
};
 
 
 
template<class K, int n>
struct AdolCConversion<std::array<K,n>>
{
  using PassiveType = std::array<K, n>;
  using ActiveType = std::array<typename AdolCConversion<K>::ActiveType, n>;
 
  static void toActiveType(const PassiveType& px, ActiveType& x) {
    for(auto i : range(px.size()))
      AdolCConversion<K>::toActiveType(px[i], x[i]);
  }
 
  static void fromActiveType(ActiveType& y, PassiveType& py) {
    for(auto i : range(y.size()))
      AdolCConversion<K>::fromActiveType(y[i], py[i]);
  }
 
};
 
 
 
template<class K>
struct AdolCConversion<std::vector<K>>
{
  using PassiveType = std::vector<K>;
  using ActiveType = std::vector<typename AdolCConversion<K>::ActiveType>;
 
  static void toActiveType(const PassiveType& px, ActiveType& x) {
    x.resize(px.size());
    for(auto i : range(px.size()))
      AdolCConversion<K>::toActiveType(px[i], x[i]);
  }
 
  static void fromActiveType(ActiveType& y, PassiveType& py) {
    py.resize(y.size());
    for(auto i : range(y.size()))
      AdolCConversion<K>::fromActiveType(y[i], py[i]);
  }
 
};
 
 
 
template<int n>
struct AdolCConversion<Dune::FieldVector<double,n>>
{
  using PassiveType = Dune::FieldVector<double, n>;
  using ActiveType = Dune::FieldVector<typename AdolCConversion<double>::ActiveType, n>;
 
  static void toActiveType(const PassiveType& px, ActiveType& x) {
    for(auto i : range(px.size()))
      AdolCConversion<double>::toActiveType(px[i], x[i]);
  }
 
  static void fromActiveType(ActiveType& y, PassiveType& py) {
    for(auto i : range(y.size()))
      AdolCConversion<double>::fromActiveType(y[i], py[i]);
  }
 
};
 
 
 
template<class K, int n, int m>
struct AdolCConversion<Dune::FieldMatrix<K,n, m>>
{
  using PassiveType = Dune::FieldMatrix<K, n, m>;
  using ActiveType = Dune::FieldMatrix<typename AdolCConversion<K>::ActiveType, n, m>;
 
  static void toActiveType(const PassiveType& px, ActiveType& x) {
    for(auto i : range(px.N()))
      for(auto j : range(px.M()))
        AdolCConversion<K>::toActiveType(px[i][j], x[i][j]);
  }
 
  static void fromActiveType(ActiveType& y, PassiveType& py) {
    for(auto i : range(y.N()))
      for(auto j : range(y.M()))
        AdolCConversion<K>::fromActiveType(y[i][j], py[i][j]);
  }
 
};
 
 
 
template<class K>
struct AdolCConversion<Dune::Matrix<K>>
{
  using PassiveType = Dune::Matrix<K>;
  using ActiveType = Dune::Matrix<typename AdolCConversion<K>::ActiveType>;
 
  static void toActiveType(const PassiveType& px, ActiveType& x) {
    x.resize(px.N(), px.M());
    for(auto row : range(px.N()))
      for(auto col : range(px.M()))
        AdolCConversion<K>::toActiveType(px[row][col], x[row][col]);
  }
 
  static void fromActiveType(const ActiveType& y, PassiveType& py) {
    py.resize(y.N(), y.M());
    for(auto row : range(y.N()))
      for(auto col : range(y.M()))
        AdolCConversion<K>::fromActiveType(y[row][col], py[row][col]);
  }
 
};
 
 
 
// *****************************************************************************
// Spezializations of AdolCConstVectorDataHandle
// *****************************************************************************
 
template<class T>
struct AdolCConstVectorDataHandle
{
  AdolCConstVectorDataHandle(const T& t) : t_(t) {}
  int size() const { return t_.size(); }
  const double* data() const { return t_.data(); }
  const T& t_;
};
 
template<>
struct AdolCConstVectorDataHandle<const double>
{
  AdolCConstVectorDataHandle(const double& t) : tp_(&t) {}
  int size() const { return 1; }
  const double* data() const { return tp_; }
  const double* tp_;
};
 
 
 
// *****************************************************************************
// Spezializations of AdolCMutableMatrixDataHandle
// *****************************************************************************
 
template<class T>
struct AdolCMutableMatrixDataHandle
{
  AdolCMutableMatrixDataHandle(T& t) : t_(t), data_(t_.data()) {}
  int N() const { return 1; }
  int M() const { return t_.size(); }
  double** data() { return &data_; }
 
  template<class RowVector, class ColVector>
  void resizeTo(const RowVector& rowVector, const ColVector& colVector) {}
 
  T& t_;
  double* data_;
};
 
template<>
struct AdolCMutableMatrixDataHandle<double>
{
  AdolCMutableMatrixDataHandle(double& t) : data_(&t) {}
  int N() const { return 1; }
  int M() const { return 1; }
  double** data() { return &data_; }
 
  template<class RowVector, class ColVector>
  void resizeTo(const RowVector& rowVector, const ColVector& colVector) {}
 
  double* data_;
};
 
template<int n, int m>
struct AdolCMutableMatrixDataHandle<Dune::FieldMatrix<double,n,m>>
{
  AdolCMutableMatrixDataHandle(Dune::FieldMatrix<double,n,m>& t)
  {
    for(auto i : range(n))
      data_[i] = t[i].data();
  }
  int N() const { return n; }
  int M() const { return m; }
  double** data() { return data_.data(); }
 
  template<class RowVector, class ColVector>
  void resizeTo(const RowVector& rowVector, const ColVector& colVector) {}
 
  std::array<double*, n> data_;
};
 
template<>
struct AdolCMutableMatrixDataHandle<std::vector<double>>
{
  using T = std::vector<double>;
  AdolCMutableMatrixDataHandle(T& t) : t_(t), data_(t_.data()) {}
  int N() const { return 1; }
  int M() const { return t_.size(); }
  double** data() { return &data_; }
 
  template<class RowVector>
  void resizeTo(const RowVector& rowVector, const T& colVector) {
    t_.resize(colVector.size());
    data_ = t_.data();
  }
 
  T& t_;
  double* data_;
};
 
template<>
struct AdolCMutableMatrixDataHandle<Dune::Matrix<double>>
{
  AdolCMutableMatrixDataHandle(Dune::Matrix<double>& matrix) :
    matrix_(matrix)
  {
    data_.resize(matrix_.N());
    for(auto i : range(matrix.N()))
      data_[i] = &matrix_[i][0];
  }
  int N() const { return matrix_.N(); }
  int M() const { return matrix_.M(); }
  double** data() { return data_.data(); }
 
  template<class RowVector, class ColVector>
  void resizeTo(const RowVector& rowVector, const ColVector& colVector) {
    matrix_.setSize(rowVector.size(), colVector.size());
    data_.resize(matrix_.N());
    for(auto i : range(matrix_.N()))
      data_[i] = &matrix_[i][0];
  }
 
  Dune::Matrix<double>& matrix_;
  std::vector<double*> data_;
};
 
 
} // namespace Impl (withinin Dune::Fufem::)
 
 
 
 
template<class Signature, template<class> class DerivativeTraits, class F>
auto makeAdolCFunction(const Dune::Functions::SignatureTag<Signature, DerivativeTraits>& signatureTag, F f, short int tapeTag=0)
{
  using Domain = typename Dune::Functions::SignatureTraits<Signature>::Domain;
 
  using DerivativeRange = typename DerivativeTraits<Signature>::Range;
  using DerivativeSignature = DerivativeRange(Domain);
 
  using SecondDerivativeRange = typename DerivativeTraits<DerivativeSignature>::Range;
  using SecondDerivativeSignature = SecondDerivativeRange(Domain);
 
  auto df = [f, tapeTag](const Domain& x) {
    auto y = Dune::Fufem::Impl::adolCTapeFunction<Signature>(tapeTag, f, x);
    return Dune::Fufem::Impl::adolCJacobian<DerivativeSignature>(tapeTag, x, y);
  };
 
  if constexpr (std::is_same_v<SecondDerivativeRange, Dune::Functions::InvalidRange>)
  {
    return Dune::Functions::makeDifferentiableFunctionFromCallables(signatureTag, std::move(f), std::move(df));
  }
  else
  {
    auto ddf = [f, tapeTag](const Domain& x) {
      [[maybe_unused]] auto y = Dune::Fufem::Impl::adolCTapeFunction<Signature>(tapeTag, f, x);
      return Dune::Fufem::Impl::adolCHessian<SecondDerivativeSignature>(tapeTag, x);
    };
    return Dune::Functions::makeDifferentiableFunctionFromCallables(signatureTag, std::move(f), std::move(df), std::move(ddf));
  }
}
 
 
 
} // namespace Dune::Fufem
 
#endif // HAVE_ADOLC
 
#endif // DUNE_FUFEM_FUNCTIONS_ADOLCFUNCTION_HH