mirror of
https://github.com/NGSolve/netgen.git
synced 2024-11-11 16:49:16 +05:00
1132 lines
30 KiB
C++
1132 lines
30 KiB
C++
#ifndef FILE_AUTODIFF
|
|
#define FILE_AUTODIFF
|
|
|
|
/**************************************************************************/
|
|
/* File: autodiff.hpp */
|
|
/* Author: Joachim Schoeberl */
|
|
/* Date: 24. Oct. 02 */
|
|
/**************************************************************************/
|
|
|
|
namespace ngcore
|
|
{
|
|
using ngcore::IfPos;
|
|
|
|
// Automatic differentiation datatype
|
|
|
|
template <int D, typename SCAL = double> class AutoDiffRec;
|
|
|
|
|
|
/**
|
|
Datatype for automatic differentiation.
|
|
Contains function value and D derivatives. Algebraic
|
|
operations are overloaded by using product-rule etc. etc.
|
|
**/
|
|
template <int D, typename SCAL = double>
|
|
class AutoDiffVec
|
|
{
|
|
SCAL val;
|
|
SCAL dval[D?D:1];
|
|
public:
|
|
|
|
typedef AutoDiffVec<D,SCAL> TELEM;
|
|
typedef SCAL TSCAL;
|
|
|
|
|
|
/// elements are undefined
|
|
// NETGEN_INLINE AutoDiffVec () throw() { };
|
|
AutoDiffVec() = default;
|
|
// { val = 0; for (int i = 0; i < D; i++) dval[i] = 0; } // !
|
|
|
|
/// copy constructor
|
|
AutoDiffVec (const AutoDiffVec & ad2) = default;
|
|
/*
|
|
NETGEN_INLINE AutoDiffVec (const AutoDiffVec & ad2) throw()
|
|
{
|
|
val = ad2.val;
|
|
for (int i = 0; i < D; i++)
|
|
dval[i] = ad2.dval[i];
|
|
}
|
|
*/
|
|
/// initial object with constant value
|
|
NETGEN_INLINE AutoDiffVec (SCAL aval) throw()
|
|
{
|
|
val = aval;
|
|
for (int i = 0; i < D; i++)
|
|
dval[i] = 0;
|
|
}
|
|
|
|
/// init object with (val, e_diffindex)
|
|
NETGEN_INLINE AutoDiffVec (SCAL aval, int diffindex) throw()
|
|
{
|
|
val = aval;
|
|
for (int i = 0; i < D; i++)
|
|
dval[i] = 0;
|
|
dval[diffindex] = 1;
|
|
}
|
|
|
|
NETGEN_INLINE AutoDiffVec (SCAL aval, const SCAL * grad)
|
|
{
|
|
val = aval;
|
|
LoadGradient (grad);
|
|
}
|
|
|
|
/// assign constant value
|
|
NETGEN_INLINE AutoDiffVec & operator= (SCAL aval) throw()
|
|
{
|
|
val = aval;
|
|
for (int i = 0; i < D; i++)
|
|
dval[i] = 0;
|
|
return *this;
|
|
}
|
|
|
|
AutoDiffVec & operator= (const AutoDiffVec & ad2) = default;
|
|
|
|
/// returns value
|
|
NETGEN_INLINE SCAL Value() const throw() { return val; }
|
|
|
|
/// returns partial derivative
|
|
NETGEN_INLINE SCAL DValue (int i) const throw() { return dval[i]; }
|
|
|
|
///
|
|
NETGEN_INLINE void StoreGradient (SCAL * p) const
|
|
{
|
|
for (int i = 0; i < D; i++)
|
|
p[i] = dval[i];
|
|
}
|
|
|
|
NETGEN_INLINE void LoadGradient (const SCAL * p)
|
|
{
|
|
for (int i = 0; i < D; i++)
|
|
dval[i] = p[i];
|
|
}
|
|
|
|
/// access value
|
|
NETGEN_INLINE SCAL & Value() throw() { return val; }
|
|
|
|
/// accesses partial derivative
|
|
NETGEN_INLINE SCAL & DValue (int i) throw() { return dval[i]; }
|
|
};
|
|
|
|
|
|
//@{ AutoDiffVec helper functions.
|
|
|
|
/// prints AutoDiffVec
|
|
template<int D, typename SCAL>
|
|
inline ostream & operator<< (ostream & ost, const AutoDiffVec<D,SCAL> & x)
|
|
{
|
|
ost << x.Value() << ", D = ";
|
|
for (int i = 0; i < D; i++)
|
|
ost << x.DValue(i) << " ";
|
|
return ost;
|
|
}
|
|
|
|
/// AutoDiffVec plus AutoDiffVec
|
|
template<int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> operator+ (const AutoDiffVec<D,SCAL> & x, const AutoDiffVec<D,SCAL> & y) throw()
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value () = x.Value()+y.Value();
|
|
// AutoDiffVec<D,SCAL> res(x.Value()+y.Value());
|
|
for (int i = 0; i < D; i++)
|
|
res.DValue(i) = x.DValue(i) + y.DValue(i);
|
|
return res;
|
|
}
|
|
|
|
|
|
/// AutoDiffVec minus AutoDiffVec
|
|
template<int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> operator- (const AutoDiffVec<D,SCAL> & x, const AutoDiffVec<D,SCAL> & y) throw()
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = x.Value()-y.Value();
|
|
// AutoDiffVec<D,SCAL> res (x.Value()-y.Value());
|
|
for (int i = 0; i < D; i++)
|
|
res.DValue(i) = x.DValue(i) - y.DValue(i);
|
|
return res;
|
|
}
|
|
|
|
/// double plus AutoDiffVec
|
|
template<int D, typename SCAL, typename SCAL2,
|
|
typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> operator+ (SCAL2 x, const AutoDiffVec<D,SCAL> & y) throw()
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = x+y.Value();
|
|
for (int i = 0; i < D; i++)
|
|
res.DValue(i) = y.DValue(i);
|
|
return res;
|
|
}
|
|
|
|
/// AutoDiffVec plus double
|
|
template<int D, typename SCAL, typename SCAL2,
|
|
typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> operator+ (const AutoDiffVec<D,SCAL> & y, SCAL2 x) throw()
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = x+y.Value();
|
|
for (int i = 0; i < D; i++)
|
|
res.DValue(i) = y.DValue(i);
|
|
return res;
|
|
}
|
|
|
|
|
|
/// minus AutoDiffVec
|
|
template<int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> operator- (const AutoDiffVec<D,SCAL> & x) throw()
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = -x.Value();
|
|
for (int i = 0; i < D; i++)
|
|
res.DValue(i) = -x.DValue(i);
|
|
return res;
|
|
}
|
|
|
|
/// AutoDiffVec minus double
|
|
template<int D, typename SCAL, typename SCAL2,
|
|
typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> operator- (const AutoDiffVec<D,SCAL> & x, SCAL2 y) throw()
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = x.Value()-y;
|
|
for (int i = 0; i < D; i++)
|
|
res.DValue(i) = x.DValue(i);
|
|
return res;
|
|
}
|
|
|
|
///
|
|
template<int D, typename SCAL, typename SCAL2,
|
|
typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> operator- (SCAL2 x, const AutoDiffVec<D,SCAL> & y) throw()
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = x-y.Value();
|
|
for (int i = 0; i < D; i++)
|
|
res.DValue(i) = -y.DValue(i);
|
|
return res;
|
|
}
|
|
|
|
|
|
/// double times AutoDiffVec
|
|
template<int D, typename SCAL, typename SCAL2,
|
|
typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> operator* (SCAL2 x, const AutoDiffVec<D,SCAL> & y) throw()
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = x*y.Value();
|
|
for (int i = 0; i < D; i++)
|
|
res.DValue(i) = x*y.DValue(i);
|
|
return res;
|
|
}
|
|
|
|
/// AutoDiffVec times double
|
|
template<int D, typename SCAL, typename SCAL2,
|
|
typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
|
|
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> operator* (const AutoDiffVec<D,SCAL> & y, SCAL2 x) throw()
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = x*y.Value();
|
|
for (int i = 0; i < D; i++)
|
|
res.DValue(i) = x*y.DValue(i);
|
|
return res;
|
|
}
|
|
|
|
/// AutoDiffVec times AutoDiffVec
|
|
template<int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> operator* (const AutoDiffVec<D,SCAL> & x, const AutoDiffVec<D,SCAL> & y) throw()
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
SCAL hx = x.Value();
|
|
SCAL hy = y.Value();
|
|
|
|
res.Value() = hx*hy;
|
|
for (int i = 0; i < D; i++)
|
|
res.DValue(i) = hx*y.DValue(i) + hy*x.DValue(i);
|
|
|
|
return res;
|
|
}
|
|
|
|
/// AutoDiffVec times AutoDiffVec
|
|
using ngcore::sqr;
|
|
template<int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> sqr (const AutoDiffVec<D,SCAL> & x) throw()
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
SCAL hx = x.Value();
|
|
res.Value() = hx*hx;
|
|
hx *= 2;
|
|
for (int i = 0; i < D; i++)
|
|
res.DValue(i) = hx*x.DValue(i);
|
|
return res;
|
|
}
|
|
|
|
/// Inverse of AutoDiffVec
|
|
template<int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> Inv (const AutoDiffVec<D,SCAL> & x)
|
|
{
|
|
AutoDiffVec<D,SCAL> res(1.0 / x.Value());
|
|
for (int i = 0; i < D; i++)
|
|
res.DValue(i) = -x.DValue(i) / (x.Value() * x.Value());
|
|
return res;
|
|
}
|
|
|
|
|
|
/// AutoDiffVec div AutoDiffVec
|
|
template<int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> operator/ (const AutoDiffVec<D,SCAL> & x, const AutoDiffVec<D,SCAL> & y)
|
|
{
|
|
return x * Inv (y);
|
|
}
|
|
|
|
/// AutoDiffVec div double
|
|
template<int D, typename SCAL, typename SCAL2,
|
|
typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> operator/ (const AutoDiffVec<D,SCAL> & x, SCAL2 y)
|
|
{
|
|
return (1.0/y) * x;
|
|
}
|
|
|
|
/// double div AutoDiffVec
|
|
template<int D, typename SCAL, typename SCAL2,
|
|
typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> operator/ (SCAL2 x, const AutoDiffVec<D,SCAL> & y)
|
|
{
|
|
return x * Inv(y);
|
|
}
|
|
|
|
|
|
|
|
|
|
template <int D, typename SCAL, typename SCAL2>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> & operator+= (AutoDiffVec<D,SCAL> & x, SCAL2 y) throw()
|
|
{
|
|
x.Value() += y;
|
|
return x;
|
|
}
|
|
|
|
|
|
///
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> & operator+= (AutoDiffVec<D,SCAL> & x, AutoDiffVec<D,SCAL> y)
|
|
{
|
|
x.Value() += y.Value();
|
|
for (int i = 0; i < D; i++)
|
|
x.DValue(i) += y.DValue(i);
|
|
return x;
|
|
}
|
|
|
|
///
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> & operator-= (AutoDiffVec<D,SCAL> & x, AutoDiffVec<D,SCAL> y)
|
|
{
|
|
x.Value() -= y.Value();
|
|
for (int i = 0; i < D; i++)
|
|
x.DValue(i) -= y.DValue(i);
|
|
return x;
|
|
|
|
}
|
|
|
|
template <int D, typename SCAL, typename SCAL2>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> & operator-= (AutoDiffVec<D,SCAL> & x, SCAL2 y)
|
|
{
|
|
x.Value() -= y;
|
|
return x;
|
|
}
|
|
|
|
///
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> & operator*= (AutoDiffVec<D,SCAL> & x, AutoDiffVec<D,SCAL> y)
|
|
{
|
|
for (int i = 0; i < D; i++)
|
|
x.DValue(i) = x.DValue(i)*y.Value() + x.Value() * y.DValue(i);
|
|
x.Value() *= y.Value();
|
|
return x;
|
|
}
|
|
|
|
///
|
|
template <int D, typename SCAL, typename SCAL2>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> & operator*= (AutoDiffVec<D,SCAL> & x, SCAL2 y)
|
|
{
|
|
x.Value() *= y;
|
|
for (int i = 0; i < D; i++)
|
|
x.DValue(i) *= y;
|
|
return x;
|
|
}
|
|
|
|
///
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> & operator/= (AutoDiffVec<D,SCAL> & x, SCAL y)
|
|
{
|
|
SCAL iy = 1.0 / y;
|
|
x.Value() *= iy;
|
|
for (int i = 0; i < D; i++)
|
|
x.DValue(i) *= iy;
|
|
return x;
|
|
}
|
|
|
|
|
|
|
|
|
|
///
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE bool operator== (AutoDiffVec<D,SCAL> x, SCAL val2)
|
|
{
|
|
return x.Value() == val2;
|
|
}
|
|
|
|
///
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE bool operator!= (AutoDiffVec<D,SCAL> x, SCAL val2) throw()
|
|
{
|
|
return x.Value() != val2;
|
|
}
|
|
|
|
///
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE bool operator< (AutoDiffVec<D,SCAL> x, SCAL val2) throw()
|
|
{
|
|
return x.Value() < val2;
|
|
}
|
|
|
|
///
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE bool operator> (AutoDiffVec<D,SCAL> x, SCAL val2) throw()
|
|
{
|
|
return x.Value() > val2;
|
|
}
|
|
|
|
|
|
|
|
|
|
template<int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> fabs (const AutoDiffVec<D,SCAL> & x)
|
|
{
|
|
double abs = fabs (x.Value());
|
|
AutoDiffVec<D,SCAL> res( abs );
|
|
if (abs != 0.0)
|
|
for (int i = 0; i < D; i++)
|
|
res.DValue(i) = x.Value()*x.DValue(i) / abs;
|
|
else
|
|
for (int i = 0; i < D; i++)
|
|
res.DValue(i) = 0.0;
|
|
return res;
|
|
}
|
|
|
|
using std::sqrt;
|
|
template<int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> sqrt (const AutoDiffVec<D,SCAL> & x)
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = sqrt(x.Value());
|
|
for (int j = 0; j < D; j++)
|
|
res.DValue(j) = 0.5 / res.Value() * x.DValue(j);
|
|
return res;
|
|
}
|
|
|
|
using std::log;
|
|
template <int D, typename SCAL>
|
|
AutoDiffVec<D,SCAL> log (AutoDiffVec<D,SCAL> x)
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = log(x.Value());
|
|
for (int k = 0; k < D; k++)
|
|
res.DValue(k) = x.DValue(k) / x.Value();
|
|
return res;
|
|
}
|
|
|
|
using std::exp;
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> exp (AutoDiffVec<D,SCAL> x)
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = exp(x.Value());
|
|
for (int k = 0; k < D; k++)
|
|
res.DValue(k) = x.DValue(k) * res.Value();
|
|
return res;
|
|
}
|
|
|
|
using std::pow;
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> pow (AutoDiffVec<D,SCAL> x, AutoDiffVec<D,SCAL> y )
|
|
{
|
|
return exp(log(x)*y);
|
|
}
|
|
|
|
using std::sin;
|
|
/*
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> sin (AutoDiffVec<D,SCAL> x)
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = sin(x.Value());
|
|
SCAL c = cos(x.Value());
|
|
for (int k = 0; k < D; k++)
|
|
res.DValue(k) = x.DValue(k) * c;
|
|
return res;
|
|
}
|
|
*/
|
|
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> sin (AutoDiffVec<D,SCAL> x)
|
|
{
|
|
return sin(AutoDiffRec<D,SCAL>(x));
|
|
}
|
|
|
|
using std::cos;
|
|
/*
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> cos (AutoDiffVec<D,SCAL> x)
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = cos(x.Value());
|
|
SCAL ms = -sin(x.Value());
|
|
for (int k = 0; k < D; k++)
|
|
res.DValue(k) = x.DValue(k) * ms;
|
|
return res;
|
|
}
|
|
*/
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> cos (AutoDiffVec<D,SCAL> x)
|
|
{
|
|
return cos(AutoDiffRec<D,SCAL>(x));
|
|
}
|
|
|
|
using std::tan;
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> tan (AutoDiffVec<D,SCAL> x)
|
|
{ return sin(x) / cos(x); }
|
|
|
|
using std::sinh;
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> sinh (AutoDiffVec<D,SCAL> x)
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = sinh(x.Value());
|
|
SCAL ch = cosh(x.Value());
|
|
for (int k = 0; k < D; k++)
|
|
res.DValue(k) = x.DValue(k) * ch;
|
|
return res;
|
|
}
|
|
|
|
using std::cosh;
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> cosh (AutoDiffVec<D,SCAL> x)
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = cosh(x.Value());
|
|
SCAL sh = sinh(x.Value());
|
|
for (int k = 0; k < D; k++)
|
|
res.DValue(k) = x.DValue(k) * sh;
|
|
return res;
|
|
}
|
|
|
|
using std::erf;
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> erf (AutoDiffVec<D,SCAL> x)
|
|
{
|
|
return erf(AutoDiffRec<D,SCAL>(x));
|
|
}
|
|
|
|
using std::floor;
|
|
template<int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> floor (const AutoDiffVec<D,SCAL> & x)
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = floor(x.Value());
|
|
for (int j = 0; j < D; j++)
|
|
res.DValue(j) = 0.0;
|
|
return res;
|
|
}
|
|
|
|
using std::ceil;
|
|
template<int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> ceil (const AutoDiffVec<D,SCAL> & x)
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = ceil(x.Value());
|
|
for (int j = 0; j < D; j++)
|
|
res.DValue(j) = 0.0;
|
|
return res;
|
|
}
|
|
|
|
|
|
using std::atan;
|
|
/*
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> atan (AutoDiffVec<D,SCAL> x)
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
SCAL a = atan(x.Value());
|
|
res.Value() = a;
|
|
for (int k = 0; k < D; k++)
|
|
res.DValue(k) = x.DValue(k)/(1+x.Value()*x.Value()) ;
|
|
return res;
|
|
}
|
|
*/
|
|
template <int D, typename SCAL>
|
|
AutoDiffVec<D,SCAL> atan (AutoDiffVec<D,SCAL> x)
|
|
{
|
|
return atan (AutoDiffRec<D,SCAL> (x));
|
|
}
|
|
|
|
using std::atan2;
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> atan2 (AutoDiffVec<D,SCAL> x, AutoDiffVec<D,SCAL> y)
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
SCAL a = atan2(x.Value(), y.Value());
|
|
res.Value() = a;
|
|
for (int k = 0; k < D; k++)
|
|
res.DValue(k) = (x.Value()*y.DValue(k)-y.Value()*x.DValue(k))/(y.Value()*y.Value()+x.Value()*x.Value());
|
|
return res;
|
|
}
|
|
|
|
|
|
using std::acos;
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> acos (AutoDiffVec<D,SCAL> x)
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
SCAL a = acos(x.Value());
|
|
res.Value() = a;
|
|
SCAL da = -1 / sqrt(1-x.Value()*x.Value());
|
|
for (int k = 0; k < D; k++)
|
|
res.DValue(k) = x.DValue(k)*da;
|
|
return res;
|
|
}
|
|
|
|
|
|
using std::asin;
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> asin (AutoDiffVec<D,SCAL> x)
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
SCAL a = asin(x.Value());
|
|
res.Value() = a;
|
|
SCAL da = 1 / sqrt(1-x.Value()*x.Value());
|
|
for (int k = 0; k < D; k++)
|
|
res.DValue(k) = x.DValue(k)*da;
|
|
return res;
|
|
}
|
|
|
|
|
|
|
|
|
|
template <int D, typename SCAL, typename TB, typename TC>
|
|
auto IfPos (AutoDiffVec<D,SCAL> a, TB b, TC c) // -> decltype(IfPos (a.Value(), b, c))
|
|
{
|
|
return IfPos (a.Value(), b, c);
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffVec<D,SCAL> b, AutoDiffVec<D,SCAL> c)
|
|
{
|
|
AutoDiffVec<D,SCAL> res;
|
|
res.Value() = IfPos (a, b.Value(), c.Value());
|
|
for (int j = 0; j < D; j++)
|
|
res.DValue(j) = IfPos (a, b.DValue(j), c.DValue(j));
|
|
return res;
|
|
}
|
|
|
|
template <int D, typename SCAL, typename TC>
|
|
NETGEN_INLINE AutoDiffVec<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffVec<D,SCAL> b, TC c)
|
|
{
|
|
return IfPos (a, b, AutoDiffVec<D,SCAL> (c));
|
|
}
|
|
|
|
//@}
|
|
|
|
|
|
|
|
template <int D, typename SCAL>
|
|
class AutoDiffRec
|
|
{
|
|
AutoDiffRec<D-1, SCAL> rec;
|
|
SCAL last;
|
|
|
|
public:
|
|
NETGEN_INLINE AutoDiffRec () = default;
|
|
NETGEN_INLINE AutoDiffRec (const AutoDiffRec &) = default;
|
|
NETGEN_INLINE AutoDiffRec (AutoDiffRec<D-1,SCAL> _rec, SCAL _last) : rec(_rec), last(_last) { ; }
|
|
NETGEN_INLINE AutoDiffRec & operator= (const AutoDiffRec &) = default;
|
|
|
|
NETGEN_INLINE AutoDiffRec (SCAL aval) : rec(aval), last(0.0) { ; }
|
|
NETGEN_INLINE AutoDiffRec (SCAL aval, int diffindex) : rec(aval, diffindex), last((diffindex==D-1) ? 1.0 : 0.0) { ; }
|
|
NETGEN_INLINE AutoDiffRec (SCAL aval, const SCAL * grad)
|
|
: rec(aval, grad), last(grad[D-1]) { }
|
|
|
|
NETGEN_INLINE AutoDiffRec (const AutoDiffVec<D,SCAL> & ad)
|
|
{
|
|
Value() = ad.Value();
|
|
for (int i = 0; i < D; i++)
|
|
DValue(i) = ad.DValue(i);
|
|
}
|
|
|
|
NETGEN_INLINE AutoDiffRec & operator= (SCAL aval) { rec = aval; last = 0.0; return *this; }
|
|
NETGEN_INLINE SCAL Value() const { return rec.Value(); }
|
|
NETGEN_INLINE SCAL DValue(int i) const { return (i == D-1) ? last : rec.DValue(i); }
|
|
NETGEN_INLINE SCAL & Value() { return rec.Value(); }
|
|
NETGEN_INLINE SCAL & DValue(int i) { return (i == D-1) ? last : rec.DValue(i); }
|
|
NETGEN_INLINE auto Rec() const { return rec; }
|
|
NETGEN_INLINE auto Last() const { return last; }
|
|
NETGEN_INLINE auto & Rec() { return rec; }
|
|
NETGEN_INLINE auto & Last() { return last; }
|
|
NETGEN_INLINE operator AutoDiffVec<D,SCAL> () const
|
|
{
|
|
AutoDiffVec<D,SCAL> res(Value());
|
|
for (int i = 0; i < D; i++)
|
|
res.DValue(i) = DValue(i);
|
|
return res;
|
|
}
|
|
};
|
|
|
|
template<int D, typename SCAL>
|
|
ostream & operator<< (ostream & ost, AutoDiffRec<D,SCAL> ad)
|
|
{
|
|
return ost << AutoDiffVec<D,SCAL> (ad);
|
|
}
|
|
|
|
template <typename SCAL>
|
|
class AutoDiffRec<0,SCAL>
|
|
{
|
|
SCAL val;
|
|
public:
|
|
NETGEN_INLINE AutoDiffRec () = default;
|
|
NETGEN_INLINE AutoDiffRec (const AutoDiffRec &) = default;
|
|
NETGEN_INLINE AutoDiffRec (SCAL _val) : val(_val) { ; }
|
|
NETGEN_INLINE AutoDiffRec (SCAL _val, SCAL /* _dummylast */) : val(_val) { ; }
|
|
NETGEN_INLINE AutoDiffRec (SCAL aval, const SCAL * /* grad */)
|
|
: val(aval) { }
|
|
|
|
NETGEN_INLINE AutoDiffRec & operator= (const AutoDiffRec &) = default;
|
|
NETGEN_INLINE AutoDiffRec & operator= (SCAL aval) { val = aval; return *this; }
|
|
|
|
NETGEN_INLINE SCAL Value() const { return val; }
|
|
NETGEN_INLINE SCAL DValue(int /* i */) const { return SCAL(0); }
|
|
NETGEN_INLINE SCAL & Value() { return val; }
|
|
// SCAL & DValue(int i) { return val; }
|
|
NETGEN_INLINE auto Rec() const { return val; }
|
|
NETGEN_INLINE auto Last() const { return SCAL(0); }
|
|
NETGEN_INLINE auto & Rec() { return val; }
|
|
NETGEN_INLINE auto & Last() { return val; }
|
|
NETGEN_INLINE operator AutoDiffVec<0,SCAL> () const { return AutoDiffVec<0,SCAL>(); }
|
|
};
|
|
|
|
|
|
template <typename SCAL>
|
|
class AutoDiffRec<1,SCAL>
|
|
{
|
|
SCAL val;
|
|
SCAL last;
|
|
public:
|
|
NETGEN_INLINE AutoDiffRec () = default;
|
|
NETGEN_INLINE AutoDiffRec (const AutoDiffRec &) = default;
|
|
NETGEN_INLINE AutoDiffRec (SCAL _val) : val(_val), last(0.0) { ; }
|
|
NETGEN_INLINE AutoDiffRec (SCAL _val, SCAL _last) : val(_val), last(_last) { ; }
|
|
NETGEN_INLINE AutoDiffRec (SCAL aval, int diffindex) : val(aval), last((diffindex==0) ? 1.0 : 0.0) { ; }
|
|
NETGEN_INLINE AutoDiffRec (SCAL aval, const SCAL * grad)
|
|
: val(aval), last(grad[0]) { }
|
|
|
|
NETGEN_INLINE AutoDiffRec (const AutoDiffVec<1,SCAL> & ad)
|
|
{
|
|
Value() = ad.Value();
|
|
DValue(0) = ad.DValue(0);
|
|
}
|
|
|
|
NETGEN_INLINE AutoDiffRec & operator= (const AutoDiffRec &) = default;
|
|
NETGEN_INLINE AutoDiffRec & operator= (SCAL aval) { val = aval; last = 0.0; return *this; }
|
|
|
|
NETGEN_INLINE SCAL Value() const { return val; }
|
|
NETGEN_INLINE SCAL DValue(int /* i */) const { return last; }
|
|
NETGEN_INLINE SCAL & Value() { return val; }
|
|
NETGEN_INLINE SCAL & DValue(int /* i */) { return last; }
|
|
NETGEN_INLINE auto Rec() const { return val; }
|
|
NETGEN_INLINE auto Last() const { return last; }
|
|
NETGEN_INLINE auto & Rec() { return val; }
|
|
NETGEN_INLINE auto & Last() { return last; }
|
|
|
|
NETGEN_INLINE operator AutoDiffVec<1,SCAL> () const
|
|
{
|
|
AutoDiffVec<1,SCAL> res(Value());
|
|
res.DValue(0) = DValue(0);
|
|
return res;
|
|
}
|
|
};
|
|
|
|
template <int D, typename SCAL, typename SCAL2,
|
|
typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> operator+ (SCAL2 a, AutoDiffRec<D,SCAL> b)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (a+b.Rec(), b.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL, typename SCAL2,
|
|
typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> operator+ (AutoDiffRec<D,SCAL> a, SCAL2 b)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (a.Rec()+b, a.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> operator+ (AutoDiffRec<D,SCAL> a, AutoDiffRec<D,SCAL> b)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (a.Rec()+b.Rec(), a.Last()+b.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL, typename SCAL2,
|
|
typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> operator- (SCAL2 b, AutoDiffRec<D,SCAL> a)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (b-a.Rec(), -a.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL, typename SCAL2,
|
|
typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> operator- (AutoDiffRec<D,SCAL> a, SCAL2 b)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (a.Rec()-b, a.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> operator- (AutoDiffRec<D,SCAL> a, AutoDiffRec<D,SCAL> b)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (a.Rec()-b.Rec(), a.Last()-b.Last());
|
|
}
|
|
|
|
/// minus AutoDiff
|
|
template<int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> operator- (AutoDiffRec<D,SCAL> a)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (-a.Rec(), -a.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> operator* (AutoDiffRec<D,SCAL> a, AutoDiffRec<D,SCAL> b)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (a.Rec()*b.Rec(), a.Value()*b.Last()+b.Value()*a.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL, typename SCAL1,
|
|
typename std::enable_if<std::is_convertible<SCAL1,SCAL>::value, int>::type = 0>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> operator* (AutoDiffRec<D,SCAL> b, SCAL1 a)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (a*b.Rec(), a*b.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL, typename SCAL1,
|
|
typename std::enable_if<std::is_convertible<SCAL1,SCAL>::value, int>::type = 0>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> operator* (SCAL1 a, AutoDiffRec<D,SCAL> b)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (a*b.Rec(), a*b.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> & operator+= (AutoDiffRec<D,SCAL> & a, AutoDiffRec<D,SCAL> b)
|
|
{
|
|
a.Rec() += b.Rec();
|
|
a.Last() += b.Last();
|
|
return a;
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> & operator-= (AutoDiffRec<D,SCAL> & a, double b)
|
|
{
|
|
a.Rec() -= b;
|
|
return a;
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> & operator-= (AutoDiffRec<D,SCAL> & a, AutoDiffRec<D,SCAL> b)
|
|
{
|
|
a.Rec() -= b.Rec();
|
|
a.Last() -= b.Last();
|
|
return a;
|
|
}
|
|
|
|
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> & operator*= (AutoDiffRec<D,SCAL> & a, AutoDiffRec<D,SCAL> b)
|
|
{
|
|
a = a*b;
|
|
return a;
|
|
}
|
|
|
|
|
|
template <int D, typename SCAL, typename SCAL2>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> & operator*= (AutoDiffRec<D,SCAL> & b, SCAL2 a)
|
|
{
|
|
b.Rec() *= a;
|
|
b.Last() *= a;
|
|
return b;
|
|
}
|
|
|
|
/// Inverse of AutoDiffRec
|
|
|
|
template <typename SCAL>
|
|
auto Inv1 (SCAL x) { return 1.0/x; }
|
|
|
|
template<int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> Inv1 (AutoDiffRec<D,SCAL> x)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (Inv1(x.Rec()), (-sqr(1.0/x.Value())) * x.Last());
|
|
}
|
|
|
|
/// AutoDiffRec div AutoDiffRec
|
|
template<int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> operator/ (const AutoDiffRec<D,SCAL> & x, const AutoDiffRec<D,SCAL> & y)
|
|
{
|
|
return x * Inv1 (y);
|
|
}
|
|
|
|
|
|
/// AutoDiffVec div double
|
|
template<int D, typename SCAL, typename SCAL2,
|
|
typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> operator/ (const AutoDiffRec<D,SCAL> & x, SCAL2 y)
|
|
{
|
|
return (1.0/y) * x;
|
|
}
|
|
|
|
|
|
/// double div AutoDiffVec
|
|
template<int D, typename SCAL, typename SCAL2,
|
|
typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> operator/ (SCAL2 x, const AutoDiffRec<D,SCAL> & y)
|
|
{
|
|
return x * Inv1(y);
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
///
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE bool operator== (AutoDiffRec<D,SCAL> x, SCAL val2)
|
|
{
|
|
return x.Value() == val2;
|
|
}
|
|
|
|
///
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE bool operator!= (AutoDiffRec<D,SCAL> x, SCAL val2) throw()
|
|
{
|
|
return x.Value() != val2;
|
|
}
|
|
|
|
///
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE bool operator< (AutoDiffRec<D,SCAL> x, SCAL val2) throw()
|
|
{
|
|
return x.Value() < val2;
|
|
}
|
|
|
|
///
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE bool operator> (AutoDiffRec<D,SCAL> x, SCAL val2) throw()
|
|
{
|
|
return x.Value() > val2;
|
|
}
|
|
|
|
using std::fabs;
|
|
template<int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> fabs (const AutoDiffRec<D,SCAL> & x)
|
|
{
|
|
auto sign = IfPos(x.Value(), SCAL(1.0), IfPos(-x.Value(), SCAL(-1.0), SCAL(0.0)));
|
|
return AutoDiffRec<D,SCAL> (fabs(x.Rec()), sign*x.Last());
|
|
// return fabs (AutoDiffVec<D,SCAL>(x));
|
|
/*
|
|
double abs = fabs (x.Value());
|
|
AutoDiffVec<D,SCAL> res( abs );
|
|
if (abs != 0.0)
|
|
for (int i = 0; i < D; i++)
|
|
res.DValue(i) = x.Value()*x.DValue(i) / abs;
|
|
else
|
|
for (int i = 0; i < D; i++)
|
|
res.DValue(i) = 0.0;
|
|
return res;
|
|
*/
|
|
}
|
|
|
|
|
|
template<int D, typename SCAL>
|
|
NETGEN_INLINE auto sqrt (const AutoDiffRec<D,SCAL> & x)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (sqrt(x.Rec()), (0.5/sqrt(x.Value()))*x.Last());
|
|
}
|
|
|
|
|
|
|
|
template <int D, typename SCAL>
|
|
auto log (AutoDiffRec<D,SCAL> x)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (log(x.Rec()), (1.0/x.Value())*x.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
auto exp (AutoDiffRec<D,SCAL> x)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (exp(x.Rec()), exp(x.Value())*x.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> pow (AutoDiffRec<D,SCAL> x, AutoDiffRec<D,SCAL> y )
|
|
{
|
|
return exp(log(x)*y);
|
|
}
|
|
|
|
|
|
template <int D, typename SCAL>
|
|
auto sin (AutoDiffRec<D,SCAL> x)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (sin(x.Rec()), cos(x.Value())*x.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
auto cos (AutoDiffRec<D,SCAL> x)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (cos(x.Rec()), -sin(x.Value())*x.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
auto tan (AutoDiffRec<D,SCAL> x)
|
|
{
|
|
return sin(x) / cos(x);
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
auto sinh (AutoDiffRec<D,SCAL> x)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (sinh(x.Rec()), cosh(x.Value())*x.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
auto cosh (AutoDiffRec<D,SCAL> x)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (cosh(x.Rec()), sinh(x.Value())*x.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
auto erf (AutoDiffRec<D,SCAL> x)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (erf(x.Rec()), 2. / sqrt(M_PI) * exp(- x.Value() * x.Value())*x.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
auto floor (AutoDiffRec<D,SCAL> x)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (floor(x.Rec()), 0.0);
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
auto ceil (AutoDiffRec<D,SCAL> x)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (ceil(x.Rec()), 0.0);
|
|
}
|
|
|
|
|
|
|
|
template <int D, typename SCAL>
|
|
auto atan (AutoDiffRec<D,SCAL> x)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (atan(x.Rec()), (1./(1.+x.Value()*x.Value()))*x.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
auto atan2 (AutoDiffRec<D,SCAL> x, AutoDiffRec<D,SCAL> y)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (atan2(x.Rec(), y.Rec()),
|
|
(1./(x.Value()*x.Value()+y.Value()*y.Value()))*(x.Value()*y.Last()-y.Value()*x.Last()));
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
auto acos (AutoDiffRec<D,SCAL> x)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (acos(x.Rec()), (-1./sqrt(1.-x.Value()*x.Value()))*x.Last());
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
auto asin (AutoDiffRec<D,SCAL> x)
|
|
{
|
|
return AutoDiffRec<D,SCAL> (asin(x.Rec()), (1./sqrt(1.-x.Value()*x.Value()))*x.Last());
|
|
}
|
|
|
|
|
|
template <int D, typename SCAL, typename TB, typename TC>
|
|
auto IfPos (AutoDiffRec<D,SCAL> a, TB b, TC c) // -> decltype(IfPos (a.Value(), b, c))
|
|
{
|
|
return IfPos (a.Value(), b, c);
|
|
}
|
|
|
|
template <int D, typename SCAL>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffRec<D,SCAL> b, AutoDiffRec<D,SCAL> c)
|
|
{
|
|
/*
|
|
AutoDiffRec<D,SCAL> res;
|
|
res.Value() = IfPos (a, b.Value(), c.Value());
|
|
for (int j = 0; j < D; j++)
|
|
res.DValue(j) = IfPos (a, b.DValue(j), c.DValue(j));
|
|
return res;
|
|
*/
|
|
return AutoDiffRec<D,SCAL> (IfPos(a, b.Rec(), c.Rec()), IfPos(a, b.Last(), c.Last()));
|
|
}
|
|
|
|
template <int D, typename SCAL, typename TC>
|
|
NETGEN_INLINE AutoDiffRec<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffRec<D,SCAL> b, TC c)
|
|
{
|
|
return IfPos (a, b, AutoDiffRec<D,SCAL> (c));
|
|
}
|
|
|
|
|
|
|
|
template <int D, typename SCAL = double>
|
|
using AutoDiff = AutoDiffRec<D,SCAL>;
|
|
|
|
}
|
|
|
|
|
|
|
|
namespace ngbla
|
|
{
|
|
template <typename T> struct is_scalar_type;
|
|
template <int D, typename T>
|
|
struct is_scalar_type<ngcore::AutoDiff<D,T>> { static constexpr bool value = true; };
|
|
|
|
// not meaningful for AutoDiff<D,Complex>, since this is
|
|
// not (complex) differentiable anyway
|
|
template<int D, typename SCAL>
|
|
inline auto L2Norm2 (const ngcore::AutoDiff<D,SCAL> & x)
|
|
{
|
|
return x*x;
|
|
}
|
|
|
|
template<int D, typename SCAL>
|
|
inline auto L2Norm (const ngcore::AutoDiff<D,SCAL> & x) throw()
|
|
{
|
|
return IfPos(x.Value(), x, -x);
|
|
}
|
|
|
|
|
|
|
|
template<int D, typename TAD>
|
|
NETGEN_INLINE auto Conj (const ngcore::AutoDiff<D,TAD> & a)
|
|
{
|
|
ngcore::AutoDiff<D,TAD> b;
|
|
b.Value() = conj(a.Value());
|
|
|
|
for(int i=0;i<D;i++)
|
|
b.DValue(i) = conj(a.DValue(i));
|
|
|
|
return b;
|
|
}
|
|
|
|
|
|
}
|
|
|
|
|
|
#endif
|
|
|
|
|