Commit 3a8eb3e4 authored by Michael Tesch's avatar Michael Tesch
Browse files

Merge branch 'feature/multidual-eigen' into 'master' 

Add multidual_eigen: Eigen integration for dual<T,N>

See merge request !50
parents 8fa88ba2 3103e20f

CMakeLists.txt

+1 −1
Original line number Diff line number Diff line
@@ -3,7 +3,7 @@ 
#

cmake_minimum_required (VERSION 3.14)

project (cppduals

  VERSION 0.8.0

  VERSION 0.8.1

  LANGUAGES C CXX

  )

include (GNUInstallDirs)

README.md

+9 −1
Original line number Diff line number Diff line
@@ -80,7 +80,8 @@ Headers 
|--------|---------|

| `duals/dual` | Core `dual<T>` class, math functions, IO, complex overloads |

| `duals/multidual` | Multivariate `dual<T,N>` — N partials in one evaluation |

| `duals/dual_eigen` | Eigen NumTraits, type promotion, SIMD packet ops |

| `duals/dual_eigen` | Eigen integration for `dual<T>`: NumTraits, type promotion, SIMD packet ops |

| `duals/multidual_eigen` | Eigen integration for `dual<T,N>`: NumTraits, type promotion, rpart/dpart functors |

| `duals/dual_fmt.h` | Optional fmt library formatters |



Installation
@@ -257,6 +258,13 @@ also licensed under MPL-2.0. 
ChangeLog

=========



v0.8.1

------



- add `duals/multidual_eigen`: Eigen integration for `dual<T,N>`.

- fix C++20 build error: remove illegal `std::is_compound` specializations.

- move coverage generation to local-only, simplify CI pages job.



v0.8.0

------

duals/dual

+2 −3
Original line number Diff line number Diff line
@@ -319,9 +319,8 @@ struct is_arithmetic<duals::dual<T>> : is_arithmetic<T> {}; 


#endif // CPPDUALS_ENABLE_IS_ARITHMETIC



/// Duals are compound types.

template <class T>

struct is_compound<duals::dual<T>> : true_type {};

// is_compound specialization removed — C++20 forbids specializing

// standard type traits for program-defined types.



// Modification of std::numeric_limits<> per

// C++03 17.4.3.1/1, and C++11 18.3.2.3/1.

duals/multidual

+2 −2
Original line number Diff line number Diff line
@@ -223,8 +223,8 @@ template <class T, int N> 
struct is_arithmetic<duals::dual<T,N>> : is_arithmetic<T> {};

#endif



template <class T, int N>

struct is_compound<duals::dual<T,N>> : true_type {};

// is_compound specialization removed — C++20 forbids specializing

// standard type traits for program-defined types.



#define NOMACRO

template <class T, int N>

duals/multidual_eigen

0 → 100644
+217 −0
Original line number Diff line number Diff line 
//===-- duals/multidual_eigen - Eigen integration for dual<T,N> -*- C++ -*-===//

//

// Part of the cppduals project.

// https://gitlab.com/tesch1/cppduals

//

// (c)2026 Michael Tesch. tesch1@gmail.com

//

// This Source Code Form is subject to the terms of the Mozilla

// Public License v. 2.0. If a copy of the MPL was not distributed

// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.



#ifndef CPPDUALS_MULTIDUAL_EIGEN

#define CPPDUALS_MULTIDUAL_EIGEN



#include "multidual"



#ifndef PARSED_BY_DOXYGEN

#include <complex>

#include <Eigen/Core>

#endif



#if !EIGEN_VERSION_AT_LEAST(3, 3, 0)

#error "Eigen too old for cppduals.  Upgrade."

#endif



/** \file       multidual_eigen

    \brief      Eigen integration for multivariate dual<T, int N>



    Include this file to use `duals::dual<T,N>` as a scalar type in

    Eigen matrices.  Provides NumTraits, ScalarBinaryOpTraits (type

    promotion), and matrix-level rpart/dpart extraction functors.



    No SIMD vectorization for N>1 — Eigen falls back to scalar loops.

    For N=1, prefer `duals/dual_eigen` which includes SIMD packet

    specializations.

 */



namespace duals {



/// Unary functor: extract rpart from each element of a dual<T,N> matrix.

template<typename T, int N>

struct CwiseMDRpartOp {

  typedef T result_type;

  EIGEN_DEVICE_FUNC

  EIGEN_STRONG_INLINE T operator()(const dual<T,N> & x) const { return x.rpart(); }

};



/// Unary functor: extract dpart(i) from each element of a dual<T,N> matrix.

template<typename T, int N>

struct CwiseMDDpartOp {

  int idx;

  typedef T result_type;

  CwiseMDDpartOp(int i = 0) : idx(i) {}

  EIGEN_DEVICE_FUNC

  EIGEN_STRONG_INLINE T operator()(const dual<T,N> & x) const { return x.dpart(idx); }

};



/// Unary functor: dual-conjugate each element of a dual<T,N> matrix.

template<typename T, int N>

struct CwiseMDDconjOp {

  typedef dual<T,N> result_type;

  EIGEN_DEVICE_FUNC

  EIGEN_STRONG_INLINE dual<T,N> operator()(const dual<T,N> & x) const { return dconj(x); }

};



/// Extract the "real part" of a dual<T,N>-valued matrix.

template <typename Derived>

auto rpart(const Eigen::EigenBase<Derived> & x)

  -> decltype(x.derived().unaryExpr(CwiseMDRpartOp<

    typename Derived::Scalar::value_type, Derived::Scalar::num_vars>()))

{

  using S = typename Derived::Scalar;

  return x.derived().unaryExpr(CwiseMDRpartOp<typename S::value_type, S::num_vars>());

}



/// Extract the i-th dual part of a dual<T,N>-valued matrix.

template <typename Derived>

auto dpart(const Eigen::EigenBase<Derived> & x, int i = 0)

  -> decltype(x.derived().unaryExpr(CwiseMDDpartOp<

    typename Derived::Scalar::value_type, Derived::Scalar::num_vars>(i)))

{

  using S = typename Derived::Scalar;

  return x.derived().unaryExpr(CwiseMDDpartOp<typename S::value_type, S::num_vars>(i));

}



/// Dual-conjugate a dual<T,N>-valued matrix.

template <typename Derived>

auto dconj(const Eigen::EigenBase<Derived> & x)

  -> decltype(x.derived().unaryExpr(CwiseMDDconjOp<

    typename Derived::Scalar::value_type, Derived::Scalar::num_vars>()))

{

  using S = typename Derived::Scalar;

  return x.derived().unaryExpr(CwiseMDDconjOp<typename S::value_type, S::num_vars>());

}



} // namespace duals



namespace Eigen {



// ===================================================================

// NumTraits<dual<T,N>> — tells Eigen how to work with dual scalars

// ===================================================================



template<typename T, int N>

struct NumTraits<duals::dual<T,N>> : GenericNumTraits<T>

{

  typedef duals::dual<T,N> Real;

  typedef duals::dual<T,N> Literal;

  typedef duals::dual<T,N> Nested;



  enum {

    IsInteger           = NumTraits<T>::IsInteger,

    IsSigned            = NumTraits<T>::IsSigned,

    IsComplex           = 0,

    RequireInitialization = 1,

    ReadCost            = (N + 1) * NumTraits<T>::ReadCost,

    AddCost             = (N + 1) * NumTraits<T>::AddCost,

    MulCost             = (2 * N + 1) * NumTraits<T>::MulCost + N * NumTraits<T>::AddCost

  };



  EIGEN_DEVICE_FUNC

  static inline Real epsilon()          { return Real(NumTraits<T>::epsilon()); }

  EIGEN_DEVICE_FUNC

  static inline Real dummy_precision()  { return NumTraits<T>::dummy_precision(); }

  EIGEN_DEVICE_FUNC

  static inline Real highest()          { return Real(NumTraits<T>::highest()); }

  EIGEN_DEVICE_FUNC

  static inline Real lowest()           { return Real(NumTraits<T>::lowest()); }

  EIGEN_DEVICE_FUNC

  static inline int digits10()          { return NumTraits<T>::digits10(); }

};



// ===================================================================

// ScalarBinaryOpTraits — type promotion for mixed operations

// ===================================================================



#if !defined(CPPDUALS_NO_EIGEN_PROMOTION)



// dual<T,N> op dual<T,N>

template<typename T, int N, typename BinaryOp>

struct ScalarBinaryOpTraits<duals::dual<T,N>, duals::dual<T,N>, BinaryOp>

  : public duals::can_promote<duals::dual<T,N>, duals::dual<T,N>>::wrap {};



// dual<T,N> op T

template<typename T, int N, typename BinaryOp>

struct ScalarBinaryOpTraits<duals::dual<T,N>, T, BinaryOp>

  : public duals::can_promote<duals::dual<T,N>, T>::wrap {};



// T op dual<T,N>

template<typename T, int N, typename BinaryOp>

struct ScalarBinaryOpTraits<T, duals::dual<T,N>, BinaryOp>

  : public duals::can_promote<duals::dual<T,N>, T>::wrap {};



#ifndef PARSED_BY_DOXYGEN



// complex<dual<T,N>> op scalar promotions

#define CPPDUALS_MDN_SBOTS_REALS(U)                                             \

  template<typename T, int N, typename BinaryOp>                                \

  struct ScalarBinaryOpTraits<std::complex<duals::dual<T,N>>, U, BinaryOp>      \

    : public duals::can_promote<std::complex<duals::dual<T,N>>, U>::wrap {};    \

  template<typename T, int N, typename BinaryOp>                                \

  struct ScalarBinaryOpTraits<U, std::complex<duals::dual<T,N>>, BinaryOp>      \

    : public duals::can_promote<std::complex<duals::dual<T,N>>, U>::wrap {}

CPPDUALS_MDN_SBOTS_REALS(int);

CPPDUALS_MDN_SBOTS_REALS(float);

CPPDUALS_MDN_SBOTS_REALS(double);



// complex<dual<T,N>> op complex<T>

template<typename T, int N, typename BinaryOp>

struct ScalarBinaryOpTraits<std::complex<duals::dual<T,N>>, std::complex<T>, BinaryOp>

  : public duals::can_promote<std::complex<duals::dual<T,N>>, std::complex<T>>::wrap {};

template<typename T, int N, typename BinaryOp>

struct ScalarBinaryOpTraits<std::complex<T>, std::complex<duals::dual<T,N>>, BinaryOp>

  : public duals::can_promote<std::complex<duals::dual<T,N>>, std::complex<T>>::wrap {};



#endif // PARSED_BY_DOXYGEN



#endif // CPPDUALS_NO_EIGEN_PROMOTION



// ===================================================================

// Eigen::internal — real_impl, random, numext

// ===================================================================



#ifndef PARSED_BY_DOXYGEN



namespace numext {

using duals::rpart;

using duals::dpart;

using duals::dconj;

}



namespace internal {



template<typename T, int N>

struct real_impl<duals::dual<T,N>>

{

  typedef T RealScalar;

  EIGEN_DEVICE_FUNC

  static inline T run(const duals::dual<T,N>& x) { return x.rpart(); }

};



template<typename T, int N>

struct scalar_random_op<duals::dual<T,N>>

{

  EIGEN_EMPTY_STRUCT_CTOR(scalar_random_op)

  inline const duals::dual<T,N> operator() () const {

    return duals::random<duals::dual<T,N>>();

  }

};



} // namespace internal



#endif // PARSED_BY_DOXYGEN



} // namespace Eigen



#endif // CPPDUALS_MULTIDUAL_EIGEN