source: trunk/cppad/local/forward_sweep.hpp @ 2625

Last change on this file since 2625 was 2625, checked in by bradbell, 7 years ago

Change comment so end of group is included in doxygen processing
(do not know why works without this command at end).

ad_assign.hpp: fix doxygen brief description.
ad_ctor.hpp: fix doxygen brief description.

  • Property svn:keywords set to Id
File size: 18.1 KB
Line 
1/* $Id: forward_sweep.hpp 2625 2012-12-23 14:34:12Z bradbell $ */
2# ifndef CPPAD_FORWARD_SWEEP_INCLUDED
3# define CPPAD_FORWARD_SWEEP_INCLUDED
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 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
16CPPAD_BEGIN_NAMESPACE
17/*!
18\defgroup forward_sweep_hpp forward_sweep.hpp
19\{
20\file forward_sweep.hpp
21Compute zero order forward mode Taylor coefficients.
22*/
23
24/*!
25\def CPPAD_FORWARD_SWEEP_TRACE
26This value is either zero or one.
27Zero is the normal operational value.
28If it is one, a trace of every forward_sweep computation is printed.
29*/
30# define CPPAD_FORWARD_SWEEP_TRACE 0
31
32/*!
33Compute arbitrary order forward mode Taylor coefficients.
34
35\tparam Base
36base type for the operator; i.e., this operation sequence was recorded
37using AD< \a Base > and computations by this routine are done using type
38\a Base.
39
40\param s_out
41Is the stream where output corresponding to \c PriOp operations will
42be written.
43
44\param print
45If \a print is false,
46suppress the output that is otherwise generated by the PriOp instructions
47(must be false when \c d is nonzero).
48
49\param d
50is the order of the Taylor coefficients that are computed during this call.
51
52\param n
53is the number of independent variables on the tape.
54
55\param numvar
56is the total number of variables on the tape.
57This is also equal to the number of rows in the matrix \a Taylor; i.e.,
58Rec->num_rec_var().
59
60\param Rec
61The information stored in \a Rec
62is a recording of the operations corresponding to the function
63\f[
64        F : {\bf R}^n \rightarrow {\bf R}^m
65\f]
66where \f$ n \f$ is the number of independent variables and
67\f$ m \f$ is the number of dependent variables.
68\n
69\n
70The object \a Rec is effectly constant.
71There are two exceptions to this.
72The first exception is that while palying back the tape
73the object \a Rec holds information about the current location
74with in the tape and this changes during palyback.
75The second exception is the fact that the
76zero order ( \a d = 0 ) versions of the VecAD operators LdpOp and LdvOp
77modify the corresponding \a op_arg values returned by
78\ref player::next_forward and \ref player::next_reverse; see the
79\link load_op.hpp LdpOp and LdvOp \endlink operations.
80
81\param J
82Is the number of columns in the coefficient matrix \a Taylor.
83This must be greater than or equal \a d + 1.
84
85\param Taylor
86\b Input: For j = 1 , ... , \a n, and for k = 0 , ... , \a d,
87\a Taylor [ j * J + k ]
88is the k-th order Taylor coefficient corresponding to
89variable with index j on the tape
90(independent variable with index (j-1) in the independent variable vector).
91\n
92\n
93\b Output: For i = \a n + 1, ... , \a numvar - 1, and for k = 0 , ... , \a d,
94\a Taylor [ i * J + k ]
95is the k-th order Taylor coefficient for the variable with
96index i on the tape.
97
98\a return
99If \a d is not zero, the return value is zero.
100If \a d is zero,
101the return value is equal to the number of ComOp operations
102that have a different result from when the information in
103\a Rec was recorded.
104(Note that if NDEBUG is true, there are no ComOp operations
105in Rec and hence this return value is always zero.)
106*/
107
108template <class Base>
109size_t forward_sweep(
110        std::ostream&         s_out,
111        bool                  print,
112        size_t                d,
113        size_t                n,
114        size_t                numvar,
115        player<Base>         *Rec,
116        size_t                J,
117        Base                 *Taylor
118)
119{       CPPAD_ASSERT_UNKNOWN( J >= d + 1 );
120
121        // op code for current instruction
122        OpCode           op;
123
124        // index for current instruction
125        size_t         i_op;
126
127        // next variables
128        size_t        i_var;
129
130        CPPAD_ASSERT_UNKNOWN( d == 0 || ! print );
131# if CPPAD_USE_FORWARD0SWEEP
132        CPPAD_ASSERT_UNKNOWN( d > 0 );
133# else
134        addr_t*         non_const_arg;
135# endif
136        const addr_t*   arg = 0;
137
138        // temporary indices
139        size_t i, ell;
140
141        // initialize the comparision operator (ComOp) counter
142        size_t compareCount = 0;
143
144        // if this is an order zero calculation, initialize vector indices
145        pod_vector<size_t> VectorInd;  // address for each element
146        pod_vector<bool>   VectorVar;  // is element a variable
147        i = Rec->num_rec_vecad_ind();
148        if( i > 0 )
149        {       VectorInd.extend(i);
150                VectorVar.extend(i);
151                while(i--)
152                {       VectorInd[i] = Rec->GetVecInd(i);
153                        VectorVar[i] = false;
154                }
155        }
156
157        // work space used by UserOp.
158        const size_t user_k  = d;    // order of this forward mode calculation
159        const size_t user_k1 = d+1;  // number of orders for this calculation
160        vector<Base> user_tx;        // argument vector Taylor coefficients
161        vector<Base> user_ty;        // result vector Taylor coefficients
162        vector<size_t> user_iy;      // variable indices for results vector
163        size_t user_index = 0;       // indentifier for this user_atomic operation
164        size_t user_id    = 0;       // user identifier for this call to operator
165        size_t user_i     = 0;       // index in result vector
166        size_t user_j     = 0;       // index in argument vector
167        size_t user_m     = 0;       // size of result vector
168        size_t user_n     = 0;       // size of arugment vector
169        // next expected operator in a UserOp sequence
170        enum { user_start, user_arg, user_ret, user_end } user_state = user_start;
171
172        // check numvar argument
173        CPPAD_ASSERT_UNKNOWN( Rec->num_rec_var() == numvar );
174
175        // length of the parameter vector (used by CppAD assert macros)
176        const size_t num_par = Rec->num_rec_par();
177
178        // pointer to the beginning of the parameter vector
179        const Base* parameter = 0;
180        if( num_par > 0 )
181                parameter = Rec->GetPar();
182
183# if ! CPPAD_USE_FORWARD0SWEEP
184        // length of the text vector (used by CppAD assert macros)
185        const size_t num_text = Rec->num_rec_text();
186
187        // pointer to the beginning of the text vector
188        const char* text = 0;
189        if( num_text > 0 )
190                text = Rec->GetTxt(0);
191# endif
192
193        // skip the BeginOp at the beginning of the recording
194        Rec->start_forward(op, arg, i_op, i_var);
195        CPPAD_ASSERT_UNKNOWN( op == BeginOp );
196# if CPPAD_FORWARD_SWEEP_TRACE
197        std::cout << std::endl;
198# endif
199        while(op != EndOp)
200        {
201                // this op
202                Rec->next_forward(op, arg, i_op, i_var);
203                CPPAD_ASSERT_UNKNOWN( (i_op > n)  | (op == InvOp) ); 
204                CPPAD_ASSERT_UNKNOWN( (i_op <= n) | (op != InvOp) ); 
205
206                // action depends on the operator
207                switch( op )
208                {
209                        case AbsOp:
210                        forward_abs_op(d, i_var, arg[0], J, Taylor);
211                        break;
212                        // -------------------------------------------------
213
214                        case AddvvOp:
215                        forward_addvv_op(d, i_var, arg, parameter, J, Taylor);
216                        break;
217                        // -------------------------------------------------
218
219                        case AddpvOp:
220                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
221                        forward_addpv_op(d, i_var, arg, parameter, J, Taylor);
222                        break;
223                        // -------------------------------------------------
224
225                        case AcosOp:
226                        // sqrt(1 - x * x), acos(x)
227                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
228                        forward_acos_op(d, i_var, arg[0], J, Taylor);
229                        break;
230                        // -------------------------------------------------
231
232                        case AsinOp:
233                        // sqrt(1 - x * x), asin(x)
234                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
235                        forward_asin_op(d, i_var, arg[0], J, Taylor);
236                        break;
237                        // -------------------------------------------------
238
239                        case AtanOp:
240                        // 1 + x * x, atan(x)
241                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
242                        forward_atan_op(d, i_var, arg[0], J, Taylor);
243                        break;
244                        // -------------------------------------------------
245
246                        case CSumOp:
247                        // CSumOp has a variable number of arguments and
248                        // next_forward thinks it one has one argument.
249                        // we must inform next_forward of this special case.
250                        Rec->forward_csum(op, arg, i_op, i_var);
251                        forward_csum_op(
252                                d, i_var, arg, num_par, parameter, J, Taylor
253                        );
254                        break;
255                        // -------------------------------------------------
256
257                        case CExpOp:
258                        forward_cond_op(
259                                d, i_var, arg, num_par, parameter, J, Taylor
260                        );
261                        break;
262                        // ---------------------------------------------------
263
264                        case ComOp:
265# if ! USE_FORWARD0SWEEP
266                        if( d == 0 ) forward_comp_op_0(
267                        compareCount, arg, num_par, parameter, J, Taylor
268                        );
269# endif
270                        break;
271                        // ---------------------------------------------------
272
273                        case CosOp:
274                        // sin(x), cos(x)
275                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
276                        forward_cos_op(d, i_var, arg[0], J, Taylor);
277                        break;
278                        // ---------------------------------------------------
279
280                        case CoshOp:
281                        // sinh(x), cosh(x)
282                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
283                        forward_cosh_op(d, i_var, arg[0], J, Taylor);
284                        break;
285                        // -------------------------------------------------
286
287                        case DisOp:
288# if ! CPPAD_USE_FORWARD0SWEEP
289                        if( d == 0 )
290                                forward_dis_op_0(i_var, arg, J, Taylor);
291                        else
292# endif
293                        {       Taylor[ i_var * J + d] = Base(0);
294                        }
295                        break;
296                        // -------------------------------------------------
297
298                        case DivvvOp:
299                        forward_divvv_op(d, i_var, arg, parameter, J, Taylor);
300                        break;
301                        // -------------------------------------------------
302
303                        case DivpvOp:
304                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
305                        forward_divpv_op(d, i_var, arg, parameter, J, Taylor);
306                        break;
307                        // -------------------------------------------------
308
309                        case DivvpOp:
310                        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < num_par );
311                        forward_divvp_op(d, i_var, arg, parameter, J, Taylor);
312                        break;
313                        // -------------------------------------------------
314
315                        case EndOp:
316                        CPPAD_ASSERT_NARG_NRES(op, 0, 0);
317                        break;
318                        // -------------------------------------------------
319
320                        case ExpOp:
321                        forward_exp_op(d, i_var, arg[0], J, Taylor);
322                        break;
323                        // -------------------------------------------------
324
325                        case InvOp:
326                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 0 );
327                        break;
328                        // -------------------------------------------------
329
330                        case LdpOp:
331# if ! CPPAD_USE_FORWARD0SWEEP
332                        if( d == 0 )
333                        {       non_const_arg = Rec->forward_non_const_arg();
334                                forward_load_p_op_0(
335                                        i_var, 
336                                        non_const_arg, 
337                                        num_par, 
338                                        parameter, 
339                                        J, 
340                                        Taylor,
341                                        Rec->num_rec_vecad_ind(),
342                                        VectorVar.data(),
343                                        VectorInd.data()
344                                );
345                        }
346                        else
347# endif
348                        {       forward_load_op( op, d, i_var, arg, J, Taylor);
349                        }
350                        break;
351                        // -------------------------------------------------
352
353                        case LdvOp:
354# if ! CPPAD_USE_FORWARD0SWEEP
355                        if( d == 0 )
356                        {       non_const_arg = Rec->forward_non_const_arg();
357                                forward_load_v_op_0(
358                                        i_var, 
359                                        non_const_arg, 
360                                        num_par, 
361                                        parameter, 
362                                        J, 
363                                        Taylor,
364                                        Rec->num_rec_vecad_ind(),
365                                        VectorVar.data(),
366                                        VectorInd.data()
367                                );
368                        }
369                        else
370# endif
371                        {       forward_load_op( op, d, i_var, arg, J, Taylor);
372                        }
373                        break;
374                        // -------------------------------------------------
375
376                        case LogOp:
377                        forward_log_op(d, i_var, arg[0], J, Taylor);
378                        break;
379                        // -------------------------------------------------
380
381                        case MulvvOp:
382                        forward_mulvv_op(d, i_var, arg, parameter, J, Taylor);
383                        break;
384                        // -------------------------------------------------
385
386                        case MulpvOp:
387                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
388                        forward_mulpv_op(d, i_var, arg, parameter, J, Taylor);
389                        break;
390                        // -------------------------------------------------
391
392                        case ParOp:
393# if ! CPPAD_USE_FORWARD0SWEEP
394                        if( d == 0 ) forward_par_op_0(
395                                i_var, arg, num_par, parameter, J, Taylor
396                        );
397                        else
398# endif
399                        {       Taylor[ i_var * J + d] = Base(0); 
400                        }
401                        break;
402                        // -------------------------------------------------
403
404                        case PowvpOp:
405                        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < num_par );
406                        forward_powvp_op(d, i_var, arg, parameter, J, Taylor);
407                        break;
408                        // -------------------------------------------------
409
410                        case PowpvOp:
411                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
412                        forward_powpv_op(d, i_var, arg, parameter, J, Taylor);
413                        break;
414                        // -------------------------------------------------
415
416                        case PowvvOp:
417                        forward_powvv_op(d, i_var, arg, parameter, J, Taylor);
418                        break;
419                        // -------------------------------------------------
420
421                        case PriOp:
422# if ! CPPAD_USE_FORWARD0SWEEP
423                        if( print ) forward_pri_0(s_out,
424                                i_var, arg, num_text, text, num_par, parameter, J, Taylor
425                        );
426# endif
427                        break;
428                        // -------------------------------------------------
429
430                        case SignOp:
431                        // cos(x), sin(x)
432                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
433                        forward_sign_op(d, i_var, arg[0], J, Taylor);
434                        break;
435
436                        case SinOp:
437                        // cos(x), sin(x)
438                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
439                        forward_sin_op(d, i_var, arg[0], J, Taylor);
440                        break;
441                        // -------------------------------------------------
442
443                        case SinhOp:
444                        // cosh(x), sinh(x)
445                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
446                        forward_sinh_op(d, i_var, arg[0], J, Taylor);
447                        break;
448                        // -------------------------------------------------
449
450                        case SqrtOp:
451                        forward_sqrt_op(d, i_var, arg[0], J, Taylor);
452                        break;
453                        // -------------------------------------------------
454
455                        case StppOp:
456# if ! CPPAD_USE_FORWARD0SWEEP
457                        if( d == 0 )
458                        {       forward_store_pp_op_0(
459                                        i_var, 
460                                        arg, 
461                                        num_par, 
462                                        J, 
463                                        Taylor,
464                                        Rec->num_rec_vecad_ind(),
465                                        VectorVar.data(),
466                                        VectorInd.data()
467                                );
468                        }
469# endif
470                        break;
471                        // -------------------------------------------------
472
473                        case StpvOp:
474# if ! CPPAD_USE_FORWARD0SWEEP
475                        if( d == 0 )
476                        {       forward_store_pv_op_0(
477                                        i_var, 
478                                        arg, 
479                                        num_par, 
480                                        J, 
481                                        Taylor,
482                                        Rec->num_rec_vecad_ind(),
483                                        VectorVar.data(),
484                                        VectorInd.data()
485                                );
486                        }
487# endif
488                        break;
489                        // -------------------------------------------------
490
491                        case StvpOp:
492# if ! CPPAD_USE_FORWARD0SWEEP
493                        if( d == 0 )
494                        {       forward_store_vp_op_0(
495                                        i_var, 
496                                        arg, 
497                                        num_par, 
498                                        J, 
499                                        Taylor,
500                                        Rec->num_rec_vecad_ind(),
501                                        VectorVar.data(),
502                                        VectorInd.data()
503                                );
504                        }
505# endif
506                        break;
507                        // -------------------------------------------------
508
509                        case StvvOp:
510# if ! CPPAD_USE_FORWARD0SWEEP
511                        if( d == 0 )
512                        {       forward_store_vv_op_0(
513                                        i_var, 
514                                        arg, 
515                                        num_par, 
516                                        J, 
517                                        Taylor,
518                                        Rec->num_rec_vecad_ind(),
519                                        VectorVar.data(),
520                                        VectorInd.data()
521                                );
522                        }
523# endif
524                        break;
525                        // -------------------------------------------------
526
527                        case SubvvOp:
528                        forward_subvv_op(d, i_var, arg, parameter, J, Taylor);
529                        break;
530                        // -------------------------------------------------
531
532                        case SubpvOp:
533                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
534                        forward_subpv_op(d, i_var, arg, parameter, J, Taylor);
535                        break;
536                        // -------------------------------------------------
537
538                        case SubvpOp:
539                        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < num_par );
540                        forward_subvp_op(d, i_var, arg, parameter, J, Taylor);
541                        break;
542                        // -------------------------------------------------
543
544                        case TanOp:
545                        // tan(x)^2, tan(x)
546                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
547                        forward_tan_op(d, i_var, arg[0], J, Taylor);
548                        break;
549                        // -------------------------------------------------
550
551                        case TanhOp:
552                        // tanh(x)^2, tanh(x)
553                        CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
554                        forward_tanh_op(d, i_var, arg[0], J, Taylor);
555                        break;
556                        // -------------------------------------------------
557
558                        case UserOp:
559                        // start or end an atomic operation sequence
560                        CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 );
561                        CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 );
562                        if( user_state == user_start )
563                        {       user_index = arg[0];
564                                user_id    = arg[1];
565                                user_n     = arg[2];
566                                user_m     = arg[3];
567                                if(user_tx.size() < user_n * user_k1)
568                                        user_tx.resize(user_n * user_k1);
569                                if(user_ty.size() < user_m * user_k1)
570                                        user_ty.resize(user_m * user_k1);
571                                if(user_iy.size() < user_m)
572                                        user_iy.resize(user_m);
573                                user_j     = 0;
574                                user_i     = 0;
575                                user_state = user_arg;
576                        }
577                        else
578                        {       CPPAD_ASSERT_UNKNOWN( user_state == user_end );
579                                CPPAD_ASSERT_UNKNOWN( user_index == size_t(arg[0]) );
580                                CPPAD_ASSERT_UNKNOWN( user_id    == size_t(arg[1]) );
581                                CPPAD_ASSERT_UNKNOWN( user_n     == size_t(arg[2]) );
582                                CPPAD_ASSERT_UNKNOWN( user_m     == size_t(arg[3]) );
583                                user_state = user_start;
584
585                                // call users function for this operation
586                                user_atomic<Base>::forward(user_index, user_id,
587                                        user_k, user_n, user_m, user_tx, user_ty
588                                );
589                                for(i = 0; i < user_m; i++) if( user_iy[i] > 0 )
590                                        Taylor[ user_iy[i] * J + user_k ] = 
591                                                user_ty[ i * user_k1 + user_k ];
592                        }
593                        break;
594
595                        case UsrapOp:
596                        // parameter argument in an atomic operation sequence
597                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
598                        CPPAD_ASSERT_UNKNOWN( user_j < user_n );
599                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
600                        user_tx[user_j * user_k1 + 0] = parameter[ arg[0]];
601                        for(ell = 1; ell < user_k1; ell++)
602                                user_tx[user_j * user_k1 + ell] = Base(0);
603                        ++user_j;
604                        if( user_j == user_n )
605                                user_state = user_ret;
606                        break;
607
608                        case UsravOp:
609                        // variable argument in an atomic operation sequence
610                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
611                        CPPAD_ASSERT_UNKNOWN( user_j < user_n );
612                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) <= i_var );
613                        for(ell = 0; ell < user_k1; ell++)
614                                user_tx[user_j * user_k1 + ell] = Taylor[ arg[0] * J + ell];
615                        ++user_j;
616                        if( user_j == user_n )
617                                user_state = user_ret;
618                        break;
619
620                        case UsrrpOp:
621                        // parameter result in an atomic operation sequence
622                        CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
623                        CPPAD_ASSERT_UNKNOWN( user_i < user_m );
624                        user_iy[user_i] = 0;
625                        user_ty[user_i * user_k1 + 0] = parameter[ arg[0]];
626                        for(ell = 1; ell < user_k; ell++)
627                                user_ty[user_i * user_k1 + ell] = Base(0);
628                        user_i++;
629                        if( user_i == user_m )
630                                user_state = user_end;
631                        break;
632
633                        case UsrrvOp:
634                        // variable result in an atomic operation sequence
635                        CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
636                        CPPAD_ASSERT_UNKNOWN( user_i < user_m );
637                        user_iy[user_i] = i_var;
638                        for(ell = 0; ell < user_k; ell++)
639                                user_ty[user_i * user_k1 + ell] = Taylor[ i_var * J + ell];
640                        user_i++;
641                        if( user_i == user_m )
642                                user_state = user_end;
643                        break;
644                        // -------------------------------------------------
645
646                        default:
647                        CPPAD_ASSERT_UNKNOWN(0);
648                }
649# if CPPAD_FORWARD_SWEEP_TRACE
650                size_t       i_tmp  = i_var;
651                Base*        Z_tmp  = Taylor + J * i_var;
652                printOp(
653                        std::cout, 
654                        Rec,
655                        i_tmp,
656                        op, 
657                        arg,
658                        d + 1, 
659                        Z_tmp, 
660                        0, 
661                        (Base *) CPPAD_NULL
662                );
663        }
664        std::cout << std::endl;
665# else
666        }
667# endif
668        CPPAD_ASSERT_UNKNOWN( user_state == user_start );
669        CPPAD_ASSERT_UNKNOWN( i_var + 1 == Rec->num_rec_var() );
670
671        return compareCount;
672}
673
674# undef CPPAD_FORWARD_SWEEP_TRACE
675
676/*! \} */
677CPPAD_END_NAMESPACE
678# endif
Note: See TracBrowser for help on using the repository browser.