Coin-All: /home/coin/SVN-release/CoinAll-1.1.0/cppad/cppad/local/rev_hes

00001 # ifndef CPPAD_REV_HES_SWEEP_INCLUDED
00002 # define CPPAD_REV_HES_SWEEP_INCLUDED
00003 
00004 /* --------------------------------------------------------------------------
00005 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell
00006 
00007 CppAD is distributed under multiple licenses. This distribution is under
00008 the terms of the 
00009                     Common Public License Version 1.0.
00010 
00011 A copy of this license is included in the COPYING file of this distribution.
00012 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00013 -------------------------------------------------------------------------- */
00014 
00015 /*
00016 $begin RevHesSweep$$ $comment CppAD Developer Documentation$$
00017 $spell
00018         Hessian
00019         RevHesSweep
00020         npv
00021         numvar
00022         Num
00023         Var
00024         Op
00025         Taylor
00026         Inv
00027 $$
00028 
00029 $section Reverse Computation of Hessian Sparsity Pattern$$
00030 $index RevHesSweep$$
00031 $index sparsity, reverse Hessian$$
00032 $index reverse, Hessian sparsity$$
00033 $index pattern, reverse Hessian$$
00034 $index bit pattern, Hessian$$
00035 
00036 $head Syntax$$
00037 $syntax%void RevHesSweep(
00038         size_t %npv%,
00039         size_t %numvar%,
00040         const TapeRec<%Base%> *%Rec%,
00041         size_t %TaylorColDim%,
00042         const %Base% *%Taylor%,
00043         const %Pack% *%ForJac%,
00044         const %Pack% *%RevJac%,
00045         %Pack% *%RevHes%
00046 )%$$
00047 
00048 
00049 $head Rec$$
00050 The information stored in $italic Rec$$
00051 is a recording of the operations corresponding to a function
00052 $latex \[
00053         F : B^n \rightarrow B^m
00054 \] $$
00055 
00056 $head Description$$
00057 Given the forward Jacobian sparsity pattern for all the variables
00058 $code RevHesSweep$$ computes the Hessian sparsity pattern
00059 for the dependent variables with respect to the independent variables.
00060 
00061 $head numvar$$
00062 is the number of rows in the entire sparsity patterns
00063 $italic ForJac$$ and $italic RevHes$$.
00064 It must also be equal to $syntax%%Rec%->TotNumVar()%$$.
00065 
00066 $head npv$$
00067 Is the number of elements of type $italic Pack$$
00068 (per variable) in the sparsity patterns 
00069 $itlaic ForJac$$ and $italic RevJac$$.
00070 (The sparsity pattern $italic RevJac$$ only has one element per variables
00071 and all of the bits in that element are the same.)
00072  
00073 $head TaylorColDim$$
00074 Is the number of columns currently stored in the matrix $italic Taylor$$.
00075 
00076 $head Taylor$$
00077 For $latex i = 1 , \ldots , numvar$$,
00078 $syntax%%Taylor%[%i% * %TaylorColDim%]%$$
00079 is the value of the variable with index $italic i$$.
00080 
00081 
00082 $head ForJac$$
00083 For $latex i = 0, \ldots , numvar-1$$ (all the variables),
00084 and $latex j = 0, \ldots , npv-1$$,
00085 $syntax%%ForJac%[%i% * %npv% + %j%]%$$ 
00086 is the $th j$$ subset of the 
00087 the forward mode Jacobian sparsity pattern for 
00088 the variable with index $italic i$$.   
00089 
00090 $head RevJac$$
00091 For $latex i = 0, \ldots , numvar-1$$ (all the variables),
00092 $syntax%%RevJac%[%i%]%$$ 
00093 is all true (ones complement of 0) if the function we are computing
00094 the Hessian for has non-zero Jacobian with respect to variable
00095 with index $italic i$$.
00096 
00097 
00098 $head RevHes$$
00099 For $latex i = 1, \ldots , numvar-1$$ 
00100 and $latex j = 0 , \ldots , npv-1$$,
00101 $syntax%%RevHes%[%i% * %npv% + %j%]%$$ is the 
00102 Hessian sparsity pattern with for the partial of the 
00103 function we are computing the Hessian for 
00104 (usually one variable) with respect to
00105 the $th i$$ variable and with respect to all the variables
00106 corresponding to the columns of $italic ForJac$$.
00107 
00108 $end
00109 ------------------------------------------------------------------------------
00110 */
00111 # define CPPAD_REV_HES_SWEEP_TRACE 0
00112 
00113 
00114 // BEGIN CppAD namespace
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         // used by CExp operator 
00150         const Base *left, *right;
00151         Pack  *trueCaseh;
00152         Pack  *falseCaseh;
00153         Pack  zero(0);
00154 
00155         size_t             j;
00156 
00157         // check numvar argument
00158         CPPAD_ASSERT_UNKNOWN( Rec->TotNumVar() == numvar );
00159         CPPAD_ASSERT_UNKNOWN( numvar > 0 );
00160 
00161         // Initialize
00162         i_op   = Rec->NumOp();
00163         i_var  = Rec->TotNumVar();
00164         i_ind  = Rec->NumInd();
00165         op     = NonOp;         // phony operator that is not there
00166 
00167         while(i_op > 1)
00168         {       --i_op;
00169 
00170                 // next op
00171                 op     = Rec->GetOp(i_op);
00172 
00173                 // corresponding varable
00174                 n_var  = NumVar(op);
00175                 CPPAD_ASSERT_UNKNOWN( i_var >= n_var );
00176                 i_var -= n_var;
00177 
00178                 // corresponding index values
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                 // sparsity for z corresponding to this op
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                 // rest of information depends on the case
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                         // acos(x) and sqrt(1 - x * x) are computed in pairs
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                         // asin(x) and sqrt(1 - x * x) are computed in pairs
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                         // atan(x) and 1 + x * x must be computed in pairs
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                         // cosine and sine must come in pairs
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                         // hyperbolic cosine and sine must come in pairs
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                         // Z is already defined
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                         // ind[2] is variable corresponding to this load
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                         // ind[2] is variable corresponding to this load
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                         // nvar should be zero for this case
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                         // sine and cosine must come in pairs
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                         // sine and cosine must come in pairs
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 } // END CppAD namespace
00731 
00732 # undef CPPAD_REV_HES_SWEEP_TRACE
00733 
00734 # endif
/home/coin/SVN-release/CoinAll-1.1.0/cppad/cppad/local/rev_hes_sweep.hpp
