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