libtrilinos-sacado-dev 12.12.1-5 / usr / include / trilinos / Sacado_rad2.hpp

// @HEADER
// ***********************************************************************
//
//                           Sacado Package
//                 Copyright (2007) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
// USA
// Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
// (etphipp@sandia.gov).
//
// ***********************************************************************
// @HEADER

// Extension of the RAD package (Reverse Automatic Differentiation) --
// a package specialized for function and gradient evaluations -- to
// Hessian-vector products.
// This variant is for Hessian-vector products of "double" variables, and
// Sacado::Rad2d::ADvar should be equivalent to Sacado::Rad2::ADvar<double>,
// but this nontemplated code may easier to understand.  It relies on ops
// implemented in Sacado_radops2.cpp.
// Written in 2007 by David M. Gay at Sandia National Labs, Albuquerque, NM.

#ifndef SACADO_RAD2_H
#define SACADO_RAD2_H

#include <stddef.h>
#include <Sacado_cmath.hpp>

#include "Sacado_ConfigDefs.h"

#if defined(RAD_DEBUG_BLOCKKEEP) && !defined(HAVE_SACADO_UNINIT)
#undef RAD_DEBUG_BLOCKKEEP
#endif

#ifndef SACADO_NO_NAMESPACE
namespace Sacado {
namespace Rad2d { // "2" for 2nd derivatives, "d" for "double"
#endif

// -DNO_USING_STDCC is needed, e.g., with Sun CC 5.7
#ifndef RAD_NO_USING_STDCC
  // Bring math functions into scope
  using std::exp;
  using std::log;
  using std::log10;
  using std::sqrt;
  using std::cos;
  using std::sin;
  using std::tan;
  using std::acos;
  using std::asin;
  using std::atan;
  using std::cosh;
  using std::sinh;
  using std::tanh;
  using std::abs;
  using std::fabs;
  using std::atan2;
  using std::pow;
#endif //NO_USING_STDCC

 class ADvar;
 class ADvari;
 class ADvar1;
 class ADvar2;
 class ConstADvar;
 class Derp;
 class IndepADvar;

 struct
ADmemblock {	// We get memory in ADmemblock chunks and never give it back,
		// but reuse it once computations start anew after call(s) on
		// ADcontext::Gradcomp() or ADcontext::Weighted_Gradcomp().
	ADmemblock *next;
	double memblk[2000];
	};

 struct
ADvari_block {
	enum { Gulp = 1021 };
	ADvari_block *next, *prev;
	ADvari **limit;
	ADvari *pADvari[Gulp];
	};

 class
ADcontext {	// A singleton class: one instance in radops.c
	ADmemblock *Busy, *Free;
	char *Mbase;
	size_t Mleft;
	ADvari **Ailimit, **Ainext;
	ADvari_block *Aibusy, *Aifree;
	ADmemblock First;
	ADvari_block AiFirst;
	void *new_ADmemblock(size_t);
	void new_ADvari_block();
 public:
	ADcontext();
	void *Memalloc(size_t len);
	static void Gradcomp(int);
	static inline void Gradcomp() { Gradcomp(1); }
	static void Hvprod(int, ADvar**, double*, double*);
	static void Weighted_Gradcomp(int, ADvar**, double*);
	inline void ADvari_record(ADvari *x) {
		if (Ainext >= Ailimit)
			new_ADvari_block();
		*Ainext++ = x;
		}
	};

inline void *ADcontext::Memalloc(size_t len) {
		if (Mleft >= len)
			return Mbase + (Mleft -= len);
		return new_ADmemblock(len);
		}

 class
CADcontext: public ADcontext {
 protected:
	bool fpval_implies_const;
 public:
	friend class ADvar;
	CADcontext(): ADcontext() { fpval_implies_const = false; }
	static const double One, negOne;
	};

 class
Derp {		// one derivative-propagation operation
 public:
	friend class ADvarn;
	static Derp *LastDerp;
	Derp *next;
	const double *a;
	const ADvari *b;
	mutable ADvari *c;
	void *operator new(size_t);
	void operator delete(void*) {} /*Should never be called.*/
	inline Derp(){};
	Derp(const ADvari *);
	Derp(const double *, const ADvari *);
	Derp(const double *, const ADvari *, const ADvari *);
	/* c->aval += a * b->aval; */
	};

inline Derp::Derp(const ADvari *c1): c((ADvari*)c1) {
		next = LastDerp;
		LastDerp = this;
		}

inline Derp::Derp(const double *a1, const ADvari *c1):
		a(a1), c((ADvari*)c1) {
		next = LastDerp;
		LastDerp = this;
		}

inline Derp::Derp(const double *a1, const ADvari *b1, const ADvari *c1):
		a(a1), b(b1), c((ADvari*)c1) {
		next = LastDerp;
		LastDerp = this;
		}

 enum Advari_Opclass {
	Hv_const,
	Hv_copy,
	Hv_binary,
	Hv_unary,
	Hv_negate,
	Hv_plusLR,
	Hv_minusLR,
	Hv_timesL,
	Hv_timesLR,
	Hv_quotLR,
	Hv_nary
	};

 extern ADvari& ADf1(double f, double g, double h, const ADvari &x);
 extern ADvari& ADf2(double f, double gx, double gy, double hxx,
			double hxy, double hyy, const ADvari &x, const ADvari &y);
 extern ADvari& ADfn(double f, int n, const ADvar *x, const double *g, const double *h);

 class
ADvari {	// implementation of an ADvar
 public:
	static ADcontext adc;
	Advari_Opclass opclass;
	double Val;		// result of this operation
	mutable double aval;	// adjoint -- partial of final result w.r.t. this Val
	mutable double dO;	// deriv of op w.r.t. t in x + t*p
	mutable double aO;	// adjoint (in Hv computation) of op
	mutable double adO;	// adjoint (in Hv computation) of dO
	void *operator new(size_t len) { return ADvari::adc.Memalloc(len); }
	void operator delete(void*) {} /*Should never be called.*/
	inline ADvari(Advari_Opclass oc, double t):
		opclass(oc), Val(t), aval(0.), dO(0.)
		{ if (oc != Hv_const) ADvari::adc.ADvari_record(this); }
	inline ADvari(Advari_Opclass oc, double t, double ta):
		opclass(oc), Val(t), aval(ta), dO(0.)
		{ if (oc != Hv_const) ADvari::adc.ADvari_record(this); }
 private:
	inline ADvari(): Val(0.), aval(0.), dO(0.) {}	// prevent construction without value (?)
 public:
	friend class ConstADvari;
#ifdef RAD_AUTO_AD_Const
	friend class ADcontext;
	friend class ADvar1;
	friend class ADvar;
	friend class ConstADvar;
	friend class IndepADvar;
 private:
	ADvari *Next;
	static ADvari *First_ADvari, **Last_ADvari;
 protected:
	IndepADvar *padv;
 public:
	ADvari(const IndepADvar *, double);
#endif
#define F friend
#define R ADvari&
#define Ai const ADvari&
#define T1(r,f) F r f(Ai);
#define T2(r,f) \
F r f(Ai,Ai); \
F r f(double,Ai); \
F r f(Ai,double);
	T1(R,operator+)
	T2(R,operator+)
	T1(R,operator-)
	T2(R,operator-)
	T2(R,operator*)
	T2(R,operator/)
	T1(R,abs)
	T1(R,acos)
	T1(R,acosh)
	T1(R,asin)
	T1(R,asinh)
	T1(R,atan)
	T1(R,atanh)
	T2(R,atan2)
	T2(R,max)
	T2(R,min)
	T1(R,cos)
	T1(R,cosh)
	T1(R,exp)
	T1(R,log)
	T1(R,log10)
	T2(R,pow)
	T1(R,sin)
	T1(R,sinh)
	T1(R,sqrt)
	T1(R,tan)
	T1(R,tanh)
	T1(R,fabs)
	T1(R,copy)
	T2(int,operator<)
	T2(int,operator<=)
	T2(int,operator==)
	T2(int,operator!=)
	T2(int,operator>=)
	T2(int,operator>)
#undef T2
#undef T1
#undef Ai
#undef R
#undef F

	friend ADvari& ADf1(double f, double g, double h, const ADvari &x);
	friend ADvari& ADf2(double f, double gx, double gy, double hxx,
			double hxy, double hyy, const ADvari &x, const ADvari &y);
	friend ADvari& ADfn(double f, int n, const ADvar *x, const double *g, const double *h);
	};

 inline void* Derp::operator new(size_t len) { return ADvari::adc.Memalloc(len); }


 class
ADvar1: public ADvari {	// simplest unary ops
 public:
	Derp d;
	ADvar1(Advari_Opclass oc, double val1, const double *a1, const ADvari *c1):
		ADvari(oc,val1), d(a1,this,c1) {}
#ifdef RAD_AUTO_AD_Const
	ADvar1(const IndepADvar *, const IndepADvar &);
	ADvar1(const IndepADvar *, const ADvari &);
#endif
	};

 class
ConstADvari: public ADvari {
 private:
	ConstADvari *prevcad;
	ConstADvari() {};	// prevent construction without value (?)
	static ConstADvari *lastcad;
 public:
	static CADcontext cadc;
	inline void *operator new(size_t len) { return ConstADvari::cadc.Memalloc(len); }
	inline ConstADvari(double t): ADvari(Hv_copy, t) { prevcad = lastcad; lastcad = this; }
	static void aval_reset(void);
	};

 class
IndepADvar
{
 private:
	inline IndepADvar& operator=(const IndepADvar &x) {
		/* private to prevent assignment */
#ifdef RAD_AUTO_AD_Const
		if (cv)
			cv->padv = 0;
		cv = new ADvar1(this,x);
		return *this;
#else
#ifdef RAD_EQ_ALIAS
		cv = x.cv;
		return *this;
#else
		return ADvar_operatoreq(this,*x.cv);
#endif
#endif /* RAD_AUTO_AD_Const */
		}
 protected:
	static void AD_Const(const IndepADvar&);
	ADvari *cv;
 public:
	typedef double value_type;
	friend class ADvar;
	friend class ADvar1;
	friend class ADvarn;
	friend class ADcontext;
	IndepADvar(double);
	IndepADvar(int);
	IndepADvar(long);
	IndepADvar& operator=(double);
#ifdef RAD_AUTO_AD_Const
	inline IndepADvar(const IndepADvar &x) { cv = x.cv ? new ADvar1(this, x) : 0; };
	inline IndepADvar() { cv = 0; }
	inline ~IndepADvar() {
			if (cv)
				cv->padv = 0;
		}
#else
	inline IndepADvar() {
#ifndef RAD_EQ_ALIAS
		cv = 0;
#endif
		}
	inline ~IndepADvar() {}
	friend IndepADvar& ADvar_operatoreq(IndepADvar*, const ADvari&);
#endif

	friend void AD_Const(const IndepADvar&);

	inline operator ADvari&() { return *cv; }
	inline operator ADvari&() const { return *(ADvari*)cv; }

	inline double val() const { return cv->Val; }
	inline double adj() const { return cv->aval; }
	static inline void Gradcomp(int wantgrad)
				{ ADcontext::Gradcomp(wantgrad); }
	static inline void Gradcomp()
				{ ADcontext::Gradcomp(1); }
	static inline void Hvprod(int n, ADvar **vp, double *v, double *hv)
				{ ADcontext::Hvprod(n, vp, v, hv); }
	static inline void aval_reset() { ConstADvari::aval_reset(); }
	static inline void Weighted_Gradcomp(int n, ADvar **v, double *w)
				{ ADcontext::Weighted_Gradcomp(n, v, w); }


#define Ai const ADvari&
#define AI const IndepADvar&
#define T2(r,f) \
 r f(AI,AI);\
 r f(Ai,AI);\
 r f(AI,Ai);\
 r f(double,AI);\
 r f(AI,double);
#define T1(f) friend ADvari& f(AI);

#define F friend ADvari&
T2(F, operator+)
T2(F, operator-)
T2(F, operator*)
T2(F, operator/)
T2(F, atan2)
T2(F, max)
T2(F, min)
T2(F, pow)
#undef F
#define F friend int
T2(F, operator<)
T2(F, operator<=)
T2(F, operator==)
T2(F, operator!=)
T2(F, operator>=)
T2(F, operator>)

T1(operator+)
T1(operator-)
T1(abs)
T1(acos)
T1(acosh)
T1(asin)
T1(asinh)
T1(atan)
T1(atanh)
T1(cos)
T1(cosh)
T1(exp)
T1(log)
T1(log10)
T1(sin)
T1(sinh)
T1(sqrt)
T1(tan)
T1(tanh)
T1(fabs)
T1(copy)

#undef T1
#undef T2
#undef F
#undef Ai
#undef AI

	};

 class
ADvar: public IndepADvar {		// an "active" variable
	void ADvar_ctr(double d);
 public:
	inline ADvar() { /* cv = 0; */ }
	inline ADvar(double d) { ADvar_ctr(d); }
	inline ADvar(int i)	{ ADvar_ctr((double)i); }
	inline ADvar(long i)	{ ADvar_ctr((double)i); }
	inline ~ADvar() {}
#ifdef RAD_AUTO_AD_Const
	friend class ADvar1;
	inline ADvar(const IndepADvar &x) { cv = x.cv ? new ADvar1(this, x) : 0; }
	inline ADvar(ADvari &x) { cv = &x; x.padv = this; }
	inline ADvar& operator=(const IndepADvar &x) {
		if (cv)
			cv->padv = 0;
		cv = new ADvar1(this,x);
		return *this;
		}
	inline ADvar& operator=(const ADvari &x) {
		if (cv)
			cv->padv = 0;
		cv = new ADvar1(this, x);
		return *this;
		}
#else
	friend ADvar& ADvar_operatoreq(ADvar*, const ADvari&);
#ifdef RAD_EQ_ALIAS
	/* allow aliasing v and w after "v = w;" */
	inline ADvar(const IndepADvar &x) { cv = x.cv; }
	inline ADvar(const ADvari &x) { cv = (ADvari*)&x; }
	inline ADvar& operator=(const ADvari &x) { cv = (ADvari*)&x;   return *this; }
	inline ADvar& operator=(const IndepADvar &x)  { cv = (ADvari*)x.cv; return *this; }
#else
	ADvar(const IndepADvar &x) { cv = x.cv ?
		new ADvar1(Hv_copy, x.cv->Val, &CADcontext::One, x.cv) : 0; }
	ADvar(const ADvari &x) { cv = new ADvar1(Hv_copy, x.Val, &CADcontext::One, &x); }
	inline ADvar& operator=(const ADvari &x) { return ADvar_operatoreq(this,x); }
	inline ADvar& operator=(const IndepADvar &x)    { return ADvar_operatoreq(this,*x.cv); }
#endif /* RAD_EQ_ALIAS */
#endif /* RAD_AUTO_AD_Const */
	ADvar& operator=(double);
	ADvar& operator+=(const ADvari&);
	ADvar& operator+=(double);
	ADvar& operator-=(const ADvari&);
	ADvar& operator-=(double);
	ADvar& operator*=(const ADvari&);
	ADvar& operator*=(double);
	ADvar& operator/=(const ADvari&);
	ADvar& operator/=(double);
	inline static bool get_fpval_implies_const(void)
		{ return ConstADvari::cadc.fpval_implies_const; }
	inline static void set_fpval_implies_const(bool newval)
		{ ConstADvari::cadc.fpval_implies_const = newval; }
	inline static bool setget_fpval_implies_const(bool newval) {
		bool oldval = ConstADvari::cadc.fpval_implies_const;
		ConstADvari::cadc.fpval_implies_const = newval;
		return oldval;
		}
	static inline void Gradcomp(int wantgrad)
				{ ADcontext::Gradcomp(wantgrad); }
	static inline void Gradcomp()
				{ ADcontext::Gradcomp(1); }
	static inline void Hvprod(int n, ADvar **vp, double *v, double *hv)
				{ ADcontext::Hvprod(n, vp, v, hv); }
	static inline void aval_reset() { ConstADvari::aval_reset(); }
	static inline void Weighted_Gradcomp(int n, ADvar **v, double *w)
				{ ADcontext::Weighted_Gradcomp(n, v, w); }
	};

 inline void AD_Const(const IndepADvar&v) { IndepADvar::AD_Const(v); }

 class
ConstADvar: public ADvar {
 private: // disable op=
	ConstADvar& operator+=(const ADvari&);
	ConstADvar& operator+=(double);
	ConstADvar& operator-=(const ADvari&);
	ConstADvar& operator-=(double);
	ConstADvar& operator*=(const ADvari&);
	ConstADvar& operator*=(double);
	ConstADvar& operator/=(const ADvari&);
	ConstADvar& operator/=(double);
	void ConstADvar_ctr(double);
 public:
	inline ConstADvar(double d)	{ ConstADvar_ctr(d); }
	inline ConstADvar(int i)	{ ConstADvar_ctr((double)i); }
	inline ConstADvar(long i)	{ ConstADvar_ctr((double)i); }
	ConstADvar(const ADvari &x);
#ifdef RAD_AUTO_AD_Const
	ConstADvar(const IndepADvar &x) { cv = new ADvar1(this,x); }
#endif
	inline ~ConstADvar(){}
#ifdef RAD_NO_CONST_UPDATE
 private:
#endif
	ConstADvar();
	inline ConstADvar& operator=(double d) { cv->Val = d; return *this; }
	inline ConstADvar& operator=(const IndepADvar& d) { cv->Val = d.val(); return *this; }
 };

 class
ADvar1s: public ADvar1 { // unary ops with partials
 public:
	double pL;	// deriv of op w.r.t. left operand L
	ADvar1s(double val1, double d1, const ADvari *c1):
		ADvar1(Hv_timesL,val1,&pL,c1), pL(d1) {}
	};

 class
ADvar1g: public ADvar1 { // unary ops with partials
 public:
	double pL;	// deriv of op w.r.t. left operand L
	double pL2;	// partial of op w.r.t. L,L
	ADvar1g(double val1, double d1, double d2, const ADvari *c1):
		ADvar1(Hv_unary,val1,&pL,c1), pL(d1), pL2(d2) {}
	};

 class
ADvar2: public ADvari {	// basic binary ops
 public:
	Derp dL, dR;
	ADvar2(Advari_Opclass oc, double val1, const ADvari *Lcv, const double *Lc,
			const ADvari *Rcv, const double *Rc): ADvari(oc,val1) {
		dR.next = Derp::LastDerp;
		dL.next = &dR;
		Derp::LastDerp = &dL;
		dL.a = Lc;
		dL.c = (ADvari*)Lcv;
		dR.a = Rc;
		dR.c = (ADvari*)Rcv;
		dL.b = dR.b = this;
		}
	};

 class
ADvar2q: public ADvar2 { // binary ops with partials
 public:
	double pL;	// deriv of op w.r.t. left operand L
	double pR;	// deriv of op w.r.t. right operand R
	double pLR;	// second partial w.r.t. L,R
	double pR2;	// second partial w.r.t. R,R
	ADvar2q(double val1, double Lp, double Rp, double LR, double R2,
		const ADvari *Lcv, const ADvari *Rcv);
	};

 class
ADvar2g: public ADvar2 { // general binary ops with partials
 public:
	double pL;	// deriv of op w.r.t. left operand L
	double pR;	// deriv of op w.r.t. right operand R
	double pL2;	// second partial w.r.t. L,L
	double pLR;	// second partial w.r.t. L,R
	double pR2;	// second partial w.r.t. R,R
	ADvar2g(double val1, double Lp, double Rp, double L2, double LR, double R2,
		const ADvari *Lcv, const ADvari *Rcv);
	};

 class
ADvarn: public ADvari { // n-ary ops with partials g and 2nd partials h (lower triangle, rowwise)
 public:
	int n;
	double *G, *H;
	Derp *D;
	ADvarn(double val1, int n1, const ADvar *x, const double *g, const double *h);
	};

inline ADvari &operator+(ADvari &T) { return T; }
inline ADvari &operator+(const ADvari &T) { return (ADvari&) T; }

inline int operator<(const ADvari &L, const ADvari &R) { return L.Val < R.Val; }
inline int operator<(const ADvari &L, double R) { return L.Val < R; }
inline int operator<(double L, const ADvari &R) { return L < R.Val; }

inline int operator<=(const ADvari &L, const ADvari &R) { return L.Val <= R.Val; }
inline int operator<=(const ADvari &L, double R) { return L.Val <= R; }
inline int operator<=(double L, const ADvari &R) { return L <= R.Val; }

inline int operator==(const ADvari &L, const ADvari &R) { return L.Val == R.Val; }
inline int operator==(const ADvari &L, double R) { return L.Val == R; }
inline int operator==(double L, const ADvari &R) { return L == R.Val; }

inline int operator!=(const ADvari &L, const ADvari &R) { return L.Val != R.Val; }
inline int operator!=(const ADvari &L, double R) { return L.Val != R; }
inline int operator!=(double L, const ADvari &R) { return L != R.Val; }

inline int operator>=(const ADvari &L, const ADvari &R) { return L.Val >= R.Val; }
inline int operator>=(const ADvari &L, double R) { return L.Val >= R; }
inline int operator>=(double L, const ADvari &R) { return L >= R.Val; }

inline int operator>(const ADvari &L, const ADvari &R) { return L.Val > R.Val; }
inline int operator>(const ADvari &L, double R) { return L.Val > R; }
inline int operator>(double L, const ADvari &R) { return L > R.Val; }

inline ADvari& copy(const IndepADvar &x)
{ return *(new ADvar1(Hv_copy, x.cv->Val, &CADcontext::One, x.cv)); }

inline ADvari& copy(const ADvari &x)
{ return *(new ADvar1(Hv_copy, x.Val, &CADcontext::One, &x)); }

inline ADvari& abs(const ADvari &x)
{ return fabs(x); }

#define A (ADvari*)
#define T inline
#define F ADvari&
#define Ai const ADvari&
#define AI const IndepADvar&
#define D double
#define T2(r,f) \
 T r f(Ai L, AI R) { return f(L, *A R.cv); }\
 T r f(AI L, Ai R) { return f(*A L.cv, R); }\
 T r f(AI L, AI R) { return f(*A L.cv, *A R.cv); }\
 T r f(AI L, D R) { return f(*A L.cv, R); }\
 T r f(D L, AI R) { return f(L, *A R.cv); }

T2(F, operator+)
T2(F, operator-)
T2(F, operator*)
T2(F, operator/)
T2(F, atan2)
T2(F, pow)
T2(F, max)
T2(F, min)
T2(int, operator<)
T2(int, operator<=)
T2(int, operator==)
T2(int, operator!=)
T2(int, operator>=)
T2(int, operator>)

#undef T2
#undef D

#define T1(f)\
 T F f(AI x) { return f(*A x.cv); }

T1(operator+)
T1(operator-)
T1(abs)
T1(acos)
T1(acosh)
T1(asin)
T1(asinh)
T1(atan)
T1(atanh)
T1(cos)
T1(cosh)
T1(exp)
T1(log)
T1(log10)
T1(sin)
T1(sinh)
T1(sqrt)
T1(tan)
T1(tanh)
T1(fabs)

#undef T1
#undef AI
#undef Ai
#undef F
#undef T
#undef A

#ifndef SACADO_NO_NAMESPACE
} // namespace Rad2d
} // namespace Sacado
#endif // SACADO_NAMESPACE
#endif // SACADO_RAD2_H