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