Coin-All: /home/coin/svn-release/OptimizationSuite-1.1.0/cppad/cppad/local/forward

00001 # ifndef CPPAD_FORWARD_SWEEP_INCLUDED
00002 # define CPPAD_FORWARD_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 ForwardSweep$$ $comment CppAD Developer Documentation$$
00017 $spell
00018         Var
00019         numvar
00020         bool
00021         Prip
00022         Priv
00023         Inv
00024         Num
00025         Len
00026         const
00027         Taylor
00028         CppAD
00029         zs
00030         op
00031         Ind
00032 $$
00033 
00034 $section Forward Computation of Taylor Coefficients$$
00035 $index ForwardSweep$$
00036 $index mode, forward$$
00037 $index forward, mode$$
00038 $index derivative, forward$$
00039 $index Taylor coefficient, forward$$
00040 
00041 
00042 $head Syntax$$
00043 $syntax%size_t ForwardSweep(
00044         bool %print%,
00045         size_t %d%,
00046         size_t %numvar%,
00047         TapeRec<%Base%> *%Rec%,
00048         size_t %J%,
00049         Base *%Taylor%,
00050 )%$$
00051 
00052 
00053 $head Return Value$$
00054 The return value is equal to the number of
00055 $syntax%AD<%Base%>%$$ comparison operations have a different result
00056 from when the information in $italic Rec$$ was recorded.
00057 
00058 
00059 $head Rec$$
00060 The information stored in $italic Rec$$
00061 is a recording of the operations corresponding to a function
00062 $latex \[
00063         F : B^n \rightarrow B^m
00064 \] $$
00065 
00066 $head print$$
00067 If $italic print$$ is false,
00068 suppress the output that is otherwise generated 
00069 by the PripOp and PrivOp instructions.
00070 
00071 
00072 $head d$$
00073 Given the $th d-1$$ order Taylor coefficients matrix for all the variables,
00074 and the $th d$$ order Taylor coefficients for all the independent variables,
00075 $code ForwardSweep$$ computes the $th d$$ order Taylor coefficients 
00076 for all the other variables.
00077 
00078 
00079 $head numvar$$
00080 is the number of rows in the matrix $italic Taylor$$.
00081 It must also be equal to $syntax%%Rec%->TotNumVar()%$$.
00082 
00083 
00084 $head J$$
00085 Is the number of columns in the coefficient matrix $italic Taylor$$.
00086 This must be greater than or equal $latex d+1$$.
00087 
00088 
00089 $head On Input$$
00090 
00091 $subhead Independent Variables and Operators$$
00092 The independent variable records come first.
00093 For $latex i = 1, \ldots , n$$ and $latex j = 0 , \ldots , d$$,
00094 $table
00095         $bold field$$ $cnext 
00096         $bold Value$$          
00097 $rnext
00098         $syntax%%Taylor%[%0% * %J% + %j%]%$$      $cnext 
00099         the variable with index zero is not used
00100 $rnext
00101         $syntax%%Rec%->GetOp(0)%$$                $cnext 
00102         the operator with index zero must be a $code NonOp$$
00103 $rnext
00104         $syntax%%Taylor%[%i% * %J% + %j%]%$$      $cnext 
00105         $th j$$ order coefficient for variable with index $italic i$$   
00106 $rnext
00107         $syntax%%Rec%->GetOp(%i%)%$$              $cnext 
00108         the operator with index $italic i$$ must be a $code InvOp$$
00109 $tend
00110 
00111 $subhead Other Variables and Operators$$
00112 The other variables follow the independent variables.
00113 For $latex i = n+1, \ldots , numvar-1$$,
00114 $latex j = 0 , \ldots , d-1$$,
00115 and $latex k = n+1, \ldots ,$$ $syntax%%Rec%->NumOp() - 1%$$,
00116 $table
00117         $bold field$$ $cnext 
00118         $bold Value$$          
00119 $rnext
00120         $syntax%%Taylor%[%i% * %J% + %j%]%$$      $cnext 
00121         $th j$$ coefficient for variable with index $italic i$$     
00122 $rnext
00123         $syntax%%Rec%->GetOp(%i%)%$$              $cnext 
00124         any operator except for $code InvOp$$ 
00125 $tend
00126 
00127 $head On Output$$
00128 
00129 $subhead Rec$$
00130 None of the values stored in $italic Rec$$ are modified.
00131 
00132 $subhead Independent Variables$$
00133 For $latex i = 1, \ldots , n$$ and $latex j = 0 , \ldots , J-1$$,
00134 $syntax%%Taylor%[%i% * %J% + %j%]%$$ is not modified.
00135 
00136 
00137 $subhead Other Variables$$
00138 For $latex i = n+1, \ldots , numvar-1$$ and $latex j < d$$,
00139 $syntax%%Taylor%[%i% * %J% + %j%]%$$ is not modified.
00140 For $latex i = n+1, \ldots , numvar-1$$, 
00141 $syntax%%Taylor%[%i% * %J% + %d%]%$$ is set equal to the
00142 $th d$$ order Taylor coefficient for the variable with index $italic i$$.
00143 
00144 
00145 $end
00146 ------------------------------------------------------------------------------
00147 */
00148 # define CPPAD_FORWARD_SWEEP_TRACE 0
00149 
00150 // BEGIN CppAD namespace
00151 namespace CppAD {
00152 
00153 template <class Base>
00154 size_t ForwardSweep(
00155         bool                  print,
00156         size_t                d,
00157         size_t                numvar,
00158         TapeRec<Base>        *Rec,
00159         size_t                J,
00160         Base                 *Taylor
00161 )
00162 {
00163         size_t        numop;
00164         OpCode           op;
00165         size_t         i_op;
00166         size_t        i_var;
00167         size_t        i_ind;
00168         size_t        n_var;
00169         size_t        n_ind;
00170 
00171         const size_t   *ind;
00172         const Base       *P = 0;
00173         const Base       *X = 0;
00174         const Base       *Y = 0;
00175 
00176         // used by CExp operator (left and right also used by Com operator)
00177         const Base  *left, *right, *trueCase, *falseCase;
00178         const Base  zero = Base(0);
00179 
00180         // used by Com operator
00181         bool result;
00182 
00183         Base             *Z = 0;
00184         Base             *W = 0;
00185         Base             *U = 0;
00186 
00187         size_t            i;
00188         size_t          len;
00189 
00190 
00191         // initialize the comparision operator (ComOp) counter
00192         size_t compareCount = 0;
00193 
00194         // if this is an order zero calculation, initialize vector indices
00195         size_t *VectorInd = CPPAD_NULL;  // address for each element
00196         bool   *VectorVar = CPPAD_NULL;  // is element a variable
00197         i = Rec->NumVecInd();
00198         if( i > 0 )
00199         {       VectorInd = CPPAD_TRACK_NEW_VEC(i, VectorInd);
00200                 VectorVar = CPPAD_TRACK_NEW_VEC(i, VectorVar);
00201                 while(i--)
00202                 {       VectorInd[i] = Rec->GetVecInd(i);
00203                         VectorVar[i] = false;
00204                 }
00205         }
00206 
00207 
00208         // check numvar argument
00209         CPPAD_ASSERT_UNKNOWN( Rec->TotNumVar() == numvar );
00210 
00211         // set the number of operators
00212         numop = Rec->NumOp();
00213 
00214         // skip the NonOp at the beginning of the recording
00215         i_op  = 0;
00216         i_var = 0;
00217         i_ind = 0;
00218         op    = Rec->GetOp(i_op);
00219         n_var = NumVar(op);
00220         n_ind = NumInd(op);
00221         CPPAD_ASSERT_UNKNOWN( op == NonOp );
00222         CPPAD_ASSERT_UNKNOWN( n_var == 1 );
00223         CPPAD_ASSERT_UNKNOWN( n_ind == 0 );
00224 
00225         while(++i_op < numop)
00226         {
00227                 // increment for previous op
00228                 i_var += n_var;
00229                 i_ind += n_ind;
00230 
00231                 // this op
00232                 op     = Rec->GetOp(i_op);
00233 
00234                 // number of variables
00235                 n_var  = NumVar(op);
00236 
00237                 // index field values for this op
00238                 n_ind  = NumInd(op);
00239                 ind    = Rec->GetInd(n_ind, i_ind);
00240 
00241                 // value of z for this op
00242                 Z      = Taylor + i_var * J;
00243 
00244                 // rest of information depends on the case
00245 
00246                 switch( op )
00247                 {
00248                         case AbsOp:
00249                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00250                         CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00251                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00252                         X   = Taylor + ind[0] * J;
00253                         ForAbsOp(d, Z, X);
00254                         break;
00255                         // -------------------------------------------------
00256 
00257                         case AddvvOp:
00258                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00259                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00260                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00261                         CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00262 
00263                         X = Taylor + ind[0] * J;
00264                         Y = Taylor + ind[1] * J;
00265                         ForAddvvOp(d, Z, X, Y);
00266                         break;
00267                         // -------------------------------------------------
00268 
00269                         case AddpvOp:
00270                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00271                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00272                         CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00273 
00274                         P = Rec->GetPar( ind[0] );
00275                         Y = Taylor + ind[1] * J;
00276                         ForAddpvOp(d, Z, P, Y);
00277                         break;
00278                         // -------------------------------------------------
00279 
00280                         case AddvpOp:
00281                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00282                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00283                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00284 
00285                         X = Taylor + ind[0] * J;
00286                         P = Rec->GetPar( ind[1] );
00287                         ForAddvpOp(d, Z, X, P);
00288                         break;
00289                         // -------------------------------------------------
00290 
00291                         case AcosOp:
00292                         CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00293                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00294 
00295                         // acos(x) and sqrt(1 - x * x) are computed in pairs
00296                         CPPAD_ASSERT_UNKNOWN( n_var == 2);
00297                         CPPAD_ASSERT_UNKNOWN( (i_var+1) < numvar  );
00298 
00299                         // use W for data stored in variable record
00300                         W = Taylor + (i_var+1) * J;
00301                         X   = Taylor + ind[0] * J;
00302                         ForAcosOp(d, Z, W, X);
00303                         break;
00304                         // -------------------------------------------------
00305 
00306                         case AsinOp:
00307                         CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00308                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00309 
00310                         // asin(x) and sqrt(1 - x * x) are computed in pairs
00311                         CPPAD_ASSERT_UNKNOWN( n_var == 2);
00312                         CPPAD_ASSERT_UNKNOWN( (i_var+1) < numvar  );
00313 
00314                         // use W for data stored in variable record
00315                         W = Taylor + (i_var+1) * J;
00316                         X   = Taylor + ind[0] * J;
00317                         ForAsinOp(d, Z, W, X);
00318                         break;
00319                         // -------------------------------------------------
00320 
00321                         case AtanOp:
00322                         CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00323                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00324 
00325                         // atan(x) and 1 + x * x must be computed in pairs
00326                         CPPAD_ASSERT_UNKNOWN( n_var == 2);
00327                         CPPAD_ASSERT_UNKNOWN( (i_var+1) < numvar  );
00328 
00329                         // use W for data stored in variable record
00330                         W = Taylor + (i_var+1) * J;
00331                         X   = Taylor + ind[0] * J;
00332                         ForAtanOp(d, Z, W, X);
00333                         break;
00334                         // -------------------------------------------------
00335 
00336                         case CExpOp:
00337                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00338                         CPPAD_ASSERT_UNKNOWN( n_ind == 6);
00339                         CPPAD_ASSERT_UNKNOWN( ind[1] != 0 );
00340                         if( ind[1] & 1 )
00341                                 left = Taylor + ind[2] * J;
00342                         else    left = Rec->GetPar(ind[2]);
00343                         if( ind[1] & 2 )
00344                                 right = Taylor + ind[3] * J;
00345                         else    right = Rec->GetPar(ind[3]);
00346                         if( d == 0 )
00347                         {       if( ind[1] & 4 )
00348                                         trueCase = Taylor + ind[4] * J;
00349                                 else    trueCase = Rec->GetPar(ind[4]);
00350                                 if( ind[1] & 8 )
00351                                         falseCase = Taylor + ind[5] * J;
00352                                 else    falseCase = Rec->GetPar(ind[5]);
00353                         }
00354                         else
00355                         {       if( ind[1] & 4 )
00356                                         trueCase = Taylor + ind[4] * J + d;
00357                                 else    trueCase = &zero;
00358                                 if( ind[1] & 8 )
00359                                         falseCase = Taylor + ind[5] * J + d;
00360                                 else    falseCase = &zero;
00361                         }
00362                         Z[d] = CondExpOp(
00363                                 CompareOp( ind[0] ),
00364                                 *left,
00365                                 *right,
00366                                 *trueCase,
00367                                 *falseCase
00368                         );
00369                         break;
00370                         // ---------------------------------------------------
00371 
00372                         case ComOp:
00373                         CPPAD_ASSERT_UNKNOWN( n_var == 0);
00374                         CPPAD_ASSERT_UNKNOWN( n_ind == 4);
00375                         CPPAD_ASSERT_UNKNOWN( ind[1] > 1 );
00376                         if( d == 0 )
00377                         {       if( ind[1] & 1 )
00378                                         result = true;
00379                                 else    result = false;
00380                                 if( ind[1] & 2 )
00381                                         left = Taylor + ind[2] * J;
00382                                 else    left = Rec->GetPar(ind[2]);
00383                                 if( ind[1] & 4 )
00384                                         right = Taylor + ind[3] * J;
00385                                 else    right = Rec->GetPar(ind[3]);
00386                                 switch( CompareOp( ind[0] ) )
00387                                 {       case CompareLt:
00388                                         compareCount += ( result != 
00389                                         LessThanZero(*left - *right) );
00390                                         break;
00391 
00392                                         case CompareLe:
00393                                         compareCount += ( result !=
00394                                         LessThanOrZero(*left - *right) );
00395                                         break;
00396 
00397                                         case CompareEq:
00398                                         compareCount += ( result != 
00399                                         (*left == *right) );
00400                                         break;
00401 
00402                                         case CompareGe:
00403                                         compareCount += ( result !=
00404                                         GreaterThanOrZero(*left - *right) );
00405                                         break;
00406 
00407                                         case CompareGt:
00408                                         compareCount += ( result != 
00409                                         GreaterThanZero(*left - *right) );
00410                                         break;
00411 
00412                                         case CompareNe:
00413                                         compareCount += ( result != 
00414                                         (*left != *right) );
00415                                         break;
00416 
00417                                         default:
00418                                         CPPAD_ASSERT_UNKNOWN(0);
00419                                 }
00420                         }
00421                         break;
00422                         // ---------------------------------------------------
00423 
00424                         case CosOp:
00425                         CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00426                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00427 
00428                         // cosine and sine must come in pairs
00429                         CPPAD_ASSERT_UNKNOWN( n_var == 2);
00430                         CPPAD_ASSERT_UNKNOWN( (i_var+1) < numvar  );
00431 
00432                         // use W for data stored in variable record
00433                         W = Taylor + (i_var+1) * J;
00434                         X   = Taylor + ind[0] * J;
00435                         ForTrigSinCos(d, W, Z, X);
00436                         break;
00437                         // ---------------------------------------------------
00438 
00439                         case CoshOp:
00440                         CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00441                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00442 
00443                         // hyperbolic cosine and sine must come in pairs
00444                         CPPAD_ASSERT_UNKNOWN( n_var == 2);
00445                         CPPAD_ASSERT_UNKNOWN( (i_var+1) < numvar  );
00446 
00447                         // use W for data stored in variable record
00448                         W = Taylor + (i_var+1) * J;
00449                         X   = Taylor + ind[0] * J;
00450                         ForHypSinCos(d, W, Z, X);
00451                         break;
00452                         // -------------------------------------------------
00453 
00454                         case DisOp:
00455                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00456                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00457                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00458                         if( d == 0 ) 
00459                         {       X   = Taylor + ind[0] * J;
00460                                 Z[0] = ADDiscrete<Base>::Eval(ind[1], X[0]);
00461                         }
00462                         else    Z[d] = Base(0);
00463                         break;
00464                         // -------------------------------------------------
00465 
00466                         case DivvvOp:
00467                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00468                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00469                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00470                         CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00471 
00472                         X = Taylor + ind[0] * J;
00473                         Y = Taylor + ind[1] * J;
00474                         ForDivvvOp(d, Z, X, Y);
00475                         break;
00476                         // -------------------------------------------------
00477 
00478                         case DivpvOp:
00479                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00480                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00481                         CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00482 
00483                         Y = Taylor + ind[1] * J;
00484                         P = Rec->GetPar( ind[0] );
00485                         ForDivpvOp(d, Z, P, Y);
00486                         break;
00487                         // -------------------------------------------------
00488 
00489                         case DivvpOp:
00490                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00491                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00492                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00493 
00494                         P = Rec->GetPar( ind[1] );
00495                         X = Taylor + ind[0] * J;
00496                         ForDivvpOp(d, Z, X, P);
00497                         break;
00498                         // -------------------------------------------------
00499 
00500                         case ExpOp:
00501                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00502                         CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00503                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00504 
00505                         X = Taylor + ind[0] * J;
00506                         ForExpOp(d, Z, X);
00507                         break;
00508                         // -------------------------------------------------
00509 
00510                         case InvOp:
00511                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00512                         CPPAD_ASSERT_UNKNOWN( n_ind == 0 );
00513                         break;
00514                         // -------------------------------------------------
00515 
00516                         case LdpOp:
00517                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00518                         CPPAD_ASSERT_UNKNOWN( n_ind == 3 );
00519                         
00520                         CPPAD_ASSERT_UNKNOWN( ind[0] > 0 );
00521                         CPPAD_ASSERT_UNKNOWN( ind[0] < Rec->NumVecInd() );
00522                         CPPAD_ASSERT_UNKNOWN( VectorInd != CPPAD_NULL );
00523                         CPPAD_ASSERT_UNKNOWN( VectorVar != CPPAD_NULL );
00524 
00525                         if( d == 0 )
00526                         {       i   = ind[1];
00527                                 CPPAD_ASSERT_UNKNOWN( 
00528                                         i < VectorInd[ind[0] - 1] 
00529                                 );
00530                                 CPPAD_ASSERT_UNKNOWN( 
00531                                         i + ind[0] < Rec->NumVecInd() 
00532                                 );
00533 
00534                                 if( VectorVar[ i + ind[0] ] )
00535                                 {       i   = VectorInd[ i + ind[0] ];
00536                                         Rec->ReplaceInd(i_ind + 2, i);
00537                                         CPPAD_ASSERT_UNKNOWN(i > 0 );
00538                                         CPPAD_ASSERT_UNKNOWN( i < i_var );
00539                                         Y     = Taylor + i * J;
00540                                         Z[d]  = Y[d];
00541                                 }
00542                                 else
00543                                 {       i   = VectorInd[ i + ind[0] ];
00544                                         Rec->ReplaceInd(i_ind + 2, 0);
00545                                         Z[d] = *(Rec->GetPar(i));
00546                                         i    = 0;
00547                                 }
00548                         }
00549                         else
00550                         {       i = ind[2];
00551                                 if( i > 0 )
00552                                 {       CPPAD_ASSERT_UNKNOWN( i < i_var );
00553                                         Y     = Taylor + i * J;
00554                                         Z[d]  = Y[d];
00555                                 }
00556                                 else    Z[d]  = Base(0);
00557                         }
00558                         break;
00559                         // -------------------------------------------------
00560 
00561                         case LdvOp:
00562                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00563                         CPPAD_ASSERT_UNKNOWN( n_ind == 3 );
00564                         
00565                         CPPAD_ASSERT_UNKNOWN( ind[0] > 0 );
00566                         CPPAD_ASSERT_UNKNOWN( ind[0] < Rec->NumVecInd() );
00567                         CPPAD_ASSERT_UNKNOWN( VectorInd != CPPAD_NULL );
00568                         CPPAD_ASSERT_UNKNOWN( VectorVar != CPPAD_NULL );
00569 
00570                         if( d == 0 )
00571                         {
00572                                 X   = Taylor + ind[1] * J;
00573                                 i   = Integer( X[0] );
00574                                 len = VectorInd[ ind[0] - 1 ];
00575                                 CPPAD_ASSERT_KNOWN( 
00576                                         i < len,
00577                                         "VecAD index value >= vector length"
00578                                 );
00579                                 CPPAD_ASSERT_UNKNOWN( 
00580                                         i + ind[0] < Rec->NumVecInd() 
00581                                 );
00582 
00583                                 if( VectorVar[ i + ind[0] ] )
00584                                 {       i   = VectorInd[ i + ind[0] ];
00585                                         Rec->ReplaceInd(i_ind + 2, i);
00586                                         CPPAD_ASSERT_UNKNOWN(i > 0 );
00587                                         CPPAD_ASSERT_UNKNOWN( i < i_var );
00588                                         Y     = Taylor + i * J;
00589                                         Z[d]  = Y[d];
00590                                 }
00591                                 else
00592                                 {       i   = VectorInd[ i + ind[0] ];
00593                                         Rec->ReplaceInd(i_ind + 2, 0);
00594                                         Z[d] = *(Rec->GetPar(i));
00595                                         i    = 0;
00596                                 }
00597                         }
00598                         else
00599                         {       i = ind[2];
00600                                 if( i > 0 )
00601                                 {       CPPAD_ASSERT_UNKNOWN( i < i_var );
00602                                         Y     = Taylor + i * J;
00603                                         Z[d]  = Y[d];
00604                                 }
00605                                 else    Z[d]  = Base(0);
00606                         }
00607                         break;
00608                         // -------------------------------------------------
00609 
00610                         case LogOp:
00611                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00612                         CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00613                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00614 
00615                         X = Taylor + ind[0] * J;
00616                         ForLogOp(d, Z, X);
00617                         break;
00618                         // -------------------------------------------------
00619 
00620                         case MulvvOp:
00621                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00622                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00623                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00624                         CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00625 
00626                         X = Taylor + ind[0] * J;
00627                         Y = Taylor + ind[1] * J;
00628                         ForMulvvOp(d, Z, X, Y);
00629                         break;
00630                         // -------------------------------------------------
00631 
00632                         case MulpvOp:
00633                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00634                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00635                         CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00636 
00637                         Y = Taylor + ind[1] * J;
00638                         P = Rec->GetPar( ind[0] );
00639                         ForMulpvOp(d, Z, P, Y);
00640                         break;
00641                         // -------------------------------------------------
00642 
00643                         case MulvpOp:
00644                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00645                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00646                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00647 
00648                         X = Taylor + ind[0] * J;
00649                         P = Rec->GetPar( ind[1] );
00650                         ForMulvpOp(d, Z, X, P);
00651                         break;
00652                         // -------------------------------------------------
00653 
00654                         case NonOp:
00655                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00656                         CPPAD_ASSERT_UNKNOWN( n_ind == 0 );
00657                         break;
00658                         // -------------------------------------------------
00659 
00660                         case ParOp:
00661                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00662                         CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00663 
00664                         P = Rec->GetPar( ind[0] );
00665                         if( d == 0 )
00666                                 Z[d] = *P;
00667                         else    Z[d] = Base(0); 
00668                         break;
00669                         // -------------------------------------------------
00670 
00671                         case PowvpOp:
00672                         CPPAD_ASSERT_UNKNOWN( n_var == 3);
00673                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00674                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var);
00675                         U = Z + J;
00676                         W = U + J;
00677 
00678                         // u = log(x)
00679                         X = Taylor + ind[0] * J;
00680                         ForLogOp(d, U, X);
00681 
00682                         // w = u * y
00683                         Y = Rec->GetPar( ind[1] );
00684                         ForMulvpOp(d, W, U, Y);
00685 
00686                         // z = exp(w)
00687                         // zero order case exactly same as Base type operation
00688                         if( d == 0 )
00689                                 Z[0] = pow(X[0], Y[0]);
00690                         else    ForExpOp(d, Z, W);
00691 
00692                         break;
00693                         // -------------------------------------------------
00694 
00695                         case PowpvOp:
00696                         CPPAD_ASSERT_UNKNOWN( n_var == 3);
00697                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00698                         CPPAD_ASSERT_UNKNOWN( ind[1] < i_var);
00699                         U = Z + J;
00700                         W = U + J;
00701 
00702                         // u = log(x)
00703                         X = Rec->GetPar(ind[0]);
00704                         if( d == 0 )
00705                                 U[d] = log(X[d]);
00706                         else    U[d] = Base(0);
00707 
00708                         // w = u * y
00709                         Y   = Taylor + ind[1] * J;
00710                         ForMulpvOp(d, W, U, Y);
00711 
00712                         // z = exp(w)
00713                         // zero order case exactly same as Base type operation
00714                         if( d == 0 )
00715                                 Z[0] = pow(X[0], Y[0]);
00716                         else    ForExpOp(d, Z, W);
00717 
00718                         break;
00719                         // -------------------------------------------------
00720 
00721                         case PowvvOp:
00722                         CPPAD_ASSERT_UNKNOWN( n_var == 3);
00723                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00724                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var);
00725                         CPPAD_ASSERT_UNKNOWN( ind[1] < i_var);
00726                         U = Z + J;
00727                         W = U + J;
00728 
00729                         // u = log(x)
00730                         X = Taylor + ind[0] * J;
00731                         ForLogOp(d, U, X);
00732 
00733                         // w = u * y
00734                         Y   = Taylor + ind[1] * J;
00735                         ForMulvvOp(d, W, U, Y);
00736 
00737                         // z = exp(w)
00738                         // zero order case exactly same as Base type operation
00739                         if( d == 0 )
00740                                 Z[0] = pow(X[0], Y[0]);
00741                         else    ForExpOp(d, Z, W);
00742 
00743                         break;
00744                         // -------------------------------------------------
00745 
00746                         case PripOp:
00747                         CPPAD_ASSERT_UNKNOWN( n_var == 0 );
00748                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00749                         if( print & (d == 0) )
00750                         {       CPPAD_ASSERT_UNKNOWN( ind[0] < Rec->NumTxt() );
00751                                 std::cout << Rec->GetTxt(ind[0]);
00752                                 std::cout << *(Rec->GetPar(ind[1]));
00753                         }
00754                         break;
00755                         // -------------------------------------------------
00756 
00757                         case PrivOp:
00758                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00759                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00760                         if( print & (d == 0) )
00761                         {       CPPAD_ASSERT_UNKNOWN( ind[0] < Rec->NumTxt() );
00762                                 CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00763 
00764                                 X      = Taylor + ind[1] * J;
00765                                 std::cout << Rec->GetTxt(ind[0]);
00766                                 std::cout << X[0];
00767                         }
00768                         break;
00769                         // -------------------------------------------------
00770 
00771                         case SinOp:
00772                         CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00773                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00774 
00775                         // sine and cosine must come in pairs
00776                         CPPAD_ASSERT_UNKNOWN( n_var == 2);
00777                         CPPAD_ASSERT_UNKNOWN( (i_var+1) < numvar  );
00778 
00779                         // use W for data stored in second variable
00780                         W = Taylor + (i_var+1) * J;
00781                         X   = Taylor + ind[0] * J;
00782                         ForTrigSinCos(d, Z, W, X);
00783                         break;
00784                         // -------------------------------------------------
00785 
00786                         case SinhOp:
00787                         CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00788                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00789 
00790                         // sine and cosine must come in pairs
00791                         CPPAD_ASSERT_UNKNOWN( n_var == 2);
00792                         CPPAD_ASSERT_UNKNOWN( (i_var+1) < numvar  );
00793 
00794                         // use W for data stored in second variable
00795                         W = Taylor + (i_var+1) * J;
00796                         X   = Taylor + ind[0] * J;
00797                         ForHypSinCos(d, Z, W, X);
00798                         break;
00799                         // -------------------------------------------------
00800 
00801                         case SqrtOp:
00802                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00803                         CPPAD_ASSERT_UNKNOWN( n_ind == 1 );
00804                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00805 
00806                         X = Taylor + ind[0] * J;
00807                         ForSqrtOp(d, Z, X);
00808                         break;
00809                         // -------------------------------------------------
00810 
00811                         case StppOp:
00812                         CPPAD_ASSERT_UNKNOWN( n_var == 0);
00813                         CPPAD_ASSERT_UNKNOWN( n_ind == 3 );
00814 
00815                         if( d == 0 )
00816                         {       CPPAD_ASSERT_UNKNOWN( VectorInd != CPPAD_NULL );
00817                                 CPPAD_ASSERT_UNKNOWN( VectorVar != CPPAD_NULL );
00818                                 CPPAD_ASSERT_UNKNOWN( ind[0] < Rec->NumVecInd() );
00819 
00820                                 i   = ind[1];
00821                                 CPPAD_ASSERT_UNKNOWN(i < VectorInd[ind[0] - 1]);
00822                                 CPPAD_ASSERT_UNKNOWN( 
00823                                         i + ind[0] < Rec->NumVecInd() 
00824                                 );
00825                                 VectorInd[ i + ind[0] ] = ind[2];
00826                                 VectorVar[ i + ind[0] ] = false;
00827 
00828                                 Z[d] = *( Rec->GetPar( ind[2] ) );
00829                         }
00830                         break;
00831                         // -------------------------------------------------
00832 
00833                         case StpvOp:
00834                         CPPAD_ASSERT_UNKNOWN( n_var == 0);
00835                         CPPAD_ASSERT_UNKNOWN( n_ind == 3 );
00836 
00837                         if( d == 0 )
00838                         {       CPPAD_ASSERT_UNKNOWN( VectorInd != CPPAD_NULL );
00839                                 CPPAD_ASSERT_UNKNOWN( VectorVar != CPPAD_NULL );
00840                                 CPPAD_ASSERT_UNKNOWN( ind[0] < Rec->NumVecInd() );
00841                                 CPPAD_ASSERT_UNKNOWN( ind[2] < i_var );
00842 
00843                                 i   = ind[1];
00844                                 CPPAD_ASSERT_UNKNOWN(i < VectorInd[ind[0] - 1]);
00845                                 CPPAD_ASSERT_UNKNOWN( 
00846                                         i + ind[0] < Rec->NumVecInd() 
00847                                 );
00848                                 VectorInd[ i + ind[0] ] = ind[2];
00849                                 VectorVar[ i + ind[0] ] = true;
00850                         }
00851                         break;
00852                         // -------------------------------------------------
00853 
00854                         case StvpOp:
00855                         CPPAD_ASSERT_UNKNOWN( n_var == 0);
00856                         CPPAD_ASSERT_UNKNOWN( n_ind == 3 );
00857 
00858                         if( d == 0 )
00859                         {       CPPAD_ASSERT_UNKNOWN( VectorInd != CPPAD_NULL );
00860                                 CPPAD_ASSERT_UNKNOWN( VectorVar != CPPAD_NULL );
00861                                 CPPAD_ASSERT_UNKNOWN( ind[0] < Rec->NumVecInd() );
00862                                 CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00863 
00864                                 X   = Taylor + ind[1] * J;
00865                                 i   = Integer( X[0] );
00866                                 len = VectorInd[ ind[0] - 1 ];
00867                                 CPPAD_ASSERT_KNOWN( 
00868                                         i < len,
00869                                         "VecAD index value >= vector length"
00870                                 );
00871                                 CPPAD_ASSERT_UNKNOWN( 
00872                                         i + ind[0] < Rec->NumVecInd() 
00873                                 );
00874                                 VectorInd[ i + ind[0] ] = ind[2];
00875                                 VectorVar[ i + ind[0] ] = false;
00876 
00877                                 Z[d] = *( Rec->GetPar( ind[2] ) );
00878                         }
00879                         break;
00880                         // -------------------------------------------------
00881 
00882                         case StvvOp:
00883                         CPPAD_ASSERT_UNKNOWN( n_var == 0);
00884                         CPPAD_ASSERT_UNKNOWN( n_ind == 3 );
00885 
00886                         if( d == 0 )
00887                         {       CPPAD_ASSERT_UNKNOWN( VectorInd != CPPAD_NULL );
00888                                 CPPAD_ASSERT_UNKNOWN( VectorVar != CPPAD_NULL );
00889                                 CPPAD_ASSERT_UNKNOWN( ind[0] < Rec->NumVecInd() );
00890                                 CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00891                                 CPPAD_ASSERT_UNKNOWN( ind[2] < i_var );
00892 
00893                                 X   = Taylor + ind[1] * J;
00894                                 i   = Integer( X[0] );
00895                                 len = VectorInd[ ind[0] - 1 ];
00896                                 CPPAD_ASSERT_KNOWN( 
00897                                         i < len,
00898                                         "VecAD index value >= vector length"
00899                                 );
00900                                 CPPAD_ASSERT_UNKNOWN( 
00901                                         i + ind[0] < Rec->NumVecInd() 
00902                                 );
00903                                 VectorInd[ i + ind[0] ] = ind[2];
00904                                 VectorVar[ i + ind[0] ] = true;
00905                         }
00906                         break;
00907                         // -------------------------------------------------
00908 
00909                         case SubvvOp:
00910                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00911                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00912                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00913                         CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00914 
00915                         X = Taylor + ind[0] * J;
00916                         Y = Taylor + ind[1] * J;
00917                         ForSubvvOp(d, Z, X, Y);
00918                         break;
00919                         // -------------------------------------------------
00920 
00921                         case SubpvOp:
00922                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00923                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00924                         CPPAD_ASSERT_UNKNOWN( ind[1] < i_var );
00925 
00926                         P = Rec->GetPar( ind[0] );
00927                         Y = Taylor + ind[1] * J;
00928                         ForSubpvOp(d, Z, P, Y);
00929                         break;
00930                         // -------------------------------------------------
00931 
00932                         case SubvpOp:
00933                         CPPAD_ASSERT_UNKNOWN( n_var == 1);
00934                         CPPAD_ASSERT_UNKNOWN( n_ind == 2 );
00935                         CPPAD_ASSERT_UNKNOWN( ind[0] < i_var );
00936 
00937                         X = Taylor + ind[0] * J;
00938                         P = Rec->GetPar( ind[1] );
00939                         ForSubvpOp(d, Z, X, P);
00940                         break;
00941                         // -------------------------------------------------
00942 
00943                         default:
00944                         CPPAD_ASSERT_UNKNOWN(0);
00945                 }
00946 # if CPPAD_FORWARD_SWEEP_TRACE
00947                 printOp(
00948                         std::cout, 
00949                         Rec,
00950                         i_var,
00951                         op, 
00952                         ind,
00953                         d + 1, 
00954                         Z, 
00955                         0, 
00956                         (Base *) CPPAD_NULL
00957                 );
00958         }
00959         std::cout << std::endl;
00960 # else
00961         }
00962 # endif
00963         CPPAD_ASSERT_UNKNOWN( (i_var + n_var) == Rec->TotNumVar() );
00964         if( VectorInd != CPPAD_NULL )
00965                 CPPAD_TRACK_DEL_VEC(VectorInd);
00966         if( VectorVar != CPPAD_NULL )
00967                 CPPAD_TRACK_DEL_VEC(VectorVar);
00968 
00969         return compareCount;
00970 }
00971 
00972 } // END CppAD namespace
00973 
00974 # undef CPPAD_FORWARD_SWEEP_TRACE
00975 
00976 # endif
/home/coin/svn-release/OptimizationSuite-1.1.0/cppad/cppad/local/forward_sweep.hpp
