Changeset 1474


Ignore:
Timestamp:
Aug 2, 2009 11:51:37 PM (11 years ago)
Author:
bradbell
Message:

trunk: Compute reverse jacobian and Hessian along side each other.

sparse_vec_ad.cpp: better testing of sparsity pattern for Hessians.
rev_sparse_hes.cpp: better testing of sparsity pattern for Hessians.
prototype_op.hpp: change reverse jacobian to compute along side Hessian.
op_code.hpp: fix printing of load and store parameter indexing.
rev_sparse_hes.hpp: change reverse jacobian to compute along side Hessian.
cond_op.hpp: change reverse jacobian to compute along side Hessian.
rev_hes_sweep.hpp: change reverse jacobian to compute along side Hessian.
sparse_binary_op.hpp: change reverse jacobian to compute along side Hessian.
sparse_unary_op.hpp: change reverse jacobian to compute along side Hessian.

Location:
trunk
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • trunk/cppad/local/cond_op.hpp

    r1469 r1474  
    508508\copydetails sparse_conditional_exp_op
    509509
    510 \param z_jac
    511 is all true (ones complement of 0) if the scalar valued
    512 function we are computing the Hessian sparsity for
    513 has a non-zero partial with respect to the variable z
    514 (actually may have a non-zero partial with respect to z).
    515 Otherwise it zero.
     510\param jac_reverse
     511\a jac_reverse[i_z]
     512is all zero (ones) if the Jacobian of G with respect to z is zero (non-zero).
     513\n
     514\n
     515\a jac_reverse[ arg[4] ]
     516If y_2 is a variable,
     517\a jac_reverse[ arg[4] ]
     518is all zero (ones) if the Jacobian with respect to y_2 is zero (non-zero).
     519On input, it corresponds to the function G,
     520and on output it corresponds to the function H.
     521\n
     522\n
     523\a jac_reverse[ arg[5] ]
     524If y_3 is a variable,
     525\a jac_reverse[ arg[5] ]
     526is all zero (ones) if the Jacobian with respect to y_3 is zero (non-zero).
     527On input, it corresponds to the function G,
     528and on output it corresponds to the function H.
    516529
    517530\param hes_sparsity
     
    522535\n
    523536\n
     537If y_2 is a variable,
    524538For k = 0 , ... , nc_sparsity-1,
    525539\a sparsity [ \a arg[4] * nc_sparsity + k ]
     
    529543\n
    530544\n
     545If y_3 is a variable,
    531546For k = 0 , ... , nc_sparsity-1,
    532547\a sparsity [ \a arg[5] * nc_sparsity + k ]
     
    540555        const size_t*  arg           ,
    541556        size_t         num_par       ,
    542         Pack           z_jac         ,
     557        Pack*          jac_reverse   ,
    543558        size_t         nc_sparsity   ,
    544559        Pack*          hes_sparsity  )
     
    579594                while(k--)
    580595                        y_2_hes[k] |= z_hes[k];
     596
     597                jac_reverse[ arg[4] ] |= jac_reverse[i_z];
    581598        }
    582599        if( arg[1] & 8 )
     
    586603                while(k--)
    587604                        y_3_hes[k] |= z_hes[k];
     605
     606                jac_reverse[ arg[5] ] |= jac_reverse[i_z];
    588607        }
    589608        return;
  • trunk/cppad/local/op_code.hpp

    r1447 r1474  
    432432                CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 );
    433433                printOpField(os, "off=", ind[0], ncol);
    434                 printOpField(os, "  p=", *(Rec->GetPar(ind[1])), ncol);
     434                printOpField(os, "idx=", ind[1], ncol);
    435435                break;
    436436
     
    444444                CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 );
    445445                printOpField(os, "off=", ind[0], ncol);
    446                 printOpField(os, " pl=", *(Rec->GetPar(ind[1])), ncol);
     446                printOpField(os, "idx=", ind[1], ncol);
    447447                printOpField(os, " pr=", *(Rec->GetPar(ind[2])), ncol);
    448448                break;
     
    451451                CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 );
    452452                printOpField(os, "off=", ind[0], ncol);
    453                 printOpField(os, " pl=", *(Rec->GetPar(ind[1])), ncol);
     453                printOpField(os, "idx=", ind[1], ncol);
    454454                printOpField(os, " vr=", ind[2], ncol);
    455455                break;
  • trunk/cppad/local/prototype_op.hpp

    r1468 r1474  
    11621162i.e. the row index in sparsity corresponding to x.
    11631163
    1164 \param jac_z
    1165 is all zero (ones) if the Jacobian of G with respect to z is
    1166 zero (non-zero).
     1164\param jac_reverse
     1165\a jac_reverse[i_z]
     1166is all zero (ones) if the Jacobian of G with respect to z is zero (non-zero).
     1167\n
     1168\n
     1169\a jac_reverse[i_x]
     1170is all zero (ones) if the Jacobian with respect to x is zero (non-zero).
     1171On input, it corresponds to the function G,
     1172and on output it corresponds to the function H.
    11671173
    11681174\param nc_sparsity
     
    11701176the number of columns in the sparsity pattern matrices.
    11711177
    1172 \param jac_sparsity
     1178\param jac_forward
    11731179for k = 0 , ... , \a nc_sparsity - 1,
    1174 jac_sparsity[ \a i_x * \a nc_sparsity + k ]
     1180jac_forward[ \a i_x * \a nc_sparsity + k ]
    11751181is the forward mode Jacobian sparsity pattern for the variable x.
    11761182
    11771183\param hes_sparsity
    1178 \b Input: hes_sparsity[ \a i_z * \a nc_sparsity + k ]
     1184hes_sparsity[ \a i_z * \a nc_sparsity + k ]
    11791185for k = 0 , ... , \a nc_sparsity -1
    11801186is the Hessian sparsity pattern for the fucntion G
    11811187where one of the partials derivative is with respect to z.
    11821188\n
    1183 \b Input: hes_sparsity[ \a i_x * \a nc_sparsity + k ]
     1189\n
     1190hes_sparsity[ \a i_x * \a nc_sparsity + k ]
    11841191for k = 0 , ... , \a nc_sparsity -1
    1185 is the Hessian sparsity pattern for the fucntion G
    1186 where one of the partials derivative is with respect to y.
    1187 \n
    1188 \b Output: hes_sparsity[ \a i_x * \a nc_sparsity + k ]
    1189 for k = 0 , ... , \a nc_sparsity -1
    1190 is the Hessian sparsity pattern for the fucntion H
    1191 where one of the partials derivative is with respect to y.
     1192is the Hessian sparsity pattern
     1193where one of the partials derivative is with respect to x.
     1194On input, it corresonds to the function G,
     1195and on output it corresponds to the function H.
    11921196
    11931197\par Checked Assertions:
     
    11981202        size_t      i_z           ,
    11991203        size_t      i_x           ,
    1200         Pack        jac_z         ,
     1204        Pack*       jac_reverse   ,
    12011205        size_t      nc_sparsity   ,
    1202         const Pack* jac_sparsity  ,
     1206        const Pack* jac_forward   ,
    12031207        Pack*       hes_sparsity  )
    12041208{       
     
    12341238i.e. the row index in sparsity corresponding to y.
    12351239
    1236 \param z_jac
    1237 is all true (ones complement of 0) if the scalar valued
    1238 function we are computing the Hessian sparsity for
    1239 has a non-zero partial with respect to the variable z
    1240 (actually may have a non-zero partial with respect to z).
    1241 Otherwise it zero.
    1242 (not used by add and subtract operators).
     1240\param jac_reverse
     1241\a jac_reverse[i_z]
     1242is all zero (ones) if the Jacobian of G with respect to z is zero (non-zero).
     1243\n
     1244\n
     1245\a jac_reverse[ \a arg[0] ]
     1246is all zero (ones) if the Jacobian with respect to x is zero (non-zero).
     1247On input, it corresponds to the function G,
     1248and on output it corresponds to the function H.
     1249\n
     1250\n
     1251\a jac_reverse[ \a arg[1] ]
     1252is all zero (ones) if the Jacobian with respect to y is zero (non-zero).
     1253On input, it corresponds to the function G,
     1254and on output it corresponds to the function H.
    12431255
    12441256\param nc_sparsity
     
    12461258the number of columns in the sparsity pattern matrix.
    12471259
    1248 \param jac_sparsity
     1260\param jac_forward
    12491261For j = 0 , ... , \a nc_sparsity - 1,
    1250 jac_sparsity[ \a arg[0] * \a nc_sparsity + j ]
     1262jac_forward[ \a arg[0] * \a nc_sparsity + j ]
    12511263is the forward mode Jacobiain sparsity pattern for the variable x.
     1264\n
     1265\n
    12521266For j = 0 , ... , \a nc_sparsity - 1,
    1253 jac_sparsity[ \a arg[1] * \a nc_sparsity + j ]
     1267jac_forward[ \a arg[1] * \a nc_sparsity + j ]
    12541268is the forward mode Jacobiain sparsity pattern for the variable y.
    1255 (These values are not used by add and subtract operators)
    12561269
    12571270\param hes_sparsity
    1258 \b Input: for j = 0 , ... , \a nc_sparsity - 1,
     1271For j = 0 , ... , \a nc_sparsity - 1,
    12591272hes_sparsity [ \a i_z * \a nc_sparsity + j ]
    12601273is the Hessian sparsity pattern for the function G
    12611274where one of the partial derivatives is with respect to z.
    12621275\n
    1263 \b Input: for j = 0 , ... , \a nc_sparsity - 1,
     1276For j = 0 , ... , \a nc_sparsity - 1,
    12641277hes_sparsity [ \a arg[0] * \a nc_sparsity + j ]
    1265 is the Hessian sparsity pattern for the function G
     1278is the Hessian sparsity pattern
    12661279where one of the partial derivatives is with respect to x.
    1267 \n
    1268 \b Input: for j = 0 , ... , \a nc_sparsity - 1,
     1280On input, it corresponds to the function G,
     1281and on output it correspondst to H.
     1282\n
     1283For j = 0 , ... , \a nc_sparsity - 1,
    12691284hes_sparsity [ \a arg[1] * \a nc_sparsity + j ]
    1270 is the Hessian sparsity pattern for the function G
     1285is the Hessian sparsity pattern
    12711286where one of the partial derivatives is with respect to y.
    1272 \n
    1273 \b Output: for j = 0 , ... , \a nc_sparsity - 1,
    1274 hes_sparsity [ \a arg[0] * \a nc_sparsity + j ]
    1275 is the Hessian sparsity pattern for the function H
    1276 where one of the partial derivatives is with respect to x.
    1277 \n
    1278 \b Output: for j = 0 , ... , \a nc_sparsity - 1,
    1279 hes_sparsity [ \a arg[1] * \a nc_sparsity + j ]
    1280 is the Hessian sparsity pattern for the function H
    1281 where one of the partial derivatives is with respect to y.
     1287On input, it corresponds to the function G,
     1288and on output it correspondst to H.
    12821289
    12831290\par Checked Assertions:
     
    12891296        size_t            i_z           ,
    12901297        const size_t*     arg           ,
    1291         Pack              z_jac         ,
     1298        Pack*             jac_reverse   ,
    12921299        size_t            nc_sparsity   ,
    1293         const Pack*       jac_sparsity  ,
     1300        const Pack*       jac_forward   ,
    12941301        Pack*             hes_sparsity  )
    12951302{       
  • trunk/cppad/local/rev_hes_sweep.hpp

    r1473 r1474  
    9090
    9191$head RevJac$$
    92 For $latex i = 0, \ldots , numvar-1$$ (all the variables),
    93 $syntax%%RevJac%[%i%]%$$
     92On input,
     93for $latex i = 1, \ldots , m$$ (the number of dependent variables),
     94$syntax%%RevJac%[%numvar-i%]%$$
    9495is all true (ones complement of 0) if the function we are computing
    95 the Hessian for has non-zero Jacobian with respect to variable
    96 with index $italic i$$.
     96the Hessian for has non-zero Jacobian with respect to the
     97dependent variable with index $italic i$$.
     98On output, this vecotr has been used for temporary storage.
    9799
    98100
     
    124126        const Base           *Taylor,
    125127        const Pack           *ForJac,
    126         const Pack           *RevJac,
     128        Pack                 *RevJac,
    127129        Pack                 *RevHes
    128130)
     
    219221                        CPPAD_ASSERT_NARG_NRES(op, 1, 1)
    220222                        reverse_sparse_hessian_linear_unary_op(
    221                                 i_var, arg[0], *Zr, npv, ForJac, RevHes
     223                                i_var, arg[0], RevJac, npv, ForJac, RevHes
    222224                        );
    223225                        break;
     
    227229                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
    228230                        reverse_sparse_hessian_addsub_op(
    229                                 i_var, arg, *Zr, npv, ForJac, RevHes
     231                                i_var, arg, RevJac, npv, ForJac, RevHes
    230232                        );
    231233                        break;
     
    235237                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
    236238                        reverse_sparse_hessian_linear_unary_op(
    237                                 i_var, arg[1], *Zr, npv, ForJac, RevHes
     239                                i_var, arg[1], RevJac, npv, ForJac, RevHes
    238240                        );
    239241                        break;
     
    243245                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
    244246                        reverse_sparse_hessian_linear_unary_op(
    245                                 i_var, arg[0], *Zr, npv, ForJac, RevHes
     247                                i_var, arg[0], RevJac, npv, ForJac, RevHes
    246248                        );
    247249                        break;
     
    253255                        CPPAD_ASSERT_NARG_NRES(op, 1, 2)
    254256                        reverse_sparse_hessian_nonlinear_unary_op(
    255                                 i_var, arg[0], *Zr, npv, ForJac, RevHes
     257                                i_var, arg[0], RevJac, npv, ForJac, RevHes
    256258                        );
    257259                        break;
     
    263265                        CPPAD_ASSERT_NARG_NRES(op, 1, 2)
    264266                        reverse_sparse_hessian_nonlinear_unary_op(
    265                                 i_var, arg[0], *Zr, npv, ForJac, RevHes
     267                                i_var, arg[0], RevJac, npv, ForJac, RevHes
    266268                        );
    267269                        break;
     
    273275                        CPPAD_ASSERT_NARG_NRES(op, 1, 2)
    274276                        reverse_sparse_hessian_nonlinear_unary_op(
    275                                 i_var, arg[0], *Zr, npv, ForJac, RevHes
     277                                i_var, arg[0], RevJac, npv, ForJac, RevHes
    276278                        );
    277279                        break;
     
    279281                        case CExpOp:
    280282                        reverse_sparse_hessian_cond_op(
    281                                 i_var, arg, num_par, *Zr, npv, RevHes
     283                                i_var, arg, num_par, RevJac, npv, RevHes
    282284                        );
    283285                        break;
     
    295297                        CPPAD_ASSERT_NARG_NRES(op, 1, 2)
    296298                        reverse_sparse_hessian_nonlinear_unary_op(
    297                                 i_var, arg[0], *Zr, npv, ForJac, RevHes
     299                                i_var, arg[0], RevJac, npv, ForJac, RevHes
    298300                        );
    299301                        break;
     
    305307                        CPPAD_ASSERT_NARG_NRES(op, 1, 2)
    306308                        reverse_sparse_hessian_nonlinear_unary_op(
    307                                 i_var, arg[0], *Zr, npv, ForJac, RevHes
     309                                i_var, arg[0], RevJac, npv, ForJac, RevHes
    308310                        );
    309311                        break;
     
    319321                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
    320322                        reverse_sparse_hessian_div_op(
    321                                 i_var, arg, *Zr, npv, ForJac, RevHes
     323                                i_var, arg, RevJac, npv, ForJac, RevHes
    322324                        );
    323325                        break;
     
    327329                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
    328330                        reverse_sparse_hessian_nonlinear_unary_op(
    329                                 i_var, arg[1], *Zr, npv, ForJac, RevHes
     331                                i_var, arg[1], RevJac, npv, ForJac, RevHes
    330332                        );
    331333                        break;
     
    335337                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
    336338                        reverse_sparse_hessian_linear_unary_op(
    337                                 i_var, arg[0], *Zr, npv, ForJac, RevHes
     339                                i_var, arg[0], RevJac, npv, ForJac, RevHes
    338340                        );
    339341                        break;
     
    343345                        CPPAD_ASSERT_NARG_NRES(op, 1, 1)
    344346                        reverse_sparse_hessian_nonlinear_unary_op(
    345                                 i_var, arg[0], *Zr, npv, ForJac, RevHes
     347                                i_var, arg[0], RevJac, npv, ForJac, RevHes
    346348                        );
    347349                        break;
     
    381383                        CPPAD_ASSERT_NARG_NRES(op, 1, 1)
    382384                        reverse_sparse_hessian_nonlinear_unary_op(
    383                                 i_var, arg[0], *Zr, npv, ForJac, RevHes
     385                                i_var, arg[0], RevJac, npv, ForJac, RevHes
    384386                        );
    385387                        break;
     
    389391                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
    390392                        reverse_sparse_hessian_mul_op(
    391                                 i_var, arg, *Zr, npv, ForJac, RevHes
     393                                i_var, arg, RevJac, npv, ForJac, RevHes
    392394                        );
    393395                        break;
     
    397399                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
    398400                        reverse_sparse_hessian_linear_unary_op(
    399                                 i_var, arg[1], *Zr, npv, ForJac, RevHes
     401                                i_var, arg[1], RevJac, npv, ForJac, RevHes
    400402                        );
    401403                        break;
     
    405407                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
    406408                        reverse_sparse_hessian_linear_unary_op(
    407                                 i_var, arg[0], *Zr, npv, ForJac, RevHes
     409                                i_var, arg[0], RevJac, npv, ForJac, RevHes
    408410                        );
    409411                        break;
     
    427429                        CPPAD_ASSERT_NARG_NRES(op, 2, 3)
    428430                        reverse_sparse_hessian_nonlinear_unary_op(
    429                                 i_var+2, arg[1], *(Zr+2), npv, ForJac, RevHes
     431                                i_var+2, arg[1], RevJac, npv, ForJac, RevHes
    430432                        );
    431433                        break;
     
    437439                        CPPAD_ASSERT_NARG_NRES(op, 2, 3)
    438440                        reverse_sparse_hessian_nonlinear_unary_op(
    439                                 i_var+2, arg[0], *(Zr+2), npv, ForJac, RevHes
     441                                i_var+2, arg[0], RevJac, npv, ForJac, RevHes
    440442                        );
    441443                        break;
     
    477479                        CPPAD_ASSERT_NARG_NRES(op, 1, 2)
    478480                        reverse_sparse_hessian_nonlinear_unary_op(
    479                                 i_var, arg[0], *Zr, npv, ForJac, RevHes
     481                                i_var, arg[0], RevJac, npv, ForJac, RevHes
    480482                        );
    481483                        break;
     
    487489                        CPPAD_ASSERT_NARG_NRES(op, 1, 2)
    488490                        reverse_sparse_hessian_nonlinear_unary_op(
    489                                 i_var, arg[0], *Zr, npv, ForJac, RevHes
     491                                i_var, arg[0], RevJac, npv, ForJac, RevHes
    490492                        );
    491493                        break;
     
    495497                        CPPAD_ASSERT_NARG_NRES(op, 1, 1)
    496498                        reverse_sparse_hessian_nonlinear_unary_op(
    497                                 i_var, arg[0], *Zr, npv, ForJac, RevHes
     499                                i_var, arg[0], RevJac, npv, ForJac, RevHes
    498500                        );
    499501                        break;
     
    539541                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
    540542                        reverse_sparse_hessian_addsub_op(
    541                                 i_var, arg, *Zr, npv, ForJac, RevHes
     543                                i_var, arg, RevJac, npv, ForJac, RevHes
    542544                        );
    543545                        break;
     
    547549                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
    548550                        reverse_sparse_hessian_linear_unary_op(
    549                                 i_var, arg[1], *Zr, npv, ForJac, RevHes
     551                                i_var, arg[1], RevJac, npv, ForJac, RevHes
    550552                        );
    551553                        break;
     
    555557                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
    556558                        reverse_sparse_hessian_linear_unary_op(
    557                                 i_var, arg[0], *Zr, npv, ForJac, RevHes
     559                                i_var, arg[0], RevJac, npv, ForJac, RevHes
    558560                        );
    559561                        break;
  • trunk/cppad/local/rev_sparse_hes.hpp

    r1401 r1474  
    259259        }
    260260
    261         // comput the reverse mode Jacobian sparsity
    262         RevJacSweep(1, total_num_var_, &play_, taylor_col_dim_, taylor_, RevJac);
    263 
    264261        // compute the Hessian sparsity patterns
    265262        RevHesSweep(
  • trunk/cppad/local/sparse_binary_op.hpp

    r1466 r1474  
    204204        size_t            i_z           ,
    205205        const size_t*     arg           ,
    206         Pack              z_jac         ,
    207         size_t            nc_sparsity   ,
    208         const Pack*       jac_sparsity  ,
     206        Pack*             jac_reverse   ,
     207        size_t            nc_sparsity   ,
     208        const Pack*       jac_forward   ,
    209209        Pack*             hes_sparsity  )
    210210{       
     
    221221                y_hes[j] |= z_hes[j];
    222222        }
     223
     224        jac_reverse[arg[0]] |= jac_reverse[i_z];
     225        jac_reverse[arg[1]] |= jac_reverse[i_z];
    223226        return;
    224227}       
     
    239242        size_t            i_z           ,
    240243        const size_t*     arg           ,
    241         Pack              z_jac         ,
    242         size_t            nc_sparsity   ,
    243         const Pack*       jac_sparsity  ,
     244        Pack*             jac_reverse         ,
     245        size_t            nc_sparsity   ,
     246        const Pack*       jac_forward   ,
    244247        Pack*             hes_sparsity  )
    245248{       
     
    248251        CPPAD_ASSERT_UNKNOWN( arg[1] < i_z );
    249252
    250         const Pack* x_jac  = jac_sparsity + arg[0] * nc_sparsity;
    251         const Pack* y_jac  = jac_sparsity + arg[1] * nc_sparsity;
     253        const Pack* x_for  = jac_forward + arg[0] * nc_sparsity;
     254        const Pack* y_for  = jac_forward + arg[1] * nc_sparsity;
    252255
    253256        const Pack* z_hes  = hes_sparsity + i_z * nc_sparsity;
     
    257260        size_t j = nc_sparsity;
    258261        while(j--)
    259         {       x_hes[j] |= z_hes[j] | (z_jac & y_jac[j]);
    260                 y_hes[j] |= z_hes[j] | (z_jac & x_jac[j]);
    261         }
     262        {       x_hes[j] |= z_hes[j] | (jac_reverse[i_z] & y_for[j]);
     263                y_hes[j] |= z_hes[j] | (jac_reverse[i_z] & x_for[j]);
     264        }
     265
     266        jac_reverse[arg[0]] |= jac_reverse[i_z];
     267        jac_reverse[arg[1]] |= jac_reverse[i_z];
    262268        return;
    263269}       
     
    278284        size_t            i_z           ,
    279285        const size_t*     arg           ,
    280         Pack              z_jac         ,
    281         size_t            nc_sparsity   ,
    282         const Pack*       jac_sparsity  ,
     286        Pack*             jac_reverse         ,
     287        size_t            nc_sparsity   ,
     288        const Pack*       jac_forward   ,
    283289        Pack*             hes_sparsity  )
    284290{       
     
    287293        CPPAD_ASSERT_UNKNOWN( arg[1] < i_z );
    288294
    289         const Pack* x_jac  = jac_sparsity + arg[0] * nc_sparsity;
    290         const Pack* y_jac  = jac_sparsity + arg[1] * nc_sparsity;
     295        const Pack* x_for  = jac_forward + arg[0] * nc_sparsity;
     296        const Pack* y_for  = jac_forward + arg[1] * nc_sparsity;
    291297
    292298        const Pack* z_hes  = hes_sparsity + i_z * nc_sparsity;
     
    296302        size_t j = nc_sparsity;
    297303        while(j--)
    298         {       x_hes[j] |= z_hes[j] | (z_jac & y_jac[j]);
    299                 y_hes[j] |= z_hes[j] | ( z_jac & (x_jac[j] | y_jac[j]) );
    300         }
     304        {       x_hes[j] |= z_hes[j] | (jac_reverse[i_z] & y_for[j]);
     305                y_hes[j] |= z_hes[j] | ( jac_reverse[i_z] & (x_for[j] | y_for[j]) );
     306        }
     307
     308        jac_reverse[arg[0]] |= jac_reverse[i_z];
     309        jac_reverse[arg[1]] |= jac_reverse[i_z];
    301310        return;
    302311}       
     
    317326        size_t            i_z           ,
    318327        const size_t*     arg           ,
    319         Pack              z_jac         ,
    320         size_t            nc_sparsity   ,
    321         const Pack*       jac_sparsity  ,
     328        Pack*             jac_reverse         ,
     329        size_t            nc_sparsity   ,
     330        const Pack*       jac_forward   ,
    322331        Pack*             hes_sparsity  )
    323332{       
     
    326335        CPPAD_ASSERT_UNKNOWN( arg[1] < i_z );
    327336
    328         const Pack* x_jac  = jac_sparsity + arg[0] * nc_sparsity;
    329         const Pack* y_jac  = jac_sparsity + arg[1] * nc_sparsity;
     337        const Pack* x_for  = jac_forward + arg[0] * nc_sparsity;
     338        const Pack* y_for  = jac_forward + arg[1] * nc_sparsity;
    330339
    331340        const Pack* z_hes  = hes_sparsity + i_z * nc_sparsity;
     
    335344        size_t j = nc_sparsity;
    336345        while(j--)
    337         {       x_hes[j] |= z_hes[j] | ( z_jac & (x_jac[j] | y_jac[j]) );
    338                 y_hes[j] |= z_hes[j] | ( z_jac & (x_jac[j] | y_jac[j]) );
    339         }
     346        {       x_hes[j] |= z_hes[j] | ( jac_reverse[i_z] & (x_for[j] | y_for[j]) );
     347                y_hes[j] |= z_hes[j] | ( jac_reverse[i_z] & (x_for[j] | y_for[j]) );
     348        }
     349
     350        jac_reverse[arg[0]] |= jac_reverse[i_z];
     351        jac_reverse[arg[1]] |= jac_reverse[i_z];
    340352        return;
    341353}       
  • trunk/cppad/local/sparse_unary_op.hpp

    r1464 r1474  
    174174        size_t      i_z           ,
    175175        size_t      i_x           ,
    176         Pack        jac_z         ,
     176        Pack*       jac_reverse   ,
    177177        size_t      nc_sparsity   ,
    178         const Pack* jac_sparsity  ,
     178        const Pack* jac_forward   ,
    179179        Pack*       hes_sparsity  )
    180180{       
     
    182182        CPPAD_ASSERT_UNKNOWN( i_x < i_z );
    183183
    184         Pack* hes_z  = hes_sparsity + i_z * nc_sparsity;
    185         Pack* hes_x  = hes_sparsity + i_x * nc_sparsity;
    186         size_t j = nc_sparsity;
    187         while(j--)
    188                 hes_x[j] |= hes_z[j];
     184        Pack* z_hes  = hes_sparsity + i_z * nc_sparsity;
     185        Pack* x_hes  = hes_sparsity + i_x * nc_sparsity;
     186        size_t j = nc_sparsity;
     187        while(j--)
     188                x_hes[j] |= z_hes[j];
     189
     190        jac_reverse[i_x] |= jac_reverse[i_z];
    189191        return;
    190192}
     
    210212        size_t      i_z           ,
    211213        size_t      i_x           ,
    212         Pack        jac_z         ,
     214        Pack*       jac_reverse   ,
    213215        size_t      nc_sparsity   ,
    214         const Pack* jac_sparsity  ,
     216        const Pack* jac_forward   ,
    215217        Pack*       hes_sparsity  )
    216218{       
     
    218220        CPPAD_ASSERT_UNKNOWN( i_x < i_z );
    219221
    220         const Pack* hes_z  = hes_sparsity + i_z * nc_sparsity;
    221         const Pack* jac_x  = jac_sparsity + i_x * nc_sparsity;
    222         Pack* hes_x         = hes_sparsity + i_x * nc_sparsity;
    223         size_t j = nc_sparsity;
    224         while(j--)
    225                 hes_x[j] |= hes_z[j] | (jac_z & jac_x[j]);
     222        const Pack* z_hes  = hes_sparsity + i_z * nc_sparsity;
     223        const Pack* x_for  = jac_forward  + i_x * nc_sparsity;
     224        Pack* x_hes        = hes_sparsity + i_x * nc_sparsity;
     225        size_t j = nc_sparsity;
     226        while(j--)
     227                x_hes[j] |= z_hes[j] | (jac_reverse[i_z] & x_for[j]);
     228
     229        jac_reverse[i_x] |= jac_reverse[i_z];
    226230        return;
    227231}
  • trunk/test_more/rev_sparse_hes.cpp

    r1469 r1474  
    167167}
    168168
     169bool case_three(void)
     170{       bool ok = true;
     171        using CppAD::AD;
     172
     173        // domain space vector
     174        size_t n = 1;
     175        CPPAD_TEST_VECTOR< AD<double> > X(n);
     176        X[0] = 0.;
     177
     178        // declare independent variables and start recording
     179        CppAD::Independent(X);
     180
     181        // range space vector
     182        size_t m = 1;
     183        CPPAD_TEST_VECTOR< AD<double> > Y(m);
     184
     185        // make sure reverse jacobian is propogating dependency to
     186        // intermediate values (not just final ones).
     187        Y[0] = X[0] * X[0] + 2;
     188
     189        // create f: X -> Y and stop tape recording
     190        CppAD::ADFun<double> f(X, Y);
     191
     192        // sparsity pattern for the identity matrix
     193        CppAD::vector<bool> r(n * n);
     194        size_t i, j;
     195        for(i = 0; i < n; i++)
     196        {       for(j = 0; j < n; j++)
     197                        r[ i * n + j ] = false;
     198                r[ i * n + i ] = true;
     199        }
     200
     201        // compute sparsity pattern for J(x) = F^{(1)} (x)
     202        f.ForSparseJac(n, r);
     203
     204        // compute sparsity pattern for H(x) = F_0^{(2)} (x)
     205        CppAD::vector<bool> s(m);
     206        for(i = 0; i < m; i++)
     207                s[i] = false;
     208        s[0] = true;
     209        CppAD::vector<bool> h(n * n);
     210        h    = f.RevSparseHes(n, s);
     211
     212        // check values
     213        ok  &= (h[ 0 * n + 0 ] == true);  // second partial w.r.t x[0], x[0]
     214
     215        return ok;
     216}
     217
    169218} // End of empty namespace
    170219
    171220bool rev_sparse_hes(void)
    172221{       bool ok = true;
    173         // ok &= case_one();
     222        ok &= case_one();
    174223        ok &= case_two();
    175 
    176         return ok;
    177 }
     224        ok &= case_three();
     225
     226        return ok;
     227}
  • trunk/test_more/sparse_vec_ad.cpp

    r1472 r1474  
    5050                // y[i] depends on x[j] for j <= i
    5151                // (and is non-linear for j <= 1).
    52                 if( j <= 1 )
     52                if( j == 1 )
    5353                        Y[j] = Z[J] * Z[J];
    5454                else    Y[j] = Z[J];
     
    9393        }       
    9494
    95         // compute sparsity pattern for Hessian of F_m ( Identity(x) )
     95        // test sparsity pattern for Hessian of F_2 ( Identity(x) )
    9696        CPPAD_TEST_VECTOR<bool> Hy(m);
    9797        for(i = 0; i < m; i++)
    9898                Hy[i] = false;
    99         Hy[m-1] = true;
     99        Hy[2] = true;
    100100        CPPAD_TEST_VECTOR<bool> Pxx(n * n);
    101101        Pxx = F.RevSparseHes(n, Hy);
     102        for(i = 0; i < n; i++)
     103        {       for(j = 0; j < n; j++)
     104                        ok &= (Pxx[i * n + j] == false );
     105        }
    102106
    103 # if 0
    104         // This test case demonstrates a bug in the current version of
    105         // of the reverse Hessian sparsity pattern calculations. It will be
    106         // fixed in the next commit.
     107        // test sparsity pattern for Hessian of F_1 ( Identity(x) )
     108        for(i = 0; i < m; i++)
     109                Hy[i] = false;
     110        Hy[1] = true;
     111        Pxx = F.RevSparseHes(n, Hy);
    107112        for(i = 0; i < n; i++)
    108113        {       for(j = 0; j < n; j++)
    109114                        ok &= (Pxx[i * n + j] == ( (i <= 1) & (j <= 1) ) );
    110115        }
    111 # endif
     116
    112117
    113118        return ok;
Note: See TracChangeset for help on using the changeset viewer.