source: branches/cache/cppad/local/rev_hes_sweep.hpp @ 3334

Last change on this file since 3334 was 3334, checked in by bradbell, 6 years ago

Add chace2var argument to printOp (but not used yet).

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