[ADOL-C] gcc or g++

[Kharche, Sanjay]

Hi

Sorry for my previous posting with the error - I can see that the Jacobian is indeed computed. thanks.


________________________________________
From: adol-c-bounces at list.coin-or.org [adol-c-bounces at list.coin-or.org] on behalf of Kharche, Sanjay
Sent: 13 December 2013 18:07
To: Kshitij Kulshreshtha; adol-c at list.coin-or.org
Subject: Re: [ADOL-C] gcc or g++

Hi

Many thanks for the rapid reply. I worked through the odexam.cpp example. accode sets the values of 2 3D arrays (code pasted below) and 1 2D array, none of which are the Jacobian of the ODE that is defined in tracerhs. How do I get the Jacobian?

Also, if you have a simple example of how to calculate a Jacobian from ADOL-C and then feed that into CVODES, that will help in the local sensitivity analysis that I would like to do.

thanks
Sanjay

Code:

//
#include "adolc/adolc.h"
//

#include <iostream>
#include<math.h>
#include<stdio.h>

/****************************************************************************/
/*                                                          ADOUBLE ROUTINE */
void tracerhs(short int tag, double* py, double* pyprime) {
    adouble y[3];                // we left the parameters passive
    adouble yprime[3];

    trace_on(tag);
    y[0] <<= py[0]; y[1] <<= py[1];   y[2] <<= py[2];                              // initialize and mark independents
    yprime[0] = -sin(y[2]) + 1.0e8*y[2]*(1.0-1.0/y[0]);
    yprime[1] = -10.0*y[0] + 3.0e7*y[2]*(1-y[1]);
    yprime[2] = -yprime[0] - yprime[1];
    yprime[0] >>= pyprime[0]; yprime[1] >>= pyprime[1];  yprime[2] >>= pyprime[2]; // mark and pass dependents
    trace_off(1);                // write tape array onto a file
} // end tracerhs

/****************************************************************************/
/*                                                             MAIN PROGRAM */
int main() {
    int i,j,deg,k;
    int n = 3;
    double py[3];
    double pyp[3];

    deg = 5; // this is an input.

    short** nz = new short*[n];
    double  **X;
    double ***Z; // what are Z and B?
    double ***B;

    X = myalloc2(n,deg+1); // these are ADOL-C functions
    Z = myalloc3(n,n,deg);
    B = myalloc3(n,n,deg);

    for(i=0; i<n; i++) {
        py[i]   = (i==0) ? 1.0 : 0.0;       // Initialize the base point
        X[i][0] = py[i];                    // and the Taylor coefficient;
        nz[i]   = new short[n];             // set up sparsity array. So this array is sparsity,
    } // end for

    tracerhs(1,py,pyp);                   // trace RHS with tag = 1
// what does forode do?
    forode(1,n,deg,X);                    // compute deg coefficients
// what does reverse do?
    reverse(1,n,n,deg-1,Z,nz);            // U defaults to the identity
    accode(n,deg-1,Z,B,nz);

    printf("nonzero pattern:\n");
    for(i=0; i<n; i++) {
        for(j=0; j<n; j++)
            printf("%d \t",nz[i][j]);
        printf("\n");
    } // end for

// see what is in B and Z for all degrees of the Taylor polynomial.

for(k = 0; k<deg; k++){

    printf("B (at k = %d):\n",k);
    for(i=0; i<n; i++) {
        for(j=0; j<n; j++)
            printf("%f \t",B[i][j][k]);
        printf("\n");
    } // end for
printf("\n\n");

    printf("Z (at k = %d):\n",k);
    for(i=0; i<n; i++) {
        for(j=0; j<n; j++)
            printf("%f \t",Z[i][j][k]);
        printf("\n");
    } // end for

printf("\n\n");
}

return 1;
} // end main

_______________________________________________
ADOL-C mailing list
ADOL-C at list.coin-or.org
http://list.coin-or.org/mailman/listinfo/adol-c