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