CppAD: A C++ Algorithmic Differentiation Package  20130918
rev_hes_sweep.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_REV_HES_SWEEP_INCLUDED
00003 # define CPPAD_REV_HES_SWEEP_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Eclipse Public License Version 1.0.
00011 
00012 A copy of this license is included in the COPYING file of this distribution.
00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00014 -------------------------------------------------------------------------- */
00015 
00016 namespace CppAD { // BEGIN_CPPAD_NAMESPACE
00017 /*!
00018 \file rev_hes_sweep.hpp
00019 Compute Reverse mode Hessian sparsity patterns.
00020 */
00021 
00022 /*!
00023 \def CPPAD_REV_HES_SWEEP_TRACE
00024 This value is either zero or one. 
00025 Zero is the normal operational value.
00026 If it is one, a trace of every rev_hes_sweep computation is printed.
00027 */
00028 # define CPPAD_REV_HES_SWEEP_TRACE 0
00029 
00030 /*!
00031 Given the forward Jacobian sparsity pattern for all the variables,
00032 and the reverse Jacobian sparsity pattern for the dependent variables,
00033 RevHesSweep computes the Hessian sparsity pattern for all the independent 
00034 variables.
00035 
00036 \tparam Base
00037 base type for the operator; i.e., this operation sequence was recorded
00038 using AD< \a Base > and computations by this routine are done using type 
00039 \a Base.
00040 
00041 \tparam Vector_set
00042 is the type used for vectors of sets. It can be either
00043 \c sparse_pack, \c sparse_set, or \c sparse_list.
00044 
00045 \param n
00046 is the number of independent variables on the tape.
00047 
00048 \param numvar
00049 is the total number of variables on the tape; i.e.,
00050 \a play->num_var_rec().
00051 This is also the number of rows in the entire sparsity pattern 
00052 \a rev_hes_sparse.
00053 
00054 \param play
00055 The information stored in \a play
00056 is a recording of the operations corresponding to a function
00057 \f[
00058      F : {\bf R}^n \rightarrow {\bf R}^m
00059 \f]
00060 where \f$ n \f$ is the number of independent variables
00061 and \f$ m \f$ is the number of dependent variables.
00062 The object \a play is effectly constant.
00063 It is not declared const because while playing back the tape
00064 the object \a play holds information about the currentl location
00065 with in the tape and this changes during playback.
00066 
00067 \param for_jac_sparse
00068 For i = 0 , ... , \a numvar - 1, 
00069 (for all the variables on the tape),
00070 the forward Jacobian sparsity pattern for the variable with index i
00071 corresponds to the set with index i in \a for_jac_sparse.
00072 
00073 \param RevJac
00074 \b Input:
00075 For i = 0, ... , \a numvar - 1 
00076 the if the variable with index i on the tape is an dependent variable and
00077 included in the Hessian, \a RevJac[ i ] is equal to true,
00078 otherwise it is equal to false. 
00079 \n
00080 \n
00081 \b Output: The values in \a RevJac upon return are not specified; i.e.,
00082 it is used for temporary work space.
00083 
00084 \param rev_hes_sparse
00085 The reverse Hessian sparsity pattern for the variable with index i
00086 corresponds to the set with index i in \a rev_hes_sparse.
00087 \n
00088 \n
00089 \b Input: For i = 0 , ... , \a numvar - 1  
00090 the reverse Hessian sparsity pattern for the variable with index i is empty.
00091 \n
00092 \n
00093 \b Output: For j = 1 , ... , \a n,
00094 the reverse Hessian sparsity pattern for the independent dependent variable 
00095 with index (j-1) is given by the set with index j
00096 in \a rev_hes_sparse. 
00097 The values in the rest of \a rev_hes_sparse are not specified; i.e.,
00098 they are used for temporary work space.
00099 */
00100 
00101 template <class Base, class Vector_set>
00102 void RevHesSweep(
00103      size_t                n,
00104      size_t                numvar,
00105      player<Base>         *play,
00106      Vector_set&           for_jac_sparse, // should be const
00107      bool*                 RevJac,
00108      Vector_set&           rev_hes_sparse
00109 )
00110 {
00111      OpCode           op;
00112      size_t         i_op;
00113      size_t        i_var;
00114 
00115      const addr_t*   arg = CPPAD_NULL;
00116 
00117      // length of the parameter vector (used by CppAD assert macros)
00118      const size_t num_par = play->num_par_rec();
00119 
00120      size_t             i, j, k;
00121 
00122      // check numvar argument
00123      CPPAD_ASSERT_UNKNOWN( play->num_var_rec()     == numvar );
00124      CPPAD_ASSERT_UNKNOWN( for_jac_sparse.n_set() == numvar );
00125      CPPAD_ASSERT_UNKNOWN( rev_hes_sparse.n_set() == numvar );
00126      CPPAD_ASSERT_UNKNOWN( numvar > 0 );
00127 
00128      // upper limit exclusive for set elements
00129      size_t limit   = rev_hes_sparse.end();
00130      CPPAD_ASSERT_UNKNOWN( for_jac_sparse.end() == limit );
00131 
00132      // check number of sets match
00133      CPPAD_ASSERT_UNKNOWN( 
00134           for_jac_sparse.n_set() == rev_hes_sparse.n_set()
00135      );
00136 
00137      // vecad_sparsity contains a sparsity pattern for each VecAD object.
00138      // vecad_ind maps a VecAD index (beginning of the VecAD object) 
00139      // to the index for the corresponding set in vecad_sparsity.
00140      size_t num_vecad_ind   = play->num_vec_ind_rec();
00141      size_t num_vecad_vec   = play->num_vecad_vec_rec();
00142      Vector_set vecad_sparse;
00143      vecad_sparse.resize(num_vecad_vec, limit);
00144      pod_vector<size_t> vecad_ind;
00145      pod_vector<bool>   vecad_jac;
00146      if( num_vecad_vec > 0 )
00147      {    size_t length;
00148           vecad_ind.extend(num_vecad_ind);
00149           vecad_jac.extend(num_vecad_vec);
00150           j             = 0;
00151           for(i = 0; i < num_vecad_vec; i++)
00152           {    // length of this VecAD
00153                length   = play->GetVecInd(j);
00154                // set vecad_ind to proper index for this VecAD
00155                vecad_ind[j] = i; 
00156                // make all other values for this vector invalid
00157                for(k = 1; k <= length; k++)
00158                     vecad_ind[j+k] = num_vecad_vec;
00159                // start of next VecAD
00160                j       += length + 1;
00161                // initialize this vector's reverse jacobian value 
00162                vecad_jac[i] = false;
00163           }
00164           CPPAD_ASSERT_UNKNOWN( j == play->num_vec_ind_rec() );
00165      }
00166 
00167      // work space used by UserOp.
00168      vector<size_t>     user_ix;  // variable indices for argument vector x
00169      typedef std::set<size_t> size_set;
00170      size_set::iterator set_itr;  // iterator for a standard set
00171      size_set::iterator set_end;  // end of iterator sequence
00172      vector< size_set > set_r;    // forward Jacobian sparsity for x
00173      vector< size_set > set_u;    // reverse Hessian sparsity for y
00174      vector< size_set > set_v;    // reverse Hessian sparsity for x
00175      //
00176      vector<bool>       bool_r;   // bool forward Jacobian sparsity for x
00177      vector<bool>       bool_u;   // bool reverse Hessian sparsity for y
00178      vector<bool>       bool_v;   // bool reverse Hessian sparsity for x
00179      //
00180      vector<bool>       user_vx;  // which components of x are variables
00181      vector<bool>       user_s;   // reverse Jacobian sparsity for y
00182      vector<bool>       user_t;   // reverse Jacobian sparsity for x
00183      const size_t user_q = limit; // maximum element plus one
00184      size_t user_index = 0;       // indentifier for this atomic operation
00185      size_t user_id    = 0;       // user identifier for this call to operator
00186      size_t user_i     = 0;       // index in result vector
00187      size_t user_j     = 0;       // index in argument vector
00188      size_t user_m     = 0;       // size of result vector
00189      size_t user_n     = 0;       // size of arugment vector
00190      //
00191      atomic_base<Base>* user_atom = CPPAD_NULL; // user's atomic op calculator
00192      bool               user_bool = false;      // use bool or set sparsity ?
00193 # ifndef NDEBUG
00194      bool               user_ok   = false;      // atomic op return value
00195 # endif
00196      // next expected operator in a UserOp sequence
00197      enum { user_start, user_arg, user_ret, user_end } user_state = user_end;
00198 
00199 
00200      // Initialize
00201      play->reverse_start(op, arg, i_op, i_var);
00202      CPPAD_ASSERT_UNKNOWN( op == EndOp );
00203 # if CPPAD_REV_HES_SWEEP_TRACE
00204      std::cout << std::endl;
00205      CppAD::vectorBool zf_value(limit);
00206      CppAD::vectorBool zh_value(limit);
00207 # endif
00208      bool more_operators = true;
00209      while(more_operators)
00210      {
00211           // next op
00212           play->reverse_next(op, arg, i_op, i_var);
00213 # ifndef NDEBUG
00214           if( i_op <= n )
00215           {    CPPAD_ASSERT_UNKNOWN((op == InvOp) | (op == BeginOp));
00216           }
00217           else CPPAD_ASSERT_UNKNOWN((op != InvOp) & (op != BeginOp));
00218 # endif
00219 
00220           // rest of information depends on the case
00221           switch( op )
00222           {
00223                case AbsOp:
00224                CPPAD_ASSERT_NARG_NRES(op, 1, 1)
00225                reverse_sparse_hessian_linear_unary_op(
00226                i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
00227                );
00228                break;
00229                // -------------------------------------------------
00230 
00231                case AddvvOp:
00232                CPPAD_ASSERT_NARG_NRES(op, 2, 1)
00233                reverse_sparse_hessian_addsub_op(
00234                i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse
00235                );
00236                break;
00237                // -------------------------------------------------
00238 
00239                case AddpvOp:
00240                CPPAD_ASSERT_NARG_NRES(op, 2, 1)
00241                reverse_sparse_hessian_linear_unary_op(
00242                i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse
00243                );
00244                break;
00245                // -------------------------------------------------
00246 
00247                case AcosOp:
00248                // sqrt(1 - x * x), acos(x)
00249                CPPAD_ASSERT_NARG_NRES(op, 1, 2)
00250                reverse_sparse_hessian_nonlinear_unary_op(
00251                i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
00252                );
00253                break;
00254                // -------------------------------------------------
00255 
00256                case AsinOp:
00257                // sqrt(1 - x * x), asin(x)
00258                CPPAD_ASSERT_NARG_NRES(op, 1, 2)
00259                reverse_sparse_hessian_nonlinear_unary_op(
00260                i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
00261                );
00262                break;
00263                // -------------------------------------------------
00264 
00265                case AtanOp:
00266                // 1 + x * x, atan(x)
00267                CPPAD_ASSERT_NARG_NRES(op, 1, 2)
00268                reverse_sparse_hessian_nonlinear_unary_op(
00269                i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
00270                );
00271                break;
00272                // -------------------------------------------------
00273 
00274                case BeginOp:
00275                CPPAD_ASSERT_NARG_NRES(op, 1, 1)
00276                more_operators = false;
00277                break;
00278                // -------------------------------------------------
00279 
00280                case CSkipOp:
00281                // CSkipOp has a variable number of arguments and
00282                // reverse_next thinks it one has one argument.
00283                // We must inform reverse_next of this special case.
00284                play->reverse_cskip(op, arg, i_op, i_var);
00285                break;
00286                // -------------------------------------------------
00287 
00288                case CSumOp:
00289                // CSumOp has a variable number of arguments and
00290                // reverse_next thinks it one has one argument.
00291                // We must inform reverse_next of this special case.
00292                play->reverse_csum(op, arg, i_op, i_var);
00293                reverse_sparse_hessian_csum_op(
00294                     i_var, arg, RevJac, rev_hes_sparse
00295                );
00296                break;
00297                // -------------------------------------------------
00298 
00299                case CExpOp:
00300                reverse_sparse_hessian_cond_op(
00301                     i_var, arg, num_par, RevJac, rev_hes_sparse
00302                );
00303                break;
00304                // ---------------------------------------------------
00305 
00306                case ComOp:
00307                CPPAD_ASSERT_NARG_NRES(op, 4, 0)
00308                CPPAD_ASSERT_UNKNOWN( arg[1] > 1 );
00309                break;
00310                // --------------------------------------------------
00311 
00312                case CosOp:
00313                // sin(x), cos(x)
00314                CPPAD_ASSERT_NARG_NRES(op, 1, 2)
00315                reverse_sparse_hessian_nonlinear_unary_op(
00316                i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
00317                );
00318                break;
00319                // ---------------------------------------------------
00320 
00321                case CoshOp:
00322                // sinh(x), cosh(x)
00323                CPPAD_ASSERT_NARG_NRES(op, 1, 2)
00324                reverse_sparse_hessian_nonlinear_unary_op(
00325                i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
00326                );
00327                break;
00328                // -------------------------------------------------
00329 
00330                case DisOp:
00331                CPPAD_ASSERT_NARG_NRES(op, 2, 1)
00332                // derivativve is identically zero
00333                break;
00334                // -------------------------------------------------
00335 
00336                case DivvvOp:
00337                CPPAD_ASSERT_NARG_NRES(op, 2, 1)
00338                reverse_sparse_hessian_div_op(
00339                i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse
00340                );
00341                break;
00342                // -------------------------------------------------
00343 
00344                case DivpvOp:
00345                CPPAD_ASSERT_NARG_NRES(op, 2, 1)
00346                reverse_sparse_hessian_nonlinear_unary_op(
00347                i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse
00348                );
00349                break;
00350                // -------------------------------------------------
00351 
00352                case DivvpOp:
00353                CPPAD_ASSERT_NARG_NRES(op, 2, 1)
00354                reverse_sparse_hessian_linear_unary_op(
00355                i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
00356                );
00357                break;
00358                // -------------------------------------------------
00359 
00360                case ExpOp:
00361                CPPAD_ASSERT_NARG_NRES(op, 1, 1)
00362                reverse_sparse_hessian_nonlinear_unary_op(
00363                i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
00364                );
00365                break;
00366                // -------------------------------------------------
00367 
00368                case InvOp:
00369                CPPAD_ASSERT_NARG_NRES(op, 0, 1)
00370                // Z is already defined
00371                break;
00372                // -------------------------------------------------
00373 
00374                case LdpOp:
00375                reverse_sparse_hessian_load_op(
00376                     op,
00377                     i_var,
00378                     arg,
00379                     num_vecad_ind,
00380                     vecad_ind.data(),
00381                     rev_hes_sparse,
00382                     vecad_sparse,
00383                     RevJac,
00384                     vecad_jac.data()
00385                );
00386                break;
00387                // -------------------------------------------------
00388 
00389                case LdvOp:
00390                reverse_sparse_hessian_load_op(
00391                     op,
00392                     i_var,
00393                     arg,
00394                     num_vecad_ind,
00395                     vecad_ind.data(),
00396                     rev_hes_sparse,
00397                     vecad_sparse,
00398                     RevJac,
00399                     vecad_jac.data()
00400                );
00401                break;
00402                // -------------------------------------------------
00403 
00404                case LogOp:
00405                CPPAD_ASSERT_NARG_NRES(op, 1, 1)
00406                reverse_sparse_hessian_nonlinear_unary_op(
00407                i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
00408                );
00409                break;
00410                // -------------------------------------------------
00411 
00412                case MulvvOp:
00413                CPPAD_ASSERT_NARG_NRES(op, 2, 1)
00414                reverse_sparse_hessian_mul_op(
00415                i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse
00416                );
00417                break;
00418                // -------------------------------------------------
00419 
00420                case MulpvOp:
00421                CPPAD_ASSERT_NARG_NRES(op, 2, 1)
00422                reverse_sparse_hessian_linear_unary_op(
00423                i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse
00424                );
00425                break;
00426                // -------------------------------------------------
00427 
00428                case ParOp:
00429                CPPAD_ASSERT_NARG_NRES(op, 1, 1)
00430 
00431                break;
00432                // -------------------------------------------------
00433 
00434                case PowpvOp:
00435                CPPAD_ASSERT_NARG_NRES(op, 2, 3)
00436                reverse_sparse_hessian_nonlinear_unary_op(
00437                i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse
00438                );
00439                break;
00440                // -------------------------------------------------
00441 
00442                case PowvpOp:
00443                CPPAD_ASSERT_NARG_NRES(op, 2, 3)
00444                reverse_sparse_hessian_nonlinear_unary_op(
00445                i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
00446                );
00447                break;
00448                // -------------------------------------------------
00449 
00450                case PowvvOp:
00451                CPPAD_ASSERT_NARG_NRES(op, 2, 3)
00452                         reverse_sparse_hessian_pow_op(
00453                i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse
00454                );
00455                break;
00456                // -------------------------------------------------
00457 
00458                case PriOp:
00459                CPPAD_ASSERT_NARG_NRES(op, 5, 0);
00460                break;
00461                // -------------------------------------------------
00462 
00463                case SignOp:
00464                CPPAD_ASSERT_NARG_NRES(op, 1, 1);
00465                // Derivative is identiaclly zero
00466                break;
00467                // -------------------------------------------------
00468 
00469                case SinOp:
00470                // cos(x), sin(x)
00471                CPPAD_ASSERT_NARG_NRES(op, 1, 2)
00472                reverse_sparse_hessian_nonlinear_unary_op(
00473                i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
00474                );
00475                break;
00476                // -------------------------------------------------
00477 
00478                case SinhOp:
00479                // cosh(x), sinh(x)
00480                CPPAD_ASSERT_NARG_NRES(op, 1, 2)
00481                reverse_sparse_hessian_nonlinear_unary_op(
00482                i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
00483                );
00484                break;
00485                // -------------------------------------------------
00486 
00487                case SqrtOp:
00488                CPPAD_ASSERT_NARG_NRES(op, 1, 1)
00489                reverse_sparse_hessian_nonlinear_unary_op(
00490                i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
00491                );
00492                break;
00493                // -------------------------------------------------
00494 
00495                case StppOp:
00496                // sparsity cannot propagate through a parameter
00497                CPPAD_ASSERT_NARG_NRES(op, 3, 0)
00498                break;
00499                // -------------------------------------------------
00500 
00501                case StpvOp:
00502                reverse_sparse_hessian_store_op(
00503                     op,
00504                     arg,
00505                     num_vecad_ind,
00506                     vecad_ind.data(),
00507                     rev_hes_sparse,
00508                     vecad_sparse,
00509                     RevJac,
00510                     vecad_jac.data()
00511                );
00512                break;
00513                // -------------------------------------------------
00514 
00515                case StvpOp:
00516                // sparsity cannot propagate through a parameter
00517                CPPAD_ASSERT_NARG_NRES(op, 3, 0)
00518                break;
00519                // -------------------------------------------------
00520 
00521                case StvvOp:
00522                reverse_sparse_hessian_store_op(
00523                     op,
00524                     arg,
00525                     num_vecad_ind,
00526                     vecad_ind.data(),
00527                     rev_hes_sparse,
00528                     vecad_sparse,
00529                     RevJac,
00530                     vecad_jac.data()
00531                );
00532                break;
00533                // -------------------------------------------------
00534 
00535                case SubvvOp:
00536                CPPAD_ASSERT_NARG_NRES(op, 2, 1)
00537                reverse_sparse_hessian_addsub_op(
00538                i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse
00539                );
00540                break;
00541                // -------------------------------------------------
00542 
00543                case SubpvOp:
00544                CPPAD_ASSERT_NARG_NRES(op, 2, 1)
00545                reverse_sparse_hessian_linear_unary_op(
00546                i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse
00547                );
00548                break;
00549                // -------------------------------------------------
00550 
00551                case SubvpOp:
00552                CPPAD_ASSERT_NARG_NRES(op, 2, 1)
00553                reverse_sparse_hessian_linear_unary_op(
00554                i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
00555                );
00556                break;
00557                // -------------------------------------------------
00558 
00559                case TanOp:
00560                // tan(x)^2, tan(x)
00561                CPPAD_ASSERT_NARG_NRES(op, 1, 2)
00562                reverse_sparse_hessian_nonlinear_unary_op(
00563                i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
00564                );
00565                break;
00566                // -------------------------------------------------
00567 
00568                case TanhOp:
00569                // tanh(x)^2, tanh(x)
00570                CPPAD_ASSERT_NARG_NRES(op, 1, 2)
00571                reverse_sparse_hessian_nonlinear_unary_op(
00572                i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
00573                );
00574                break;
00575                // -------------------------------------------------
00576 
00577                case UserOp:
00578                // start or end an atomic operation sequence
00579                CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 );
00580                CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 );
00581                if( user_state == user_end )
00582                {    user_index = arg[0];
00583                     user_id    = arg[1];
00584                     user_n     = arg[2];
00585                     user_m     = arg[3];
00586                     user_atom  = atomic_base<Base>::class_object(user_index);
00587 # ifndef NDEBUG
00588                     if( user_atom == CPPAD_NULL )
00589                     {    std::string msg = 
00590                               atomic_base<Base>::class_name(user_index)
00591                               + ": atomic_base function has been deleted";
00592                          CPPAD_ASSERT_KNOWN(false, msg.c_str() );
00593                     }
00594 # endif
00595                     user_bool  = user_atom->sparsity() ==
00596                                    atomic_base<Base>::bool_sparsity_enum;
00597                     user_ix.resize(user_n);
00598                     user_vx.resize(user_n);
00599                     user_s.resize(user_m);
00600                     user_t.resize(user_n);
00601 
00602                     // simpler to initialize all sparsity patterns as empty
00603                     for(i = 0; i < user_m; i++)
00604                          user_s[i] = false;
00605                     for(i = 0; i < user_n; i++)
00606                          user_t[i] = false;
00607                     if( user_bool )
00608                     {    bool_r.resize(user_n * user_q);
00609                          bool_u.resize(user_m * user_q);
00610                          bool_v.resize(user_n * user_q);
00611                          // simpler to initialize all patterns as empty
00612                          for(i = 0; i < user_m; i++)
00613                          {
00614                               for(j = 0; j < user_q; j++)
00615                                    bool_u[ i * user_q + j] = false;
00616                          }
00617                          for(i = 0; i < user_n; i++)
00618                          {
00619                               for(j = 0; j < user_q; j++)
00620                               {    bool_r[ i * user_q + j] = false;
00621                                    bool_v[ i * user_q + j] = false;
00622                               }
00623                          }
00624                     }
00625                     else
00626                     {    set_r.resize(user_n);
00627                          set_u.resize(user_m);
00628                          set_v.resize(user_n);
00629                          for(i = 0; i < user_m; i++)
00630                               set_u[i].clear();
00631                          for(i = 0; i < user_n; i++)
00632                          {    set_r[i].clear();
00633                               set_v[i].clear();
00634                          }
00635                     }
00636                     user_j     = user_n;
00637                     user_i     = user_m;
00638                     user_state = user_ret;
00639                }
00640                else
00641                {    CPPAD_ASSERT_UNKNOWN( user_state == user_start );
00642                     CPPAD_ASSERT_UNKNOWN( user_index == size_t(arg[0]) );
00643                     CPPAD_ASSERT_UNKNOWN( user_id    == size_t(arg[1]) );
00644                     CPPAD_ASSERT_UNKNOWN( user_n     == size_t(arg[2]) );
00645                     CPPAD_ASSERT_UNKNOWN( user_m     == size_t(arg[3]) );
00646                     user_state = user_end;
00647 
00648                     // call users function for this operation
00649                     user_atom->set_id(user_id);
00650 # ifdef NDEBUG
00651                     if( user_bool )
00652                          user_atom->rev_sparse_hes(user_vx,
00653                               user_s, user_t, user_q, bool_r, bool_u, bool_v
00654                     );
00655                     else
00656                          user_atom->rev_sparse_hes(user_vx,
00657                               user_s, user_t, user_q, set_r, set_u, set_v
00658                     );
00659 # else
00660                     if( user_bool )
00661                          user_ok = user_atom->rev_sparse_hes(user_vx,
00662                               user_s, user_t, user_q, bool_r, bool_u, bool_v
00663                     );
00664                     else
00665                          user_ok = user_atom->rev_sparse_hes(user_vx,
00666                               user_s, user_t, user_q, set_r, set_u, set_v
00667                     );
00668                     if( ! user_ok )
00669                     {    std::string msg = 
00670                               atomic_base<Base>::class_name(user_index)
00671                               + ": atomic_base.rev_sparse_hes: returned false";
00672                          CPPAD_ASSERT_KNOWN(false, msg.c_str() );
00673                     }
00674 # endif
00675                     for(i = 0; i < user_n; i++) if( user_ix[i] > 0 )
00676                     {
00677                          size_t  i_x = user_ix[i];
00678                          if( user_t[i] )
00679                               RevJac[i_x] = true;
00680                          if( user_bool )
00681                          {
00682                               for(j = 0; j < user_q; j++)
00683                                    if( bool_v[ i * user_q + j ] )
00684                                         rev_hes_sparse.add_element(i_x, j);
00685                          }
00686                          else
00687                          {
00688                               set_itr = set_v[i].begin();
00689                               set_end = set_v[i].end();
00690                               while( set_itr != set_end )
00691                                    rev_hes_sparse.add_element(i_x, *set_itr++);
00692                          }
00693                     }
00694                }
00695                break;
00696 
00697                case UsrapOp:
00698                // parameter argument in an atomic operation sequence
00699                CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
00700                CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
00701                CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
00702                CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
00703                --user_j;
00704                user_ix[user_j] = 0;
00705                user_vx[user_j] = false;
00706                if( user_j == 0 )
00707                     user_state = user_start;
00708                break;
00709 
00710                case UsravOp:
00711                // variable argument in an atomic operation sequence
00712                CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
00713                CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
00714                CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
00715                CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) <= i_var );
00716                CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
00717                --user_j;
00718                user_ix[user_j] = arg[0];
00719                user_vx[user_j] = true;
00720                for_jac_sparse.begin(arg[0]);
00721                i = for_jac_sparse.next_element();
00722                while( i < user_q )
00723                {    if( user_bool )
00724                          bool_r[ user_j * user_q + i ] = true;
00725                     else
00726                          set_r[user_j].insert(i);
00727                     i = for_jac_sparse.next_element();
00728                }
00729                if( user_j == 0 )
00730                     user_state = user_start;
00731                break;
00732 
00733                case UsrrpOp:
00734                // parameter result in an atomic operation sequence
00735                CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
00736                CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
00737                CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
00738                CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
00739                --user_i;
00740                if( user_i == 0 )
00741                     user_state = user_arg;
00742                break;
00743 
00744                case UsrrvOp:
00745                // variable result in an atomic operation sequence
00746                CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
00747                CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
00748                --user_i;
00749                if( RevJac[i_var] )
00750                {
00751                     user_s[user_i] = true;
00752                }
00753                rev_hes_sparse.begin(i_var);
00754                j = rev_hes_sparse.next_element();
00755                while( j < user_q )
00756                {    if( user_bool )
00757                          bool_u[user_i * user_q + j] = true;
00758                     else
00759                          set_u[user_i].insert(j);
00760                     j = rev_hes_sparse.next_element();
00761                }
00762                if( user_i == 0 )
00763                     user_state = user_arg;
00764                break;
00765 
00766                // -------------------------------------------------
00767 
00768                default:
00769                CPPAD_ASSERT_UNKNOWN(0);
00770           }
00771 # if CPPAD_REV_HES_SWEEP_TRACE
00772           for(j = 0; j < limit; j++)
00773           {    zf_value[j] = false;
00774                zh_value[j] = false;
00775           }
00776           for_jac_sparse.begin(i_var);;
00777           j = for_jac_sparse.next_element();;
00778           while( j < limit )
00779           {    zf_value[j] = true;
00780                j = for_jac_sparse.next_element();
00781           }
00782           rev_hes_sparse.begin(i_var);;
00783           j = rev_hes_sparse.next_element();;
00784           while( j < limit )
00785           {    zh_value[j] = true;
00786                j = rev_hes_sparse.next_element();
00787           }
00788           // should also print RevJac[i_var], but printOp does not
00789           // yet allow for this.
00790           printOp(
00791                std::cout, 
00792                play,
00793                i_op,
00794                i_var,
00795                op, 
00796                arg,
00797                1, 
00798                &zf_value, 
00799                1, 
00800                &zh_value
00801           );
00802 # endif
00803      }
00804      // values corresponding to BeginOp
00805      CPPAD_ASSERT_UNKNOWN( i_op == 0 );
00806      CPPAD_ASSERT_UNKNOWN( i_var == 0 );
00807 
00808      return;
00809 }
00810 } // END_CPPAD_NAMESPACE
00811 
00812 // preprocessor symbols that are local to this file
00813 # undef CPPAD_REV_HES_SWEEP_TRACE
00814 
00815 # endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines