[ADOL-C] tensor_eval problem for more than one direction

[Tilman Barz]

Dear all,
I am using Adolc in Ubuntu from Eclipse and have a problem using 
<tensor_eval>.
If I compute tensors for more than one direction (simultaneously by one 
function call), in some cases I get a zero column in the tensor output 
array. I encountered this problem for tensors of first and second order.

I changed the "Test driver 'tensor_eval(..)' " -> the problem occurs in 
line 134,144.
columns 11-19 are missing/ are zero/ not computed?

In the output displays:
"############# HERE IS THE PROBLEM ##################"

I'm not sure if this is a bug or if I am using ADOLC in the wrong way.
Any help would be appreciated,
Tilman



/*----------------------------------------------------------------------------
      ADOL-C -- Automatic Differentiation by Overloading in C++
      File:     taylorexam.cpp
      Revision: $Id$
      Contents: Test driver 'tensor_eval(..)' to compute
                higher order derivatives

      Copyright (c) Andrea Walther, Andreas Griewank

      This file is part of ADOL-C. This software is provided as open source.
      Any use, reproduction, or distribution of the software constitutes
      recipient's acceptance of the terms of the accompanying license file.

---------------------------------------------------------------------------*/

/****************************************************************************/
/* INCLUDES */
     #include <adolc/adolc.h>
     #include <cstdlib>
     #include <iostream>


     using namespace std;



     void printMatrix( int nRow, int nCol, double **Mat ){
         cout << "Mat[" << nRow << "][" << nCol << "] : " << endl;
         // --- write all elements in a row
         // write column numbers
         printf("        |");
         for (int j = 0; j < nCol; ++j) {
             printf("%10i|", j);
         }
         printf("\n");
         // write elements
         for (int i = 0; i < nRow; i++) {
             printf("-%4i:  |", i); // write row numbers
             for (int j = 0; j < nCol; j++) {
                 if ( (Mat[i][j] > 1.E-15) or (Mat[i][j] < -1.E-15) )  {
                     printf("%+9.3E|", Mat[i][j]);
                 } else {
                     printf(" -------- |");
                 }
             }
             printf("\n");
         }
         // write column numbers
         printf("        |");
         for (int j = 0; j < nCol; ++j) {
             printf("%10i|", j);
         }
         printf("\n");
      }


/****************************************************************************/
/* MAIN */
     int main() {
         int i,j,m,n,d,p,dim;

/*--------------------------------------------------------------------------*/
         cout << "TAYLOREXAM (ADOL-C Example)\n\n";                      
/* inputs */
         cout << " Number of independents = 17\n";
         n=17;
         cout << " Number of dependents = (<=n) 2\n";
         m=2;
         cout << " Degree = 1\n";
         d=1;
         cout << " Number of directions = 20\n\n";
         p=20;

/*--------------------------------------------------------------------------*/
         int* multi = new int[d];                         /* allocations 
and inits */
         double* xp = new double[n];
         double* yp = new double[m];
         double** S = new double*[n];
         double* test = new double[m];
         double** tensoren;
         double** allTensoren;
         adouble* x = new adouble[n];
         adouble* y = new adouble[m];


         // alternative Seed
         for (i=0; i<n; i++) {
             xp[i] = (i+1.0)/(2.0+i);
             S[i] = new double[p];
             for (j=0; j<p; j++)
                 S[i][j] = 1.0;
         }


/*--------------------------------------------------------------------------*/
         trace_on(1);                                       /* tracing 
the function */
         // adouble* x = new adouble[n];
         // adouble* y = new adouble[m];
         y[0] = 1;

         for (i=0; i<n; i++) {
             x[i] <<= xp[i];
             y[0] *= x[i];
         }
         for (i=1; i<m; i++)
             y[i] = x[i]  * pow(y[0],3);
         for (i=0; i<m; i++)
             y[i] >>= yp[i] ;
         trace_off();

/*--------------------------------------------------------------------------*/
         dim = binomi(p+d,d);
         cout <<"TASK 1:\n";
         cout <<" degree-d = "<<d<<", dim = "<<dim<<", numDir-p = " << p 
<<", dep-m = " << m << ", indp-n = " << n << "\n";
         cout <<" dim = binomi(p+d,d) = " << dim << "\n";
         tensoren = myalloc2(m,dim);

         double* einTensor = new double[m];

         allTensoren = myalloc2(m,dim);
         for (int io = 0; io < m; ++io) {
             for (int jk = 0; jk < dim; ++jk) {
                 allTensoren[io][jk] = 0.0;
             }
         }
//        // initialize with all
         cout << "############# HERE IS THE PROBLEM ##################" 
<< endl;
         tensor_eval(1,m,n,d,19,xp,tensoren,S);
         cout << "tensoren" << endl;
         printMatrix( m, dim, tensoren );


         for (int io = 7; io <= p; ++io) {
             cout << "\n>>>>>> tensor_eval with " << io << " 
directions." << endl;

             if (io == 19) {
                 cout << "############# HERE IS THE PROBLEM 
##################" << endl;
             }

//            for (int jk = 0; jk < 1; ++jk) {
//                cout << "repeat: " << jk << endl;
                 tensor_eval(1,m,n,d,io,xp,tensoren,S);
//            }
             cout << "tensoren" << endl;
             printMatrix( m, dim, tensoren );

             multi[0] = 0;

             for (int pos = 1; pos <= io; ++pos) { // num directions
                 multi[0] = pos;
                 tensor_value(d,m,einTensor,tensoren,multi);

                 // save the direct Deriv
                 for (int hj = 0; hj < m; ++hj) { // number of dependent 
vars
                     allTensoren[hj][pos] = einTensor[hj];
                 }

             }

             cout << "allTensoren" << endl;
             printMatrix( m, dim, allTensoren );
         }


         myfree2(tensoren);
         myfree(einTensor);
         myfree2(allTensoren);


         // ++++++++++++++++++++++++++++++++++++++++++++++++++++++
         return 1;
     }


/****************************************************************************/
     /* THAT'S ALL */