00001 # ifndef CPPAD_REV_HES_SWEEP_INCLUDED
00002 # define CPPAD_REV_HES_SWEEP_INCLUDED
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111 # define CPPAD_REV_HES_SWEEP_TRACE 0
00112
00113
00114
00115 namespace CppAD {
00116
00117 template <class Base, class Pack>
00118 void RevHesSweep(
00119 size_t npv,
00120 size_t numvar,
00121 const TapeRec<Base> *Rec,
00122 size_t TaylorColDim,
00123 const Base *Taylor,
00124 const Pack *ForJac,
00125 const Pack *RevJac,
00126 Pack *RevHes
00127 )
00128 {
00129 OpCode op;
00130 size_t i_op;
00131 size_t i_var;
00132 size_t i_ind;
00133 size_t n_var;
00134 size_t n_ind;
00135
00136 const size_t *ind;
00137
00138 const Pack *Xf;
00139 Pack *Xh;
00140
00141 const Pack *Yf;
00142 Pack *Yh;
00143
00144 const Pack *Zf;
00145 const Pack *Zr;
00146 const Pack *Zh;
00147
00148
00149
00150 const Base *left, *right;
00151 Pack *trueCaseh;
00152 Pack *falseCaseh;
00153 Pack zero(0);
00154
00155 size_t j;
00156
00157
00158 CPPAD_ASSERT_UNKNOWN( Rec->TotNumVar() == numvar );
00159 CPPAD_ASSERT_UNKNOWN( numvar > 0 );
00160
00161
00162 i_op = Rec->NumOp();
00163 i_var = Rec->TotNumVar();
00164 i_ind = Rec->NumInd();
00165 op = NonOp;
00166
00167 while(i_op > 1)
00168 { --i_op;
00169
00170
00171 op = Rec->GetOp(i_op);
00172
00173
00174 n_var = NumVar(op);
00175 CPPAD_ASSERT_UNKNOWN( i_var >= n_var );
00176 i_var -= n_var;
00177
00178
00179 n_ind = NumInd(op);
00180 CPPAD_ASSERT_UNKNOWN( i_ind >= n_ind );
00181 i_ind -= n_ind;
00182 ind = Rec->GetInd(n_ind, i_ind);
00183
00184
00185 Zf = ForJac + i_var * npv;
00186 Zr = RevJac + i_var;
00187 Zh = RevHes + i_var * npv;
00188
00189 # if CPPAD_REV_HES_SWEEP_TRACE
00190 printOp(
00191 std::cout,
00192 Rec,
00193 i_var,
00194 op,
00195 ind,
00196 npv,
00197 Zf,
00198 npv,
00199 Zh
00200 );
00201 # endif
00202
00203
00204 switch( op )
00205 {
00206 case AbsOp:
00207 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00208 CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00209 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00210
00211 Xh = RevHes + ind[0] * npv;
00212 for(j = 0; j < npv; j++)
00213 { Xh[j] |= Zh[j];
00214 }
00215 break;
00216
00217
00218 case AddvvOp:
00219 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00220 CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00221 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00222 CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00223
00224 Xh = RevHes + ind[0] * npv;
00225 Yh = RevHes + ind[1] * npv;
00226 for(j = 0; j < npv; j++)
00227 { Xh[j] |= Zh[j];
00228 Yh[j] |= Zh[j];
00229 }
00230 break;
00231
00232
00233 case AddpvOp:
00234 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00235 CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00236 CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00237
00238 Yh = RevHes + ind[1] * npv;
00239 for(j = 0; j < npv; j++)
00240 { Yh[j] |= Zh[j];
00241 }
00242 break;
00243
00244
00245 case AddvpOp:
00246 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00247 CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00248 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00249
00250 Xh = RevHes + ind[0] * npv;
00251 for(j = 0; j < npv; j++)
00252 { Xh[j] |= Zh[j];
00253 }
00254 break;
00255
00256
00257 case AcosOp:
00258 CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00259 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00260
00261
00262 CPPAD_ASSERT_UNKNOWN( n_var == 2);
00263 CPPAD_ASSERT_UNKNOWN( (i_var+1) < numvar );
00264
00265 Xf = ForJac + ind[0] * npv;
00266 Xh = RevHes + ind[0] * npv;
00267 for(j = 0; j < npv; j++)
00268 { Xh[j] |= Zh[j] | (*Zr & Xf[j]);
00269 }
00270 break;
00271
00272
00273 case AsinOp:
00274 CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00275 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00276
00277
00278 CPPAD_ASSERT_UNKNOWN( n_var == 2);
00279 CPPAD_ASSERT_UNKNOWN( (i_var+1) < numvar );
00280
00281
00282 Xf = ForJac + ind[0] * npv;
00283 Xh = RevHes + ind[0] * npv;
00284 for(j = 0; j < npv; j++)
00285 { Xh[j] |= Zh[j] | (*Zr & Xf[j]);
00286 }
00287 break;
00288
00289
00290 case AtanOp:
00291 CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00292 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00293
00294
00295 CPPAD_ASSERT_UNKNOWN( n_var == 2);
00296 CPPAD_ASSERT_UNKNOWN( (i_var+1) < numvar );
00297
00298 Xf = ForJac + ind[0] * npv;
00299 Xh = RevHes + ind[0] * npv;
00300 for(j = 0; j < npv; j++)
00301 { Xh[j] |= Zh[j] | (*Zr & Xf[j]);
00302 }
00303 break;
00304
00305 case CExpOp:
00306 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00307 CPPAD_ASSERT_UNKNOWN( n_ind == 6);
00308 CPPAD_ASSERT_UNKNOWN( ind[1] != 0 );
00309
00310 if( ind[1] & 1 )
00311 left = Taylor + ind[2] * TaylorColDim;
00312 else left = Rec->GetPar(ind[2]);
00313 if( ind[1] & 2 )
00314 right = Taylor + ind[3] * TaylorColDim;
00315 else right = Rec->GetPar(ind[3]);
00316 if( ind[1] & 4 )
00317 { trueCaseh = RevHes + ind[4] * npv;
00318 for(j = 0; j < npv; j++)
00319 { trueCaseh[j] |= CondExpTemplate(
00320 CompareOp( ind[0] ),
00321 *left,
00322 *right,
00323 Zh[j],
00324 zero
00325 );
00326 }
00327 }
00328 if( ind[1] & 8 )
00329 { falseCaseh = RevHes + ind[5] * npv;
00330 for(j = 0; j < npv; j++)
00331 { falseCaseh[j] |= CondExpTemplate(
00332 CompareOp( ind[0] ),
00333 *left,
00334 *right,
00335 zero,
00336 Zh[j]
00337 );
00338 }
00339 }
00340 break;
00341
00342
00343 case ComOp:
00344 CPPAD_ASSERT_UNKNOWN( n_var == 0 );
00345 CPPAD_ASSERT_UNKNOWN( n_ind == 4 );
00346 CPPAD_ASSERT_UNKNOWN( ind[1] > 1 );
00347 break;
00348
00349
00350 case CosOp:
00351 CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00352 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00353
00354
00355 CPPAD_ASSERT_UNKNOWN( n_var == 2);
00356 CPPAD_ASSERT_UNKNOWN( (i_var+1) < numvar );
00357
00358 Xf = ForJac + ind[0] * npv;
00359 Xh = RevHes + ind[0] * npv;
00360 for(j = 0; j < npv; j++)
00361 { Xh[j] |= Zh[j] | (*Zr & Xf[j]);
00362 }
00363 break;
00364
00365
00366 case CoshOp:
00367 CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00368 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00369
00370
00371 CPPAD_ASSERT_UNKNOWN( n_var == 2);
00372 CPPAD_ASSERT_UNKNOWN( (i_var+1) < numvar );
00373
00374
00375 Xf = ForJac + ind[0] * npv;
00376 Xh = RevHes + ind[0] * npv;
00377 for(j = 0; j < npv; j++)
00378 { Xh[j] |= Zh[j] | (*Zr & Xf[j]);
00379 }
00380 break;
00381
00382
00383 case DisOp:
00384 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00385 CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00386
00387 break;
00388
00389
00390 case DivvvOp:
00391 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00392 CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00393 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00394 CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00395
00396 Xf = ForJac + ind[0] * npv;
00397 Xh = RevHes + ind[0] * npv;
00398 Yf = ForJac + ind[1] * npv;
00399 Yh = RevHes + ind[1] * npv;
00400 for(j = 0; j < npv; j++)
00401 { Xh[j] |= Zh[j] | (*Zr & Yf[j]);
00402 Yh[j] |= Zh[j] | (*Zr & (Xf[j] | Yf[j]));
00403 }
00404 break;
00405
00406
00407 case DivpvOp:
00408 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00409 CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00410 CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00411
00412 Yf = ForJac + ind[1] * npv;
00413 Yh = RevHes + ind[1] * npv;
00414 for(j = 0; j < npv; j++)
00415 { Yh[j] |= Zh[j] | (*Zr & Yf[j]);
00416 }
00417 break;
00418
00419
00420 case DivvpOp:
00421 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00422 CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00423 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00424
00425 Xh = RevHes + ind[0] * npv;
00426 for(j = 0; j < npv; j++)
00427 { Xh[j] |= Zh[j];
00428 }
00429 break;
00430
00431
00432 case ExpOp:
00433 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00434 CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00435 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00436
00437
00438 Xf = ForJac + ind[0] * npv;
00439 Xh = RevHes + ind[0] * npv;
00440 for(j = 0; j < npv; j++)
00441 { Xh[j] |= Zh[j] | (*Zr & Xf[j]);
00442 }
00443 break;
00444
00445
00446 case InvOp:
00447 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00448 CPPAD_ASSERT_UNKNOWN( n_ind == 0 );
00449
00450 break;
00451
00452
00453 case LdpOp:
00454 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00455 CPPAD_ASSERT_UNKNOWN( n_ind == 3 );
00456
00457 CPPAD_ASSERT_UNKNOWN( ind[0] > 0 );
00458 CPPAD_ASSERT_UNKNOWN( ind[0] < Rec->NumVecInd() );
00459
00460
00461 if( ind[2] > 0 )
00462 { Xh = RevHes + ind[2] * npv;
00463 for(j = 0; j < npv; j++)
00464 Xh[j] |= Zh[j];
00465 }
00466 break;
00467
00468
00469 case LdvOp:
00470 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00471 CPPAD_ASSERT_UNKNOWN( n_ind == 3 );
00472
00473 CPPAD_ASSERT_UNKNOWN( ind[0] > 0 );
00474 CPPAD_ASSERT_UNKNOWN( ind[0] < Rec->NumVecInd() );
00475
00476
00477 if( ind[2] > 0 )
00478 { Xh = RevHes + ind[2] * npv;
00479 for(j = 0; j < npv; j++)
00480 Xh[j] |= Zh[j];
00481 }
00482 break;
00483
00484
00485 case LogOp:
00486 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00487 CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00488 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00489
00490 Xf = ForJac + ind[0] * npv;
00491 Xh = RevHes + ind[0] * npv;
00492 for(j = 0; j < npv; j++)
00493 { Xh[j] |= Zh[j] | (*Zr & Xf[j]);
00494 }
00495 break;
00496
00497
00498 case MulvvOp:
00499 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00500 CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00501 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00502 CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00503
00504 Xf = ForJac + ind[0] * npv;
00505 Xh = RevHes + ind[0] * npv;
00506 Yf = ForJac + ind[1] * npv;
00507 Yh = RevHes + ind[1] * npv;
00508 for(j = 0; j < npv; j++)
00509 { Xh[j] |= Zh[j] | (*Zr & Yf[j]);
00510 Yh[j] |= Zh[j] | (*Zr & Xf[j]);
00511 }
00512 break;
00513
00514
00515 case MulpvOp:
00516 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00517 CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00518 CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00519
00520 Yh = RevHes + ind[1] * npv;
00521 for(j = 0; j < npv; j++)
00522 { Yh[j] |= Zh[j];
00523 }
00524 break;
00525
00526
00527 case MulvpOp:
00528 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00529 CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00530 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00531
00532 Xh = RevHes + ind[0] * npv;
00533 for(j = 0; j < npv; j++)
00534 { Xh[j] |= Zh[j];
00535 }
00536 break;
00537
00538
00539 case NonOp:
00540 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00541 CPPAD_ASSERT_UNKNOWN( n_ind == 0 );
00542
00543 break;
00544
00545
00546 case ParOp:
00547 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00548 CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00549
00550 break;
00551
00552
00553 case PowpvOp:
00554 CPPAD_ASSERT_UNKNOWN( n_var == 3 );
00555 CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00556 CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00557
00558 Yf = ForJac + ind[1] * npv;
00559 Yh = RevHes + ind[1] * npv;
00560 for(j = 0; j < npv; j++)
00561 { Xh[j] |= Zh[j] | (*Zr & Yf[j]);
00562 Yh[j] |= Zh[j] | (*Zr & Xf[j]);
00563 }
00564 break;
00565
00566
00567 case PowvpOp:
00568 CPPAD_ASSERT_UNKNOWN( n_var == 3 );
00569 CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00570 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00571
00572 Xf = ForJac + ind[0] * npv;
00573 Xh = RevHes + ind[0] * npv;
00574 for(j = 0; j < npv; j++)
00575 { Xh[j] |= Zh[j] | (*Zr & Xf[j]);
00576 Yh[j] |= Zh[j] | (*Zr & Xf[j]);
00577 }
00578 break;
00579
00580
00581 case PowvvOp:
00582 CPPAD_ASSERT_UNKNOWN( n_var == 3 );
00583 CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00584 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00585 CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00586
00587 Xf = ForJac + ind[0] * npv;
00588 Xh = RevHes + ind[0] * npv;
00589 Yf = ForJac + ind[1] * npv;
00590 Yh = RevHes + ind[1] * npv;
00591 for(j = 0; j < npv; j++)
00592 { Xh[j] |= Zh[j] | (*Zr & (Xf[j] | Yf[j]));
00593 Yh[j] |= Zh[j] | (*Zr & (Xf[j] | Yf[j]));
00594 }
00595 break;
00596
00597
00598 case PripOp:
00599 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00600
00601 break;
00602
00603
00604 case PrivOp:
00605
00606 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00607
00608 break;
00609
00610
00611 case SinOp:
00612 CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00613 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00614
00615
00616 CPPAD_ASSERT_UNKNOWN( n_var == 2);
00617 CPPAD_ASSERT_UNKNOWN( (i_var+1) < numvar );
00618
00619 Xf = ForJac + ind[0] * npv;
00620 Xh = RevHes + ind[0] * npv;
00621 for(j = 0; j < npv; j++)
00622 { Xh[j] |= Zh[j] | (*Zr & Xf[j]);
00623 }
00624 break;
00625
00626
00627 case SinhOp:
00628 CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00629 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00630
00631
00632 CPPAD_ASSERT_UNKNOWN( n_var == 2);
00633 CPPAD_ASSERT_UNKNOWN( (i_var+1) < numvar );
00634
00635 Xf = ForJac + ind[0] * npv;
00636 Xh = RevHes + ind[0] * npv;
00637 for(j = 0; j < npv; j++)
00638 { Xh[j] |= Zh[j] | (*Zr & Xf[j]);
00639 }
00640 break;
00641
00642
00643 case SqrtOp:
00644 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00645 CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00646 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00647
00648 Xf = ForJac + ind[0] * npv;
00649 Xh = RevHes + ind[0] * npv;
00650 for(j = 0; j < npv; j++)
00651 { Xh[j] |= Zh[j] | (*Zr & Xf[j]);
00652 }
00653 break;
00654
00655
00656 case StppOp:
00657 CPPAD_ASSERT_UNKNOWN( n_var == 0);
00658 CPPAD_ASSERT_UNKNOWN( n_ind == 3 );
00659 break;
00660
00661
00662 case StpvOp:
00663 CPPAD_ASSERT_UNKNOWN( n_var == 0);
00664 CPPAD_ASSERT_UNKNOWN( n_ind == 3 );
00665 break;
00666
00667
00668 case StvpOp:
00669 CPPAD_ASSERT_UNKNOWN( n_var == 0);
00670 CPPAD_ASSERT_UNKNOWN( n_ind == 3 );
00671 break;
00672
00673
00674 case StvvOp:
00675 CPPAD_ASSERT_UNKNOWN( n_var == 0);
00676 CPPAD_ASSERT_UNKNOWN( n_ind == 3 );
00677 break;
00678
00679
00680 case SubvvOp:
00681 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00682 CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00683 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00684 CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00685
00686 Xh = RevHes + ind[0] * npv;
00687 Yh = RevHes + ind[1] * npv;
00688 for(j = 0; j < npv; j++)
00689 { Xh[j] |= Zh[j];
00690 Yh[j] |= Zh[j];
00691 }
00692 break;
00693
00694
00695 case SubpvOp:
00696 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00697 CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00698 CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00699
00700 Yh = RevHes + ind[1] * npv;
00701 for(j = 0; j < npv; j++)
00702 { Yh[j] |= Zh[j];
00703 }
00704 break;
00705
00706
00707 case SubvpOp:
00708 CPPAD_ASSERT_UNKNOWN( n_var == 1);
00709 CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00710 CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00711
00712 Xh = RevHes + ind[0] * npv;
00713 for(j = 0; j < npv; j++)
00714 { Xh[j] |= Zh[j];
00715 }
00716 break;
00717
00718
00719 default:
00720 CPPAD_ASSERT_UNKNOWN(0);
00721 }
00722 }
00723 CPPAD_ASSERT_UNKNOWN( i_op == 1 );
00724 CPPAD_ASSERT_UNKNOWN( Rec->GetOp(i_op-1) == NonOp );
00725 CPPAD_ASSERT_UNKNOWN( i_var == NumVar(NonOp) );
00726
00727 return;
00728 }
00729
00730 }
00731
00732 # undef CPPAD_REV_HES_SWEEP_TRACE
00733
00734 # endif