[ADOL-C] sparse hessian & parameters & operator +=

Kshitij Kulshreshtha kshitij at math.upb.de
Mon Jun 27 09:50:42 EDT 2016


Hallo Andrea,

the reason seems to be in the handling of nonlinear index domain in
eq_plus_prod operation at uni5_for.c line 2315. This somehow always
results in diagonal hessian elements being part of the pattern. This
does not happen if one uses the old implementation with
extend_nonlinearity_domain_binary()  (option == 2 or 3 in hess_pat).
Writing x = x + y*z instead of x += y*z causes trace to contain
(mult_a_a, plus_a_a, assign_a) instead of eq_plus_prod, so the problem
does not occur. Since you wrote this new implementation using fod[] and
traverse_unary(), I'm not quite familiar with the logic.

Viele Grüße
Kshitij.


On 27.06.2016 12:12, Andrea Walther wrote:
> Hi,
> 
> thanks very much for the hint!!
> 
> We will have a look at this as soon as possible.
> 
> Best regards
> 
> Andrea Walther
> 
> Am 23.06.2016 um 19:22 schrieb João Leal:
>> Hello,
>>
>> I have a small program using Adol-c (adapted from the examples) to
>> determine a sparse Hessian of a linear function using parameters with
>> the new mkparam function.
>>
>> I have notice that if I use the operator += in my model then my program
>> will produce an incorrect sparsity pattern for the hessian.
>> If that operation is replaced with "x = x + ..." then it produces the
>> correct result.
>> I've checked my program with valgrind and nothing came up.
>>
>> I'm using adolc 2.6.0 in ubuntu 15.10.
>>
>> Here is the example program:
>> //*************************************************************************************************
>>
>>
>> #include <adolc/adolc.h>
>> #include <adolc/sparse/sparsedrivers.h>
>>
>> using namespace std;
>>
>> #define tag 1
>>
>> adouble feval_ad(adouble *x) {
>>     adouble res;
>>
>>     res = 0;
>>     res += x[2] * x[0];
>>     res += x[3] * x[1];
>>
>>     return res;
>> }
>>
>> int main(int argc, char **argv) {
>>     int n = 2;
>>     std::vector<double> pars(2);
>>     int nx = n + pars.size();
>>
>>     double f, x[nx];
>>     adouble fad, xad[nx];
>>     //double f, x[n];
>>     //adouble fad, xad[n];
>>
>>     int i;
>>
>> /****************************************************************************/
>>
>> /*******                function evaluation
>> ***************/
>> /****************************************************************************/
>>
>>
>>     for(i=0;i<n;i++)
>>         x[i] = log(1.0+i);
>>
>>     /* Tracing of function f(x) */
>>
>>     trace_on(tag);
>>     for(i=0;i<n;i++)
>>         xad[i] <<= x[i];
>>
>>     for (i = n; i < nx; ++i) {
>>        xad[i] = ::mkparam(pars[i - n]);
>>     }
>>
>>     fad = feval_ad(xad);
>>
>>     fad >>= f;
>>     trace_off();
>>
>>     printf("\n f = %e\n\n\n",f);
>>
>>
>> /****************************************************************************/
>>
>> /*******       sparse Hessians, complete driver
>> ***************/
>> /****************************************************************************/
>>
>>
>>     /* coordinate format for Hessian */
>>     unsigned int    *rind  = NULL;
>>     unsigned int    *cind  = NULL;
>>     double *values = NULL;
>>     int nnz;
>>     int options[2];
>>
>>     options[0] = 0;          /*                               safe mode
>> (default) */
>>     options[1] = 0;          /*                       indirect recovery
>> (default) */
>>
>>     ::set_param_vec(tag, pars.size(), pars.data());
>>
>>     sparse_hess(tag, n, 0, x, &nnz, &rind, &cind, &values, options);
>>
>>     printf("In sparse format:\n");
>>     for (i=0;i<nnz;i++)
>>         printf("%2d %2d %10.6f\n\n",rind[i],cind[i],values[i]);
>>
>>     free(rind); rind = NULL;
>>     free(cind); cind = NULL;
>>     free(values); values = NULL;
>> }
>> //*************************************************************************************************
>>
>>
>> It prints out:
>>
>>
>>  f = 0.000000e+00
>>
>>
>> In sparse format:
>>  0  0   0.000000
>>
>>  1  1   0.000000
>>
>>
>> But the hessian shouldn't have any elements.
>> Could anyone help me figure out the problem with my code?
>>
>> Thanks for your help!
>>
>> João Leal
>>
>>
>>
>> _______________________________________________
>> ADOL-C mailing list
>> ADOL-C at list.coin-or.org
>> http://list.coin-or.org/mailman/listinfo/adol-c
>>
> 

-- 
Dr. Kshitij Kulshreshtha

Institut für Mathematik,
Universität Paderborn,
Warburger Straße 100,
33098 Paderborn.

Büro: TP21.1.21
Besucheradresse:
Technologiepark 21
33098 Paderborn.

Privatanschrift:
Arnikaweg 62
33100 Paderborn.


More information about the ADOL-C mailing list