source: branches/opt_cond_exp/cppad/local/rev_hes_sweep.hpp @ 2979

Last change on this file since 2979 was 2979, checked in by bradbell, 7 years ago
  1. Add operator index to trace output.
  2. Check that new operator values have been set and are valid.
  3. Reinitialize cskip_op_ after optimization.

optimize.hpp: Add test of entirely removing an atomic call.

  • Property svn:keywords set to Id
File size: 22.5 KB
Line 
1/* $Id: rev_hes_sweep.hpp 2979 2013-10-20 18:00:28Z bradbell $ */
2# ifndef CPPAD_REV_HES_SWEEP_INCLUDED
3# define CPPAD_REV_HES_SWEEP_INCLUDED
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
7
8CppAD is distributed under multiple licenses. This distribution is under
9the terms of the
10                    Eclipse Public License Version 1.0.
11
12A copy of this license is included in the COPYING file of this distribution.
13Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
14-------------------------------------------------------------------------- */
15
16namespace CppAD { // BEGIN_CPPAD_NAMESPACE
17/*!
18\defgroup rev_hes_sweep_hpp rev_hes_sweep.hpp
19\{
20\file rev_hes_sweep.hpp
21Compute Reverse mode Hessian sparsity patterns.
22*/
23
24/*!
25\def CPPAD_REV_HES_SWEEP_TRACE
26This value is either zero or one.
27Zero is the normal operational value.
28If it is one, a trace of every rev_hes_sweep computation is printed.
29*/
30# define CPPAD_REV_HES_SWEEP_TRACE 0
31
32/*!
33Given the forward Jacobian sparsity pattern for all the variables,
34and the reverse Jacobian sparsity pattern for the dependent variables,
35RevHesSweep computes the Hessian sparsity pattern for all the independent
36variables.
37
38\tparam Base
39base type for the operator; i.e., this operation sequence was recorded
40using AD< \a Base > and computations by this routine are done using type
41\a Base.
42
43\tparam Vector_set
44is the type used for vectors of sets. It can be either
45\c sparse_pack, \c sparse_set, or \c sparse_list.
46
47\param n
48is the number of independent variables on the tape.
49
50\param numvar
51is the total number of variables on the tape; i.e.,
52\a play->num_rec_var().
53This is also the number of rows in the entire sparsity pattern
54\a rev_hes_sparse.
55
56\param play
57The information stored in \a play
58is a recording of the operations corresponding to a function
59\f[
60        F : {\bf R}^n \rightarrow {\bf R}^m
61\f]
62where \f$ n \f$ is the number of independent variables
63and \f$ m \f$ is the number of dependent variables.
64The object \a play is effectly constant.
65It is not declared const because while playing back the tape
66the object \a play holds information about the currentl location
67with in the tape and this changes during playback.
68
69\param for_jac_sparse
70For i = 0 , ... , \a numvar - 1,
71(for all the variables on the tape),
72the forward Jacobian sparsity pattern for the variable with index i
73corresponds to the set with index i in \a for_jac_sparse.
74
75\param RevJac
76\b Input:
77For i = 0, ... , \a numvar - 1
78the if the variable with index i on the tape is an dependent variable and
79included in the Hessian, \a RevJac[ i ] is equal to true,
80otherwise it is equal to false.
81\n
82\n
83\b Output: The values in \a RevJac upon return are not specified; i.e.,
84it is used for temporary work space.
85
86\param rev_hes_sparse
87The reverse Hessian sparsity pattern for the variable with index i
88corresponds to the set with index i in \a rev_hes_sparse.
89\n
90\n
91\b Input: For i = 0 , ... , \a numvar - 1 
92the reverse Hessian sparsity pattern for the variable with index i is empty.
93\n
94\n
95\b Output: For j = 1 , ... , \a n,
96the reverse Hessian sparsity pattern for the independent dependent variable
97with index (j-1) is given by the set with index j
98in \a rev_hes_sparse.
99The values in the rest of \a rev_hes_sparse are not specified; i.e.,
100they are used for temporary work space.
101*/
102
103template <class Base, class Vector_set>
104void RevHesSweep(
105        size_t                n,
106        size_t                numvar,
107        player<Base>         *play,
108        Vector_set&           for_jac_sparse, // should be const
109        bool*                 RevJac,
110        Vector_set&           rev_hes_sparse
111)
112{
113        OpCode           op;
114        size_t         i_op;
115        size_t        i_var;
116
117        const addr_t*   arg = CPPAD_NULL;
118
119        // length of the parameter vector (used by CppAD assert macros)
120        const size_t num_par = play->num_rec_par();
121
122        size_t             i, j, k;
123
124        // check numvar argument
125        CPPAD_ASSERT_UNKNOWN( play->num_rec_var()     == numvar );
126        CPPAD_ASSERT_UNKNOWN( for_jac_sparse.n_set() == numvar );
127        CPPAD_ASSERT_UNKNOWN( rev_hes_sparse.n_set() == numvar );
128        CPPAD_ASSERT_UNKNOWN( numvar > 0 );
129
130        // upper limit exclusive for set elements
131        size_t limit   = rev_hes_sparse.end();
132        CPPAD_ASSERT_UNKNOWN( for_jac_sparse.end() == limit );
133
134        // check number of sets match
135        CPPAD_ASSERT_UNKNOWN( 
136                for_jac_sparse.n_set() == rev_hes_sparse.n_set()
137        );
138
139        // vecad_sparsity contains a sparsity pattern for each VecAD object.
140        // vecad_ind maps a VecAD index (beginning of the VecAD object)
141        // to the index for the corresponding set in vecad_sparsity.
142        size_t num_vecad_ind   = play->num_rec_vecad_ind();
143        size_t num_vecad_vec   = play->num_rec_vecad_vec();
144        Vector_set vecad_sparse;
145        vecad_sparse.resize(num_vecad_vec, limit);
146        pod_vector<size_t> vecad_ind;
147        pod_vector<bool>   vecad_jac;
148        if( num_vecad_vec > 0 )
149        {       size_t length;
150                vecad_ind.extend(num_vecad_ind);
151                vecad_jac.extend(num_vecad_vec);
152                j             = 0;
153                for(i = 0; i < num_vecad_vec; i++)
154                {       // length of this VecAD
155                        length   = play->GetVecInd(j);
156                        // set vecad_ind to proper index for this VecAD
157                        vecad_ind[j] = i; 
158                        // make all other values for this vector invalid
159                        for(k = 1; k <= length; k++)
160                                vecad_ind[j+k] = num_vecad_vec;
161                        // start of next VecAD
162                        j       += length + 1;
163                        // initialize this vector's reverse jacobian value
164                        vecad_jac[i] = false;
165                }
166                CPPAD_ASSERT_UNKNOWN( j == play->num_rec_vecad_ind() );
167        }
168
169        // work space used by UserOp.
170        vector<size_t>     user_ix;  // variable indices for argument vector x
171        typedef std::set<size_t> size_set;
172        size_set::iterator set_itr;  // iterator for a standard set
173        size_set::iterator set_end;  // end of iterator sequence
174        vector< size_set > set_r;    // forward Jacobian sparsity for x
175        vector< size_set > set_u;    // reverse Hessian sparsity for y
176        vector< size_set > set_v;    // reverse Hessian sparsity for x
177        //
178        vector<bool>       bool_r;   // bool forward Jacobian sparsity for x
179        vector<bool>       bool_u;   // bool reverse Hessian sparsity for y
180        vector<bool>       bool_v;   // bool reverse Hessian sparsity for x
181        //
182        vector<bool>       user_vx;  // which components of x are variables
183        vector<bool>       user_s;   // reverse Jacobian sparsity for y
184        vector<bool>       user_t;   // reverse Jacobian sparsity for x
185        const size_t user_q = limit; // maximum element plus one
186        size_t user_index = 0;       // indentifier for this atomic operation
187        size_t user_id    = 0;       // user identifier for this call to operator
188        size_t user_i     = 0;       // index in result vector
189        size_t user_j     = 0;       // index in argument vector
190        size_t user_m     = 0;       // size of result vector
191        size_t user_n     = 0;       // size of arugment vector
192        //
193        atomic_base<Base>* user_atom = CPPAD_NULL; // user's atomic op calculator
194        bool               user_bool = false;      // use bool or set sparsity ?
195# ifndef NDEBUG
196        bool               user_ok   = false;      // atomic op return value
197# endif
198        // next expected operator in a UserOp sequence
199        enum { user_start, user_arg, user_ret, user_end } user_state = user_end;
200
201
202        // Initialize
203        play->start_reverse(op, arg, i_op, i_var);
204        CPPAD_ASSERT_UNKNOWN( op == EndOp );
205# if CPPAD_REV_HES_SWEEP_TRACE
206        std::cout << std::endl;
207        CppAD::vectorBool zf_value(limit);
208        CppAD::vectorBool zh_value(limit);
209# endif
210        while(op != BeginOp)
211        {
212                // next op
213                play->next_reverse(op, arg, i_op, i_var);
214# ifndef NDEBUG
215                if( i_op <= n )
216                {       CPPAD_ASSERT_UNKNOWN((op == InvOp) | (op == BeginOp));
217                }
218                else    CPPAD_ASSERT_UNKNOWN((op != InvOp) & (op != BeginOp));
219# endif
220
221                // rest of information depends on the case
222                switch( op )
223                {
224                        case AbsOp:
225                        CPPAD_ASSERT_NARG_NRES(op, 1, 1)
226                        reverse_sparse_hessian_linear_unary_op(
227                        i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
228                        );
229                        break;
230                        // -------------------------------------------------
231
232                        case AddvvOp:
233                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
234                        reverse_sparse_hessian_addsub_op(
235                        i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse
236                        );
237                        break;
238                        // -------------------------------------------------
239
240                        case AddpvOp:
241                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
242                        reverse_sparse_hessian_linear_unary_op(
243                        i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse
244                        );
245                        break;
246                        // -------------------------------------------------
247
248                        case AcosOp:
249                        // sqrt(1 - x * x), acos(x)
250                        CPPAD_ASSERT_NARG_NRES(op, 1, 2)
251                        reverse_sparse_hessian_nonlinear_unary_op(
252                        i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
253                        );
254                        break;
255                        // -------------------------------------------------
256
257                        case AsinOp:
258                        // sqrt(1 - x * x), asin(x)
259                        CPPAD_ASSERT_NARG_NRES(op, 1, 2)
260                        reverse_sparse_hessian_nonlinear_unary_op(
261                        i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
262                        );
263                        break;
264                        // -------------------------------------------------
265
266                        case AtanOp:
267                        // 1 + x * x, atan(x)
268                        CPPAD_ASSERT_NARG_NRES(op, 1, 2)
269                        reverse_sparse_hessian_nonlinear_unary_op(
270                        i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
271                        );
272                        break;
273                        // -------------------------------------------------
274
275                        case BeginOp:
276                        CPPAD_ASSERT_NARG_NRES(op, 0, 1)
277                        break;
278                        // -------------------------------------------------
279
280                        case CSumOp:
281                        // CSumOp has a variable number of arguments and
282                        // next_reverse thinks it one has one argument.
283                        // We must inform next_reverse of this special case.
284                        play->reverse_csum(op, arg, i_op, i_var);
285                        reverse_sparse_hessian_csum_op(
286                                i_var, arg, RevJac, rev_hes_sparse
287                        );
288                        break;
289                        // -------------------------------------------------
290
291                        case CExpOp:
292                        reverse_sparse_hessian_cond_op(
293                                i_var, arg, num_par, RevJac, rev_hes_sparse
294                        );
295                        break;
296                        // ---------------------------------------------------
297
298                        case ComOp:
299                        CPPAD_ASSERT_NARG_NRES(op, 4, 0)
300                        CPPAD_ASSERT_UNKNOWN( arg[1] > 1 );
301                        break;
302                        // --------------------------------------------------
303
304                        case CosOp:
305                        // sin(x), cos(x)
306                        CPPAD_ASSERT_NARG_NRES(op, 1, 2)
307                        reverse_sparse_hessian_nonlinear_unary_op(
308                        i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
309                        );
310                        break;
311                        // ---------------------------------------------------
312
313                        case CoshOp:
314                        // sinh(x), cosh(x)
315                        CPPAD_ASSERT_NARG_NRES(op, 1, 2)
316                        reverse_sparse_hessian_nonlinear_unary_op(
317                        i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
318                        );
319                        break;
320                        // -------------------------------------------------
321
322                        case DisOp:
323                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
324                        // derivativve is identically zero
325                        break;
326                        // -------------------------------------------------
327
328                        case DivvvOp:
329                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
330                        reverse_sparse_hessian_div_op(
331                        i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse
332                        );
333                        break;
334                        // -------------------------------------------------
335
336                        case DivpvOp:
337                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
338                        reverse_sparse_hessian_nonlinear_unary_op(
339                        i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse
340                        );
341                        break;
342                        // -------------------------------------------------
343
344                        case DivvpOp:
345                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
346                        reverse_sparse_hessian_linear_unary_op(
347                        i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
348                        );
349                        break;
350                        // -------------------------------------------------
351
352                        case ExpOp:
353                        CPPAD_ASSERT_NARG_NRES(op, 1, 1)
354                        reverse_sparse_hessian_nonlinear_unary_op(
355                        i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
356                        );
357                        break;
358                        // -------------------------------------------------
359
360                        case InvOp:
361                        CPPAD_ASSERT_NARG_NRES(op, 0, 1)
362                        // Z is already defined
363                        break;
364                        // -------------------------------------------------
365
366                        case LdpOp:
367                        reverse_sparse_hessian_load_op(
368                                op,
369                                i_var,
370                                arg,
371                                num_vecad_ind,
372                                vecad_ind.data(),
373                                rev_hes_sparse,
374                                vecad_sparse,
375                                RevJac,
376                                vecad_jac.data()
377                        );
378                        break;
379                        // -------------------------------------------------
380
381                        case LdvOp:
382                        reverse_sparse_hessian_load_op(
383                                op,
384                                i_var,
385                                arg,
386                                num_vecad_ind,
387                                vecad_ind.data(),
388                                rev_hes_sparse,
389                                vecad_sparse,
390                                RevJac,
391                                vecad_jac.data()
392                        );
393                        break;
394                        // -------------------------------------------------
395
396                        case LogOp:
397                        CPPAD_ASSERT_NARG_NRES(op, 1, 1)
398                        reverse_sparse_hessian_nonlinear_unary_op(
399                        i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
400                        );
401                        break;
402                        // -------------------------------------------------
403
404                        case MulvvOp:
405                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
406                        reverse_sparse_hessian_mul_op(
407                        i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse
408                        );
409                        break;
410                        // -------------------------------------------------
411
412                        case MulpvOp:
413                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
414                        reverse_sparse_hessian_linear_unary_op(
415                        i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse
416                        );
417                        break;
418                        // -------------------------------------------------
419
420                        case ParOp:
421                        CPPAD_ASSERT_NARG_NRES(op, 1, 1)
422
423                        break;
424                        // -------------------------------------------------
425
426                        case PowpvOp:
427                        CPPAD_ASSERT_NARG_NRES(op, 2, 3)
428                        reverse_sparse_hessian_nonlinear_unary_op(
429                        i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse
430                        );
431                        break;
432                        // -------------------------------------------------
433
434                        case PowvpOp:
435                        CPPAD_ASSERT_NARG_NRES(op, 2, 3)
436                        reverse_sparse_hessian_nonlinear_unary_op(
437                        i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
438                        );
439                        break;
440                        // -------------------------------------------------
441
442                        case PowvvOp:
443                        CPPAD_ASSERT_NARG_NRES(op, 2, 3)
444                        reverse_sparse_hessian_pow_op(
445                        i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse
446                        );
447                        break;
448                        // -------------------------------------------------
449
450                        case PriOp:
451                        CPPAD_ASSERT_NARG_NRES(op, 5, 0);
452                        break;
453                        // -------------------------------------------------
454
455                        case SignOp:
456                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
457                        // Derivative is identiaclly zero
458                        break;
459                        // -------------------------------------------------
460
461                        case SinOp:
462                        // cos(x), sin(x)
463                        CPPAD_ASSERT_NARG_NRES(op, 1, 2)
464                        reverse_sparse_hessian_nonlinear_unary_op(
465                        i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
466                        );
467                        break;
468                        // -------------------------------------------------
469
470                        case SinhOp:
471                        // cosh(x), sinh(x)
472                        CPPAD_ASSERT_NARG_NRES(op, 1, 2)
473                        reverse_sparse_hessian_nonlinear_unary_op(
474                        i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
475                        );
476                        break;
477                        // -------------------------------------------------
478
479                        case SqrtOp:
480                        CPPAD_ASSERT_NARG_NRES(op, 1, 1)
481                        reverse_sparse_hessian_nonlinear_unary_op(
482                        i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
483                        );
484                        break;
485                        // -------------------------------------------------
486
487                        case StppOp:
488                        // sparsity cannot propagate through a parameter
489                        CPPAD_ASSERT_NARG_NRES(op, 3, 0)
490                        break;
491                        // -------------------------------------------------
492
493                        case StpvOp:
494                        reverse_sparse_hessian_store_op(
495                                op,
496                                arg,
497                                num_vecad_ind,
498                                vecad_ind.data(),
499                                rev_hes_sparse,
500                                vecad_sparse,
501                                RevJac,
502                                vecad_jac.data()
503                        );
504                        break;
505                        // -------------------------------------------------
506
507                        case StvpOp:
508                        // sparsity cannot propagate through a parameter
509                        CPPAD_ASSERT_NARG_NRES(op, 3, 0)
510                        break;
511                        // -------------------------------------------------
512
513                        case StvvOp:
514                        reverse_sparse_hessian_store_op(
515                                op,
516                                arg,
517                                num_vecad_ind,
518                                vecad_ind.data(),
519                                rev_hes_sparse,
520                                vecad_sparse,
521                                RevJac,
522                                vecad_jac.data()
523                        );
524                        break;
525                        // -------------------------------------------------
526
527                        case SubvvOp:
528                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
529                        reverse_sparse_hessian_addsub_op(
530                        i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse
531                        );
532                        break;
533                        // -------------------------------------------------
534
535                        case SubpvOp:
536                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
537                        reverse_sparse_hessian_linear_unary_op(
538                        i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse
539                        );
540                        break;
541                        // -------------------------------------------------
542
543                        case SubvpOp:
544                        CPPAD_ASSERT_NARG_NRES(op, 2, 1)
545                        reverse_sparse_hessian_linear_unary_op(
546                        i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
547                        );
548                        break;
549                        // -------------------------------------------------
550
551                        case TanOp:
552                        // tan(x)^2, tan(x)
553                        CPPAD_ASSERT_NARG_NRES(op, 1, 2)
554                        reverse_sparse_hessian_nonlinear_unary_op(
555                        i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
556                        );
557                        break;
558                        // -------------------------------------------------
559
560                        case TanhOp:
561                        // tanh(x)^2, tanh(x)
562                        CPPAD_ASSERT_NARG_NRES(op, 1, 2)
563                        reverse_sparse_hessian_nonlinear_unary_op(
564                        i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
565                        );
566                        break;
567                        // -------------------------------------------------
568
569                        case UserOp:
570                        // start or end an atomic operation sequence
571                        CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 );
572                        CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 );
573                        if( user_state == user_end )
574                        {       user_index = arg[0];
575                                user_id    = arg[1];
576                                user_n     = arg[2];
577                                user_m     = arg[3];
578                                user_atom  = atomic_base<Base>::class_object(user_index);
579# ifndef NDEBUG
580                                if( user_atom == CPPAD_NULL )
581                                {       std::string msg = 
582                                                atomic_base<Base>::class_name(user_index)
583                                                + ": atomic_base function has been deleted";
584                                        CPPAD_ASSERT_KNOWN(false, msg.c_str() );
585                                }
586# endif
587                                user_bool  = user_atom->sparsity() ==
588                                                        atomic_base<Base>::bool_sparsity_enum;
589                                user_ix.resize(user_n);
590                                user_vx.resize(user_n);
591                                user_s.resize(user_m);
592                                user_t.resize(user_n);
593
594                                // simpler to initialize all sparsity patterns as empty
595                                for(i = 0; i < user_m; i++)
596                                        user_s[i] = false;
597                                for(i = 0; i < user_n; i++)
598                                        user_t[i] = false;
599                                if( user_bool )
600                                {       bool_r.resize(user_n * user_q);
601                                        bool_u.resize(user_m * user_q);
602                                        bool_v.resize(user_n * user_q);
603                                        // simpler to initialize all patterns as empty
604                                        for(i = 0; i < user_m; i++)
605                                        {
606                                                for(j = 0; j < user_q; j++)
607                                                        bool_u[ i * user_q + j] = false;
608                                        }
609                                        for(i = 0; i < user_n; i++)
610                                        {
611                                                for(j = 0; j < user_q; j++)
612                                                {       bool_r[ i * user_q + j] = false;
613                                                        bool_v[ i * user_q + j] = false;
614                                                }
615                                        }
616                                }
617                                else
618                                {       set_r.resize(user_n);
619                                        set_u.resize(user_m);
620                                        set_v.resize(user_n);
621                                        for(i = 0; i < user_m; i++)
622                                                set_u[i].clear();
623                                        for(i = 0; i < user_n; i++)
624                                        {       set_r[i].clear();
625                                                set_v[i].clear();
626                                        }
627                                }
628                                user_j     = user_n;
629                                user_i     = user_m;
630                                user_state = user_ret;
631                        }
632                        else
633                        {       CPPAD_ASSERT_UNKNOWN( user_state == user_start );
634                                CPPAD_ASSERT_UNKNOWN( user_index == size_t(arg[0]) );
635                                CPPAD_ASSERT_UNKNOWN( user_id    == size_t(arg[1]) );
636                                CPPAD_ASSERT_UNKNOWN( user_n     == size_t(arg[2]) );
637                                CPPAD_ASSERT_UNKNOWN( user_m     == size_t(arg[3]) );
638                                user_state = user_end;
639
640                                // call users function for this operation
641                                user_atom->set_id(user_id);
642# ifdef NDEBUG
643                                if( user_bool )
644                                        user_atom->rev_sparse_hes(user_vx,
645                                                user_s, user_t, user_q, bool_r, bool_u, bool_v
646                                );
647                                else
648                                        user_atom->rev_sparse_hes(user_vx,
649                                                user_s, user_t, user_q, set_r, set_u, set_v
650                                );
651# else
652                                if( user_bool )
653                                        user_ok = user_atom->rev_sparse_hes(user_vx,
654                                                user_s, user_t, user_q, bool_r, bool_u, bool_v
655                                );
656                                else
657                                        user_ok = user_atom->rev_sparse_hes(user_vx,
658                                                user_s, user_t, user_q, set_r, set_u, set_v
659                                );
660                                if( ! user_ok )
661                                {       std::string msg = 
662                                                atomic_base<Base>::class_name(user_index)
663                                                + ": atomic_base.rev_sparse_hes: returned false";
664                                        CPPAD_ASSERT_KNOWN(false, msg.c_str() );
665                                }
666# endif
667                                for(i = 0; i < user_n; i++) if( user_ix[i] > 0 )
668                                {
669                                        size_t  i_x = user_ix[i];
670                                        if( user_t[i] )
671                                                RevJac[i_x] = true;
672                                        if( user_bool )
673                                        {
674                                                for(j = 0; j < user_q; j++)
675                                                        if( bool_v[ i * user_q + j ] )
676                                                                rev_hes_sparse.add_element(i_x, j);
677                                        }
678                                        else
679                                        {
680                                                set_itr = set_v[i].begin();
681                                                set_end = set_v[i].end();
682                                                while( set_itr != set_end )
683                                                        rev_hes_sparse.add_element(i_x, *set_itr++);
684                                        }
685                                }
686               }
687                        break;
688
689                        case UsrapOp:
690                        // parameter argument in an atomic operation sequence
691                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
692                        CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
693                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
694                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
695                        --user_j;
696                        user_ix[user_j] = 0;
697                        user_vx[user_j] = false;
698                        if( user_j == 0 )
699                                user_state = user_start;
700                        break;
701
702                        case UsravOp:
703                        // variable argument in an atomic operation sequence
704                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
705                        CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
706                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
707                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) <= i_var );
708                        CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
709                        --user_j;
710                        user_ix[user_j] = arg[0];
711                        user_vx[user_j] = true;
712                        for_jac_sparse.begin(arg[0]);
713                        i = for_jac_sparse.next_element();
714                        while( i < user_q )
715                        {       if( user_bool )
716                                        bool_r[ user_j * user_q + i ] = true;
717                                else
718                                        set_r[user_j].insert(i);
719                                i = for_jac_sparse.next_element();
720                        }
721                        if( user_j == 0 )
722                                user_state = user_start;
723                        break;
724
725                        case UsrrpOp:
726                        // parameter result in an atomic operation sequence
727                        CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
728                        CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
729                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
730                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
731                        --user_i;
732                        if( user_i == 0 )
733                                user_state = user_arg;
734                        break;
735
736                        case UsrrvOp:
737                        // variable result in an atomic operation sequence
738                        CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
739                        CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
740                        --user_i;
741                        if( RevJac[i_var] )
742                        {
743                                user_s[user_i] = true;
744                        }
745                        rev_hes_sparse.begin(i_var);
746                        j = rev_hes_sparse.next_element();
747                        while( j < user_q )
748                        {       if( user_bool )
749                                        bool_u[user_i * user_q + j] = true;
750                                else
751                                        set_u[user_i].insert(j);
752                                j = rev_hes_sparse.next_element();
753                        }
754                        if( user_i == 0 )
755                                user_state = user_arg;
756                        break;
757
758                        // -------------------------------------------------
759
760                        default:
761                        CPPAD_ASSERT_UNKNOWN(0);
762                }
763# if CPPAD_REV_HES_SWEEP_TRACE
764                for(j = 0; j < limit; j++)
765                {       zf_value[j] = false;
766                        zh_value[j] = false;
767                }
768                for_jac_sparse.begin(i_var);;
769                j = for_jac_sparse.next_element();;
770                while( j < limit )
771                {       zf_value[j] = true;
772                        j = for_jac_sparse.next_element();
773                }
774                rev_hes_sparse.begin(i_var);;
775                j = rev_hes_sparse.next_element();;
776                while( j < limit )
777                {       zh_value[j] = true;
778                        j = rev_hes_sparse.next_element();
779                }
780                // should also print RevJac[i_var], but printOp does not
781                // yet allow for this.
782                printOp(
783                        std::cout, 
784                        play,
785                        i_op,
786                        i_var,
787                        op, 
788                        arg,
789                        1, 
790                        &zf_value, 
791                        1, 
792                        &zh_value
793                );
794# endif
795        }
796        // values corresponding to BeginOp
797        CPPAD_ASSERT_UNKNOWN( i_op == 0 );
798        CPPAD_ASSERT_UNKNOWN( i_var == 0 );
799
800        return;
801}
802/*! \} */
803} // END_CPPAD_NAMESPACE
804
805// preprocessor symbols that are local to this file
806# undef CPPAD_REV_HES_SWEEP_TRACE
807
808# endif
Note: See TracBrowser for help on using the repository browser.