Using ADIC to generate derivative information

Consider the ODE

(1)

\begin{eqnarray} \dot y_1 =&& - p1_1 y_1^2 - y_3 \\ \dot y_2 =&& - y_2 \\ \dot y_3 =&& -p_2^2 y_2 y_3 \end{eqnarray}

with initial conditions $y_1(0) = y_2(0) = y_3(0) = 1.0$ which depends on the two parameters $p_1 = 1.0$ and $p_2 = 2.0$.

int myrhs(double t, double *y, double *ydot, double *p)
{
  ydot[0] = -p[0]*y[0]*y[0] - y[2];
  ydot[1] = -y[1];
  ydot[2] = -p[1]*p[1]*y[1]*y[2];
 
  return(0);
}

#if !defined(AD_DERIV_H)
#define AD_DERIV_H
 
typedef double InactiveDouble;
typedef float InactiveFloat;
 
#if !defined(ad_GRAD_PTR) 
#define ad_GRAD_PTR 0
#endif
 
#if !defined(ad_GRAD_MAX) 
#define ad_GRAD_MAX 5
#endif
 
#define AD_INIT_MAP()
#define AD_CLEANUP_MAP()
#define AD_GET_DERIV_OBJ(x) ((void*)(&x.value+1))
#define AD_FREE_DERIV_OBJ(x)
typedef struct {
    double value;
    double  grad[ad_GRAD_MAX];
} DERIV_TYPE;
 
#define DERIV_val(a) ((a).value)
 
#define DERIV_grad(a) ((a).grad)
 
#undef _FLOAT_INITIALIZER_
 
#define _FLOAT_INITIALIZER_(x) { x, 0.0 }
 
void ad_AD_Init(int);
void ad_AD_Final();
#include "ad_grad.h"
 
#define nullFunc(x) 0
 
#endif

#include "ad_deriv.h"
 
int ad_myrhs(DERIV_TYPE t, DERIV_TYPE *y, DERIV_TYPE *ydot, DERIV_TYPE *p) {
  double  ad_loc_0;
  double  ad_loc_1;
  double  ad_loc_2;
  double  ad_loc_3;
  double  ad_adj_0;
  double  ad_adj_1;
  double  ad_adj_2;
  double  ad_adj_3;
  double  ad_adj_4;
  {
    ad_loc_0 =  -DERIV_val(p[0]);
    ad_loc_1 = ad_loc_0 * DERIV_val(y[0]);
    ad_loc_2 = ad_loc_1 * DERIV_val(y[0]);
    ad_loc_3 = ad_loc_2 - DERIV_val(y[2]);
    ad_adj_0 = ad_loc_0 * DERIV_val(y[0]);
    ad_adj_1 = DERIV_val(y[0]) * DERIV_val(y[0]);
    ad_adj_2 =  -ad_adj_1;
    ad_grad_axpy_4(&(ydot[0]), ad_adj_2, &(p[0]), ad_adj_0, &(y[0]), ad_loc_1, &(y[0]), -1.000000000000000e+00, &(y[2]));
    DERIV_val(ydot[0]) = ad_loc_3;
  }
  {
    ad_loc_0 =  -DERIV_val(y[1]);
    ad_grad_axpy_1(&(ydot[1]), -1.000000000000000e+00, &(y[1]));
    DERIV_val(ydot[1]) = ad_loc_0;
  }
  {
    ad_loc_0 =  -DERIV_val(p[1]);
    ad_loc_1 = ad_loc_0 * DERIV_val(p[1]);
    ad_loc_2 = ad_loc_1 * DERIV_val(y[1]);
    ad_loc_3 = ad_loc_2 * DERIV_val(y[2]);
    ad_adj_0 = ad_loc_1 * DERIV_val(y[2]);
    ad_adj_1 = DERIV_val(y[1]) * DERIV_val(y[2]);
    ad_adj_2 = ad_loc_0 * ad_adj_1;
    ad_adj_3 = DERIV_val(p[1]) * ad_adj_1;
    ad_adj_4 =  -ad_adj_3;
    ad_grad_axpy_4(&(ydot[2]), ad_adj_4, &(p[1]), ad_adj_2, &(p[1]), ad_adj_0, &(y[1]), ad_loc_2, &(y[2]));
    DERIV_val(ydot[2]) = ad_loc_3;
  }
  return 0;
}
 
void ad_AD_Init(int arg0) {
  ad_AD_GradInit(arg0); 
}
 
void ad_AD_Final() {
  ad_AD_GradFinal();
}