source: trunk/ADOL-C/src/ho_rev.c @ 270

Last change on this file since 270 was 249, checked in by awalther, 9 years ago

correct typos

  • Property svn:keywords set to Author Date Id Revision
File size: 75.1 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     ho_rev.c
4 Revision: $Id: ho_rev.c 249 2011-06-13 18:57:34Z kulshres $
5 Contents: Contains the routines :
6           hos_reverse (higher-order-scalar reverse mode):
7              define _HOS_
8           hos_ov_reverse (higher-order-scalar reverse mode on vectors):
9              define _HOS_OV_
10           hov_reverse (higher-order-vector reverse mode):
11              define _HOV_
12 
13 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz,
14               Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel
15 
16 This file is part of ADOL-C. This software is provided as open source.
17 Any use, reproduction, or distribution of the software constitutes
18 recipient's acceptance of the terms of the accompanying license file.
19 
20----------------------------------------------------------------------------*/
21
22/*****************************************************************************
23 
24  There are four basic versions of the procedure `reverse', which
25  are optimized for the cases of scalar or vector reverse sweeps
26  with first or higher derivatives, respectively. In the calling
27  sequence this distinction is apparent from the type of the
28  parameters `lagrange' and `results'. The former may be left out
29  and the integer parameters `depen', `indep', `degre', and `nrows'
30  must be set or default according to the following matrix of
31  calling cases.
32 
33           no lagrange         double* lagrange     double** lagrange
34 
35double*   gradient of scalar   weight vector times    infeasible
36results   valued function      Jacobian product       combination
37 
38          ( depen = 1 ,         ( depen > 0 ,         
39            degre = 0 ,           degre = 0 ,              ------
40            nrows = 1 )           nrows = 1 )
41 
42double**  Jacobian of vector   weight vector times     weight matrix
43results   valued function      Taylor-Jacobians        times Jacobian
44           
45          ( 0 < depen           ( depen > 0 ,          ( depen > 0 ,
46              = nrows ,           degre > 0 ,            degre = 0 ,
47            degre = 0 )           nrows = 1 )            nrows > 0 )
48 
49double*** full family of         ------------          weigth matrix x
50results   Taylor-Jacobians       ------------          Taylor Jacobians
51 
52*****************************************************************************/
53
54/****************************************************************************/
55/*                                                                   MACROS */
56#undef _ADOLC_VECTOR_
57#undef _HIGHER_ORDER_
58
59/*--------------------------------------------------------------------------*/
60#ifdef _HOS_
61#define GENERATED_FILENAME "hos_reverse"
62
63#define _HIGHER_ORDER_
64
65#define RESULTS(l,indexi,k) results[indexi][k]
66#define LAGRANGE(l,indexd,k)  lagrange[indexd][k]
67
68#define HOV_INC(T,degree) {}
69#define HOS_OV_INC(T,degree) {}
70
71#define GET_TAYL(loc,depth,p) \
72    { \
73        UPDATE_TAYLORREAD(depth) \
74        get_taylors(loc,depth); \
75    }
76
77/*--------------------------------------------------------------------------*/
78#elif _HOS_OV_
79#define GENERATED_FILENAME "hos_ov_reverse"
80
81#define _HIGHER_ORDER_
82
83#define RESULTS(l,indexi,k) results[l][indexi][k]
84#define LAGRANGE(l,indexd,k)  lagrange[indexd][k]
85
86#define HOV_INC(T,degree) T += degree;
87#define HOS_OV_INC(T,degree) T += degree;
88
89#define GET_TAYL(loc,depth,p) \
90    { \
91        UPDATE_TAYLORREAD(depth * p) \
92        get_taylors_p(loc,depth,p); \
93    }
94
95/*--------------------------------------------------------------------------*/
96#elif _HOV_
97#define GENERATED_FILENAME "hov_reverse"
98
99#define _ADOLC_VECTOR_
100#define _HIGHER_ORDER_
101
102#define RESULTS(l,indexi,k) results[l][indexi][k]
103#define LAGRANGE(l,indexd,k)  lagrange[l][indexd][k]
104
105#define IF_HOV_
106#define ENDIF_HOV_
107
108#define HOV_INC(T,degree) T += degree;
109#define HOS_OV_INC(T,degree)
110
111#define GET_TAYL(loc,depth,p) \
112    { \
113        UPDATE_TAYLORREAD(depth) \
114        get_taylors(loc,depth); \
115    }
116
117#else
118#error Error ! Define [_HOS_ | _HOS_OV_ | _HOV_]
119#endif
120
121/*--------------------------------------------------------------------------*/
122/*                                                     access to variables  */
123
124#ifdef _FOS_                                     /* why?, not in fo_rev.c ? */
125#define ARES       *Ares
126#define AARG       *Aarg
127#define AARG1      *Aarg1
128#define AARG2      *Aarg2
129#define AQO        *Aqo
130
131#define ARES_INC   *Ares
132#define AARG_INC   *Aarg
133#define AARG1_INC  *Aarg1
134#define AARG2_INC  *Aarg2
135#define AQO_INC    *Aqo
136
137#define ARES_INC_O  Ares
138#define AARG_INC_O  Aarg
139#define AARG1_INC_O Aarg1
140#define AARG2_INC_O Aarg2
141#define AQO_INC_O   Aqo
142
143#define ASSIGN_A(a,b)  a = &b;
144#define HOS_OV_ASSIGN_A(Aqo,  rp_Atemp)
145#define FOR_0_LE_l_LT_q l = 0;
146
147#elif _HOS_OV_
148#define ARES       *Ares
149#define AARG       *Aarg
150#define AARG1      *Aarg1
151#define AARG2      *Aarg2
152#define AQO        *Aqo
153
154#define ARES_INC   *Ares++
155#define AARG_INC   *Aarg++
156#define AARG1_INC  *Aarg1++
157#define AARG2_INC  *Aarg2++
158#define AQO_INC    *Aqo++
159
160#define ARES_INC_O  Ares++
161#define AARG_INC_O  Aarg++
162#define AARG1_INC_O Aarg1++
163#define AARG2_INC_O Aarg2++
164#define AQO_INC_O   Aqo++
165
166#define ASSIGN_A(a,b)  a = b;
167#define HOS_OV_ASSIGN_A(a, b) a = b;
168#define FOR_0_LE_l_LT_q for(l=0;l<q;l++)
169
170#else  /* _FOV_, _HOS_, _HOV_ */
171#define ARES       *Ares
172#define AARG       *Aarg
173#define AARG1      *Aarg1
174#define AARG2      *Aarg2
175#define AQO        *Aqo
176
177#define ARES_INC   *Ares++
178#define AARG_INC   *Aarg++
179#define AARG1_INC  *Aarg1++
180#define AARG2_INC  *Aarg2++
181#define AQO_INC    *Aqo++
182
183#define ARES_INC_O  Ares++
184#define AARG_INC_O  Aarg++
185#define AARG1_INC_O Aarg1++
186#define AARG2_INC_O Aarg2++
187#define AQO_INC_O   Aqo++
188
189#define ASSIGN_A(a,b)  a = b;
190#define HOS_OV_ASSIGN_A(Aqo, Atemp)
191#define FOR_0_LE_l_LT_q l = 0;
192#endif
193
194#ifdef _HIGHER_ORDER_
195
196#define TRES      *Tres                  /* why ? not used here */
197#define TARG      *Targ
198#define TARG1     *Targ1
199#define TARG2     *Targ2
200
201#define ASSIGN_T(a,b)  a = b;
202#else
203
204#define TRES       rpp_T[res]
205#define TARG       rpp_T[arg]
206#define TARG1      rpp_T[arg1]
207#define TARG2      rpp_T[arg2]
208
209#define ASSIGN_T(a,b)
210#endif
211
212/*--------------------------------------------------------------------------*/
213/*                                                              loop stuff  */
214#ifdef _ADOLC_VECTOR_
215#define FOR_0_LE_l_LT_p for (l=0; l<p; l++)
216#define FOR_p_GT_l_GE_0 for (l=p-1; l>=0; l--)  /* why ? not used here */
217#elif _HOS_OV_
218#define FOR_0_LE_l_LT_p for (l=0; l<p; l++)
219#define FOR_p_GT_l_GE_0                         /* why ? not used here */
220#else
221#define FOR_0_LE_l_LT_p
222#define FOR_p_GT_l_GE_0                         /* why ? not used here */
223#endif
224
225#ifdef _HIGHER_ORDER_
226#define FOR_0_LE_i_LT_k for (i=0; i<k; i++)
227#define FOR_k_GT_i_GE_0 for (i=k-1; i>=0; i--)
228#else
229#define FOR_0_LE_i_LT_k
230#define FOR_k_GT_i_GE_0
231#endif
232
233#ifdef _HOV_
234#define FOR_0_LE_l_LT_pk1 for (l=0; l<pk1; l++)
235#define FOR_0_LE_l_LT_pk for (l=0; l<k; l++)
236#elif _FOV_
237#define FOR_0_LE_l_LT_pk1 for (l=0; l<p; l++)
238#define FOR_0_LE_l_LT_pk for (l=0; l<k; l++)
239#elif _HOS_
240#define FOR_0_LE_l_LT_pk1 for (l=0; l<k1; l++)
241#define FOR_0_LE_l_LT_pk for (l=0; l<k; l++)
242#elif _HOS_OV_
243#define FOR_0_LE_l_LT_pk1 for (l=0; l<pk1; l++)
244#define FOR_0_LE_l_LT_pk for (l=0; l<p*k; l++)
245#else
246#define FOR_0_LE_l_LT_pk1
247#define FOR_0_LE_l_LT_pk
248#endif
249
250/*--------------------------------------------------------------------------*/
251/*                                                         VEC_COMPUTED_* */
252#ifdef _ADOLC_VECTOR
253#define VEC_COMPUTED_INIT   computed = 0;
254#define VEC_COMPUTED_CHECK  if (computed == 0) { computed = 1;
255#define VEC_COMPUTED_END    }
256#else
257#define VEC_COMPUTED_INIT
258#define VEC_COMPUTED_CHECK
259#define VEC_COMPUTED_END
260#endif
261
262/* END Macros */
263
264
265/****************************************************************************/
266/*                                                       NECESSARY INCLUDES */
267#include <adolc/interfaces.h>
268#include <adolc/adalloc.h>
269#include <adolc/oplate.h>
270#include "taping_p.h"
271#include <adolc/convolut.h>
272
273#include <math.h>
274
275#if defined(ADOLC_DEBUG)
276#include <string.h>
277#endif /* ADOLC_DEBUG */
278
279BEGIN_C_DECLS
280
281/*--------------------------------------------------------------------------*/
282/*                                                   Local Static Variables */
283
284#ifdef _HOS_
285/***************************************************************************/
286/* Higher Order Scalar Reverse Pass.                                       */
287/***************************************************************************/
288int hos_reverse(short   tnum,        /* tape id */
289                int     depen,       /* consistency chk on # of deps */
290                int     indep,       /* consistency chk on # of indeps */
291                int     degre,       /* highest derivative degre  */
292                double  *lagrange,   /* range weight vector       */
293                double  **results)   /* matrix of coefficient vectors */
294{ int i, j, rc;
295    double** L = myalloc2(depen,degre+1);
296    for ( i = 0; i < depen; ++i ) {
297        L[i][0] = lagrange[i];
298        for ( j = 1; j <= degre; ++j )
299            L[i][j] = 0.0;
300    }
301    rc = hos_ti_reverse(tnum,depen,indep,degre,L,results);
302    myfree2(L);
303    return rc;
304}
305
306int hos_ti_reverse(
307    short   tnum,        /* tape id */
308    int     depen,       /* consistency chk on # of deps */
309    int     indep,       /* consistency chk on # of indeps */
310    int     degre,       /* highest derivative degre  */
311    double  **lagrange,  /* range weight vectors       */
312    double  **results)   /* matrix of coefficient vectors */
313
314#elif _HOS_OV_
315
316/***************************************************************************/
317/* Higher Order Scalar Reverse Pass, Vector Keep.                          */
318/***************************************************************************/
319int hos_ov_reverse(short   tnum,       /* tape id */
320                   int     depen,       /* consistency chk on # of deps */
321                   int     indep,       /* consistency chk on # of indeps */
322                   int     degre,       /* highest derivative degre  */
323                   int     nrows,       /* # of Jacobian rows calculated */
324                   double  **lagrange,  /* range weight vector       */
325                   double  ***results)  /* matrix of coefficient vectors */
326
327#elif _HOV_
328/***************************************************************************/
329/* Higher Order Vector Reverse Pass.                                       */
330/***************************************************************************/
331int hov_reverse(short   tnum,        /* tape id */
332                int     depen,       /* consistency chk on # of deps */
333                int     indep,       /* consistency chk on # of indeps */
334                int     degre,       /* highest derivative degre */
335                int     nrows,       /* # of Jacobian rows calculated */
336                double  **lagrange,  /* domain weight vector */
337                double  ***results,  /* matrix of coefficient vectors */
338                short   **nonzero )  /* structural sparsity  pattern  */
339{ int i, j, k, rc;
340    double*** L = myalloc3(nrows,depen,degre+1);
341    for ( k = 0; k < nrows; ++k )
342        for ( i = 0; i < depen; ++i ) {
343            L[k][i][0] = lagrange[k][i];
344            for ( j = 1; j <= degre; ++j )
345                L[k][i][j] = 0.0;
346        }
347    rc = hov_ti_reverse(tnum,depen,indep,degre,nrows,L,results,nonzero);
348    myfree3(L);
349    return rc;
350}
351
352int hov_ti_reverse(
353    short   tnum,        /* tape id */
354    int     depen,       /* consistency chk on # of deps */
355    int     indep,       /* consistency chk on # of indeps */
356    int     degre,       /* highest derivative degre */
357    int     nrows,       /* # of Jacobian rows calculated */
358    double  ***lagrange, /* domain weight vectors */
359    double  ***results,  /* matrix of coefficient vectors */
360    short   **nonzero )  /* structural sparsity  pattern  */
361
362#endif
363
364{
365    /************************************************************************/
366    /*                                                       ALL VARIABLES  */
367    unsigned char operation;   /* operation code */
368    int dc, ret_c=3;
369
370    locint size = 0;
371    locint res  = 0;
372    locint arg  = 0;
373    locint arg1 = 0;
374    locint arg2 = 0;
375
376    double coval = 0, *d = 0;
377
378    int indexi = 0,  indexd = 0;
379
380    /* loop indices */
381    int i, j, l, ls;
382
383    /* other necessary variables */
384    double *x;
385    int *jj;
386    int taycheck;
387    int numdep,numind;
388    double aTmp;
389
390    /*----------------------------------------------------------------------*/
391    /* Taylor stuff */
392    revreal *Tres, *Targ, *Targ1, *Targ2, *Tqo, *rp_Ttemp, *rp_Ttemp2;
393    revreal **rpp_T;
394
395    /*----------------------------------------------------------------------*/
396    /* Adjoint stuff */
397#ifdef _FOS_
398    double Atemp;
399# define A_TEMP Atemp
400#endif
401    revreal *Ares, *Aarg=NULL, *Aarg1, *Aarg2, *Aqo, *rp_Atemp, *rp_Atemp2;
402    revreal **rpp_A, *AP1, *AP2;
403
404    /*----------------------------------------------------------------------*/
405    int k = degre + 1;
406    int k1 = k + 1;
407    revreal comp;
408
409#ifdef _ADOLC_VECTOR_
410    int p = nrows;
411#endif
412
413#ifdef _HOV_
414    int pk1 = p*k1;
415    int q = 1;
416#elif _HOS_OV_
417    int p = nrows;
418    int pk1 = p*k1;
419    int q = p;
420#else
421    int q = 1;
422#endif
423
424
425    ADOLC_OPENMP_THREAD_NUMBER;
426    ADOLC_OPENMP_GET_THREAD_NUMBER;
427
428#if defined(ADOLC_DEBUG)
429    /************************************************************************/
430    /*                                                       DEBUG MESSAGES */
431    fprintf(DIAG_OUT,"Call of %s(..) with tag: %d, n: %d, m %d,\n",
432            GENERATED_FILENAME, tnum, indep, depen);
433
434#ifdef _HIGHER_ORDER_
435    fprintf(DIAG_OUT,"                    degree: %d\n",degre);
436#endif
437#ifdef _ADOLC_VECTOR_
438    fprintf(DIAG_OUT,"                    p: %d\n\n",nrows);
439#endif
440
441#endif
442
443
444    /************************************************************************/
445    /*                                                                INITs */
446
447    /*----------------------------------------------------------------------*/
448    /* Set up stuff for the tape */
449
450    /* Initialize the Reverse Sweep */
451    init_rev_sweep(tnum);
452
453    if ( (depen != ADOLC_CURRENT_TAPE_INFOS.stats[NUM_DEPENDENTS]) ||
454            (indep != ADOLC_CURRENT_TAPE_INFOS.stats[NUM_INDEPENDENTS]) )
455        fail(ADOLC_REVERSE_COUNTS_MISMATCH);
456
457    indexi = ADOLC_CURRENT_TAPE_INFOS.stats[NUM_INDEPENDENTS] - 1;
458    indexd = ADOLC_CURRENT_TAPE_INFOS.stats[NUM_DEPENDENTS] - 1;
459
460
461    /************************************************************************/
462    /*                                              MEMORY ALLOCATION STUFF */
463
464    /*----------------------------------------------------------------------*/
465#ifdef _HOS_                                                         /* HOS */
466    rpp_A = (revreal**)malloc(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES] *
467            sizeof(revreal*));
468    if (rpp_A == NULL) fail(ADOLC_MALLOC_FAILED);
469    Aqo = (revreal*)malloc(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES] * k1 *
470            sizeof(revreal));
471    if (Aqo == NULL) fail(ADOLC_MALLOC_FAILED);
472    for (i=0; i<ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES]; i++) {
473        rpp_A[i] = Aqo;
474        Aqo += k1;
475    }
476    rpp_T = (revreal**)malloc(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES] *
477            sizeof(revreal*));
478    if (rpp_T == NULL) fail(ADOLC_MALLOC_FAILED);
479    Tqo = (revreal *)malloc(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES] *
480            k * sizeof(revreal));
481    if (Tqo ==NULL) fail(ADOLC_MALLOC_FAILED);
482    for (i=0; i<ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES]; i++) {
483        rpp_T[i] = Tqo;
484        Tqo += k;
485    }
486    rp_Atemp  = (revreal *)malloc(k1 * sizeof(revreal));
487    rp_Atemp2 = (revreal *)malloc(k1 * sizeof(revreal));
488    rp_Ttemp2 = (revreal *)malloc(k * sizeof(revreal));
489    /*----------------------------------------------------------------------*/
490#elif _HOV_                                                          /* HOV */
491    rpp_A = (revreal**)malloc(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES] *
492            sizeof(revreal*));
493    if (rpp_A == NULL) fail(ADOLC_MALLOC_FAILED);
494    Aqo = (revreal*)malloc(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES] * pk1 *
495            sizeof(revreal));
496    if (Aqo == NULL) fail(ADOLC_MALLOC_FAILED);
497    for (i=0; i<ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES]; i++) {
498        rpp_A[i] = Aqo;
499        Aqo += pk1;
500    }
501    rpp_T = (revreal**)malloc(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES] *
502            sizeof(revreal*));
503    if (rpp_T == NULL) fail(ADOLC_MALLOC_FAILED);
504    Tqo = (revreal*)malloc(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES] *
505            k * sizeof(revreal));
506    if (Tqo == NULL) fail(ADOLC_MALLOC_FAILED);
507    for (i=0; i<ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES]; i++) {
508        rpp_T[i] = Tqo;
509        Tqo += k;
510    }
511    rp_Atemp  = (revreal*) malloc(pk1 * sizeof(revreal));
512    rp_Atemp2 = (revreal*) malloc(pk1 * sizeof(revreal));
513    rp_Ttemp2 = (revreal*) malloc(k * sizeof(revreal));
514    /*----------------------------------------------------------------------*/
515#elif _HOS_OV_                                                    /* HOS_OV */
516    rpp_A = (revreal**)malloc(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES] *
517            sizeof(revreal*));
518    if (rpp_A == NULL) fail(ADOLC_MALLOC_FAILED);
519    Aqo = (revreal*)malloc(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES] * pk1 *
520            sizeof(revreal));
521    if (Aqo == NULL) fail(ADOLC_MALLOC_FAILED);
522    for (i=0; i<ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES]; i++) {
523        rpp_A[i] = Aqo;
524        Aqo += pk1;
525    }
526    rpp_T = (revreal**)malloc(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES] *
527            sizeof(revreal*));
528    if (rpp_T == NULL) fail(ADOLC_MALLOC_FAILED);
529    Tqo = (revreal*)malloc(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES] *
530            p * k * sizeof(revreal));
531    if (Tqo == NULL) fail(ADOLC_MALLOC_FAILED);
532    for (i=0; i<ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES]; i++) {
533        rpp_T[i] = Tqo;
534        Tqo += p*k;
535    }
536    rp_Atemp  = (revreal*) malloc(pk1 * sizeof(revreal));
537    rp_Atemp2 = (revreal*) malloc(pk1 * sizeof(revreal));
538    rp_Ttemp2 = (revreal*) malloc(p*k*sizeof(revreal));
539#endif
540    rp_Ttemp  = (revreal*) malloc(k*sizeof(revreal));
541    if (rp_Ttemp == NULL) fail(ADOLC_MALLOC_FAILED);
542    if (rp_Ttemp2 == NULL) fail(ADOLC_MALLOC_FAILED);
543    x = myalloc1(q);
544    jj = (int*)malloc(q*sizeof(int));
545    if (jj == NULL) fail(ADOLC_MALLOC_FAILED);
546
547    /************************************************************************/
548    /*                                                TAYLOR INITIALIZATION */
549    ADOLC_CURRENT_TAPE_INFOS.rpp_A = rpp_A;
550    ADOLC_CURRENT_TAPE_INFOS.rpp_T = rpp_T;
551    taylor_back(tnum,&numdep,&numind,&taycheck);
552
553    if(taycheck != degre) {
554        fprintf(DIAG_OUT,"\n ADOL-C error: reverse fails because it was not"
555                " preceded\nby a forward sweep with degree>%i,"
556                " keep=%i!\n",degre,degre+1);
557        exit(-2);
558    };
559
560    if((numdep != depen)||(numind != indep)) {
561        fprintf(DIAG_OUT,"\n ADOL-C error: reverse fails on tape %d because "
562                "the number of\nindependent and/or dependent variables"
563                " given to reverse are\ninconsistent with that of the"
564                "  internal taylor array.\n",tnum);
565        exit(-2);
566    }
567
568
569    /************************************************************************/
570    /*                                                        REVERSE SWEEP */
571
572#if defined(ADOLC_DEBUG)
573    int v = 0;
574    unsigned int countPerOperation[256], taylorPerOperation[256];
575    memset(countPerOperation, 0, 1024);
576    memset(taylorPerOperation, 0, 1024);
577#   define UPDATE_TAYLORREAD(X) taylorPerOperation[operation] += X;
578#else
579#   define UPDATE_TAYLORREAD(X)
580#endif /* ADOLC_DEBUG */
581
582    operation=get_op_r();
583#if defined(ADOLC_DEBUG)
584    ++countPerOperation[operation];
585#endif /* ADOLC_DEBUG */
586
587    while (operation != start_of_tape) { 
588        /* Switch statement to execute the operations in Reverse */
589        switch (operation) {
590
591
592                /************************************************************/
593                /*                                                  MARKERS */
594
595                /*----------------------------------------------------------*/
596            case end_of_op:                                    /* end_of_op */
597                get_op_block_r();
598                operation = get_op_r();
599                /* Skip next operation, it's another end_of_op */
600                break;
601
602                /*----------------------------------------------------------*/
603            case end_of_int:                                  /* end_of_int */
604                get_loc_block_r(); /* Get the next int block */
605                break;
606
607                /*----------------------------------------------------------*/
608            case end_of_val:                                  /* end_of_val */
609                get_val_block_r(); /* Get the next val block */
610                break;
611
612                /*----------------------------------------------------------*/
613            case start_of_tape:                            /* start_of_tape */
614            case end_of_tape:                                /* end_of_tape */
615                break;
616
617
618                /************************************************************/
619                /*                                               COMPARISON */
620
621                /*----------------------------------------------------------*/
622            case eq_zero  :                                      /* eq_zero */
623                arg   = get_locint_r();
624
625                ret_c = 0;
626                break;
627
628                /*----------------------------------------------------------*/
629            case neq_zero :                                     /* neq_zero */
630            case gt_zero  :                                      /* gt_zero */
631            case lt_zero :                                       /* lt_zero */
632                arg   = get_locint_r();
633                break;
634
635                /*----------------------------------------------------------*/
636            case ge_zero :                                       /* ge_zero */
637            case le_zero :                                       /* le_zero */
638                arg   = get_locint_r();
639
640                if (*rpp_T[arg] == 0)
641                    ret_c = 0;
642                break;
643
644
645                /************************************************************/
646                /*                                              ASSIGNMENTS */
647
648                /*----------------------------------------------------------*/
649            case assign_a:     /* assign an adouble variable an    assign_a */
650                /* adouble value. (=) */
651                res = get_locint_r();
652                arg = get_locint_r();
653
654                ASSIGN_A(Aarg, rpp_A[arg])
655                ASSIGN_A(Ares, rpp_A[res])
656
657                FOR_0_LE_l_LT_p
658                if  (0 == ARES) {
659                    HOV_INC(Aarg, k1)
660                    HOV_INC(Ares, k1)
661                } else {
662                    MAXDEC(AARG,ARES);
663                    AARG_INC_O;
664                    ARES_INC = 0.0;
665                    FOR_0_LE_i_LT_k
666                    { /* ! no tempory */
667                        AARG_INC += ARES;
668                        ARES_INC = 0.0;
669                    }
670                }
671
672                GET_TAYL(res,k,p)
673                break;
674
675                /*----------------------------------------------------------*/
676            case assign_d:      /* assign an adouble variable a    assign_d */
677                /* double value. (=) */
678                res   = get_locint_r();
679                coval = get_val_r();
680
681                ASSIGN_A(Ares, rpp_A[res])
682
683                FOR_0_LE_l_LT_pk1
684                ARES_INC = 0.0;
685
686                GET_TAYL(res,k,p)
687                break;
688
689                /*----------------------------------------------------------*/
690            case assign_d_zero: /* assign an adouble a        assign_d_zero */
691            case assign_d_one:  /* double value. (=)           assign_d_one */
692                res   = get_locint_r();
693
694                ASSIGN_A(Ares, rpp_A[res])
695
696                FOR_0_LE_l_LT_pk1
697                ARES_INC = 0.0;
698
699                GET_TAYL(res,k,p)
700                break;
701
702                /*--------------------------------------------------------------------------*/
703            case assign_ind:       /* assign an adouble variable an    assign_ind */
704                /* independent double value (<<=) */
705                res = get_locint_r();
706
707                ASSIGN_A(Ares, rpp_A[res])
708
709                FOR_0_LE_l_LT_p
710                {
711#ifdef _HOV_
712                    if (nonzero) /* ??? question: why here? */
713                    nonzero[l][indexi] = (int)ARES;
714#endif /* _HOV_ */
715                    ARES_INC_O;
716                    FOR_0_LE_i_LT_k
717                        RESULTS(l,indexi,i) = ARES_INC;
718                }
719
720                GET_TAYL(res,k,p)
721                    indexi--;
722                break;
723
724                /*--------------------------------------------------------------------------*/
725            case assign_dep:           /* assign a float variable a    assign_dep */
726                /* dependent adouble value. (>>=) */
727                res = get_locint_r();
728
729                ASSIGN_A(Ares, rpp_A[res])
730                ASSIGN_A(Aarg, rpp_A[res])   /* just a helpful pointers */
731
732                FOR_0_LE_l_LT_p
733                { ARES_INC_O;
734                  dc = -1;
735                  FOR_0_LE_i_LT_k
736                  { ARES_INC = LAGRANGE(l,indexd,i);
737                    if (LAGRANGE(l,indexd,i)) dc = i;
738                  }
739                  AARG = (dc < 0)? 0.0 : (dc > 0)? 2.0 : 1.0;
740                  HOV_INC(Aarg, k1)
741                }
742                indexd--;
743            break;
744
745
746            /****************************************************************************/
747            /*                                                   OPERATION + ASSIGNMENT */
748
749            /*--------------------------------------------------------------------------*/
750        case eq_plus_d:            /* Add a floating point to an    eq_plus_d */
751            /* adouble. (+=) */
752            res   = get_locint_r();
753                coval = get_val_r();
754
755                GET_TAYL(res,k,p)
756                break;
757
758                /*--------------------------------------------------------------------------*/
759            case eq_plus_a:             /* Add an adouble to another    eq_plus_a */
760                /* adouble. (+=) */
761                res = get_locint_r();
762                arg = get_locint_r();
763
764                ASSIGN_A(Ares, rpp_A[res])
765                ASSIGN_A(Aarg, rpp_A[arg]);
766
767                FOR_0_LE_l_LT_p
768                if  (0 == ARES) {
769                    HOV_INC(Ares, k1)
770                    HOV_INC(Aarg, k1)
771                } else {
772                    MAXDEC(AARG,ARES);
773                    AARG_INC_O;
774                    ARES_INC_O;
775                    FOR_0_LE_i_LT_k
776                    AARG_INC += ARES_INC;
777                }
778
779                GET_TAYL(res,k,p)
780                break;
781
782                /*--------------------------------------------------------------------------*/
783            case eq_min_d:       /* Subtract a floating point from an    eq_min_d */
784                /* adouble. (-=) */
785                res   = get_locint_r();
786                coval = get_val_r();
787
788                GET_TAYL(res,k,p)
789                break;
790
791                /*--------------------------------------------------------------------------*/
792            case eq_min_a:        /* Subtract an adouble from another    eq_min_a */
793                /* adouble. (-=) */
794                res = get_locint_r();
795                arg = get_locint_r();
796
797                ASSIGN_A(Ares, rpp_A[res])
798                ASSIGN_A(Aarg, rpp_A[arg])
799
800                FOR_0_LE_l_LT_p
801                if  (0==ARES) {
802                    HOV_INC(Ares, k1)
803                    HOV_INC(Aarg, k1)
804                } else {
805                    MAXDEC(AARG,ARES);
806                    AARG_INC_O;
807                    ARES_INC_O;
808                    FOR_0_LE_i_LT_k
809                    AARG_INC -= ARES_INC;
810                }
811
812                GET_TAYL(res,k,p)
813                break;
814
815                /*--------------------------------------------------------------------------*/
816            case eq_mult_d:              /* Multiply an adouble by a    eq_mult_d */
817                /* flaoting point. (*=) */
818                res   = get_locint_r();
819                coval = get_val_r();
820
821                ASSIGN_A(Ares, rpp_A[res])
822
823                FOR_0_LE_l_LT_p
824                if ( 0 == ARES_INC )
825                    HOV_INC(Ares, k)
826                    else
827                        FOR_0_LE_i_LT_k
828                        ARES_INC *= coval;
829
830                GET_TAYL(res,k,p)
831                break;
832
833                /*--------------------------------------------------------------------------*/
834            case eq_mult_a:       /* Multiply one adouble by another    eq_mult_a */
835                /* (*=) */
836                res = get_locint_r();
837                arg = get_locint_r();
838
839                GET_TAYL(res,k,p)
840
841                ASSIGN_A(Ares, rpp_A[res])
842                ASSIGN_A(Aarg, rpp_A[arg])
843                ASSIGN_A(Aqo,  rp_Atemp)
844                ASSIGN_T(Tres, rpp_T[res])
845                ASSIGN_T(Targ, rpp_T[arg])
846
847                FOR_0_LE_l_LT_p {
848                    if (0 == ARES) {
849                    HOV_INC(Aarg, k1)
850                        HOV_INC(Ares, k1)
851                    } else {
852                        MAXDEC(ARES,2.0);
853                        MAXDEC(AARG,ARES);
854                        AARG_INC_O;
855                        ARES_INC_O;
856                        conv(k,Ares,Targ,rp_Atemp);
857                        if(arg != res) {
858                            inconv(k,Ares,Tres,Aarg);
859                            FOR_0_LE_i_LT_k
860                            ARES_INC = AQO_INC;
861                        } else
862                            FOR_0_LE_i_LT_k
863                            ARES_INC = 2.0 * AQO_INC;
864                        HOV_INC(Aarg,k)
865                        HOS_OV_INC(Tres,k)
866                        HOS_OV_INC(Targ,k)
867                        HOS_OV_ASSIGN_A(Aqo,  rp_Atemp)
868                    }
869            }
870                break;
871
872                /*--------------------------------------------------------------------------*/
873            case incr_a:                        /* Increment an adouble    incr_a */
874            case decr_a:                        /* Increment an adouble    decr_a */
875                res   = get_locint_r();
876
877                GET_TAYL(res,k,p)
878                break;
879
880
881                /****************************************************************************/
882                /*                                                        BINARY OPERATIONS */
883
884                /*--------------------------------------------------------------------------*/
885            case plus_a_a:                 /* : Add two adoubles. (+)    plus a_a */
886                res  = get_locint_r();
887                arg2 = get_locint_r();
888                arg1 = get_locint_r();
889
890                ASSIGN_A(Ares,  rpp_A[res])
891                ASSIGN_A(Aarg1, rpp_A[arg1])
892                ASSIGN_A(Aarg2, rpp_A[arg2])
893
894                FOR_0_LE_l_LT_p
895                if  (0 == ARES) {
896                    HOV_INC(Ares,  k1)
897                    HOV_INC(Aarg1, k1)
898                    HOV_INC(Aarg2, k1)
899                } else {
900                    aTmp = ARES;
901                    ARES_INC = 0.0;
902                    MAXDEC(AARG1,aTmp);
903                    MAXDEC(AARG2,aTmp);
904                    AARG2_INC_O;
905                    AARG1_INC_O;
906                    FOR_0_LE_i_LT_k
907                    { aTmp = ARES;
908                      ARES_INC = 0.0;
909                      AARG1_INC += aTmp;
910                      AARG2_INC += aTmp;
911                    }
912                }
913
914                GET_TAYL(res,k,p)
915                break;
916
917                /*--------------------------------------------------------------------------*/
918            case plus_d_a:             /* Add an adouble and a double    plus_d_a */
919                /* (+) */
920                res   = get_locint_r();
921                arg   = get_locint_r();
922                coval = get_val_r();
923
924                ASSIGN_A(Ares, rpp_A[res])
925                ASSIGN_A(Aarg, rpp_A[arg])
926
927                FOR_0_LE_l_LT_p
928                if  (0 == ARES) {
929                    HOV_INC(Ares, k1)
930                    HOV_INC(Aarg, k1)
931                } else {
932                    aTmp = ARES;
933                    ARES_INC = 0.0;
934                    MAXDEC(AARG,aTmp);
935                    AARG_INC_O;
936                    FOR_0_LE_i_LT_k
937                    { aTmp = ARES;
938                      ARES_INC = 0.0;
939                      AARG_INC += aTmp;
940                    }
941                }
942
943                GET_TAYL(res,k,p)
944                break;
945
946                /*--------------------------------------------------------------------------*/
947            case min_a_a:              /* Subtraction of two adoubles    min_a_a */
948                /* (-) */
949                res  = get_locint_r();
950                arg2 = get_locint_r();
951                arg1 = get_locint_r();
952
953                ASSIGN_A(Ares,  rpp_A[res])
954                ASSIGN_A(Aarg1, rpp_A[arg1])
955                ASSIGN_A(Aarg2, rpp_A[arg2])
956
957                FOR_0_LE_l_LT_p
958                if  (0 == ARES) {
959                    HOV_INC(Ares,  k1)
960                    HOV_INC(Aarg1, k1)
961                    HOV_INC(Aarg2, k1)
962                } else {
963                    aTmp = ARES;
964                    ARES_INC = 0.0;
965                    MAXDEC(AARG1,aTmp);
966                    MAXDEC(AARG2,aTmp);
967                    AARG2_INC_O;
968                    AARG1_INC_O;
969                    FOR_0_LE_i_LT_k
970                    { aTmp = ARES;
971                      ARES_INC = 0.0;
972                      AARG1_INC += aTmp;
973                      AARG2_INC -= aTmp;
974                    }
975                }
976
977                GET_TAYL(res,k,p)
978                break;
979
980                /*--------------------------------------------------------------------------*/
981            case min_d_a:                /* Subtract an adouble from a    min_d_a */
982                /* double (-) */
983                res   = get_locint_r();
984                arg   = get_locint_r();
985                coval = get_val_r();
986
987                ASSIGN_A(Ares, rpp_A[res])
988                ASSIGN_A(Aarg, rpp_A[arg])
989
990                FOR_0_LE_l_LT_p
991                if (0 == ARES) {
992                    HOV_INC(Ares, k1)
993                    HOV_INC(Aarg, k1)
994                } else {
995                    aTmp = ARES;
996                    ARES_INC = 0.0;
997                    MAXDEC(AARG,aTmp);
998                    AARG_INC_O;
999                    FOR_0_LE_i_LT_k
1000                    { aTmp = ARES;
1001                      ARES_INC = 0.0;
1002                      AARG_INC -= aTmp;
1003                    }
1004                }
1005
1006                GET_TAYL(res,k,p)
1007                break;
1008
1009                /*--------------------------------------------------------------------------*/
1010            case mult_a_a:               /* Multiply two adoubles (*)    mult_a_a */
1011                res  = get_locint_r();
1012                arg2 = get_locint_r();
1013                arg1 = get_locint_r();
1014
1015                GET_TAYL(res,k,p)
1016
1017                ASSIGN_A(Ares,  rpp_A[res])
1018                ASSIGN_A(Aarg2, rpp_A[arg2])
1019                ASSIGN_A(Aarg1, rpp_A[arg1])
1020                ASSIGN_T(Targ1, rpp_T[arg1])
1021                ASSIGN_T(Targ2, rpp_T[arg2])
1022
1023                FOR_0_LE_l_LT_p
1024                if (0 == ARES) {
1025                    HOV_INC(Aarg1, k1)
1026                    HOV_INC(Aarg2, k1)
1027                    HOV_INC(Ares,  k1)
1028                } else {
1029                    comp = (ARES > 2.0) ? ARES : 2.0 ;
1030                    ARES_INC = 0.0;
1031                    MAXDEC(AARG1,comp);
1032                    MAXDEC(AARG2,comp);
1033                    AARG1_INC_O;
1034                    AARG2_INC_O;
1035
1036                    copyAndZeroset(k,Ares,rp_Atemp);
1037                    inconv(k,rp_Atemp,Targ1,Aarg2);
1038                    inconv(k,rp_Atemp,Targ2,Aarg1);
1039
1040                    HOV_INC(Ares,  k)
1041                    HOV_INC(Aarg1, k)
1042                    HOV_INC(Aarg2, k)
1043                    HOS_OV_INC(Targ1, k)
1044                    HOS_OV_INC(Targ2, k)
1045                }
1046                break;
1047
1048                /*--------------------------------------------------------------------------*/
1049                /* olvo 991122: new op_code with recomputation */
1050            case eq_plus_prod:   /* increment a product of           eq_plus_prod */
1051                /* two adoubles (*) */
1052                res  = get_locint_r();
1053                arg2 = get_locint_r();
1054                arg1 = get_locint_r();
1055
1056
1057                ASSIGN_A(Ares,  rpp_A[res])
1058                ASSIGN_A(Aarg2, rpp_A[arg2])
1059                ASSIGN_A(Aarg1, rpp_A[arg1])
1060                ASSIGN_T(Targ1, rpp_T[arg1])
1061                ASSIGN_T(Targ2, rpp_T[arg2])
1062
1063                /* RECOMPUTATION */
1064                ASSIGN_T( Tres,  rpp_T[res])
1065                deconv1(k,Targ1,Targ2,Tres);
1066
1067                FOR_0_LE_l_LT_p
1068                if (0 == ARES) {
1069                    HOV_INC(Aarg1, k1)
1070                    HOV_INC(Aarg2, k1)
1071                    HOV_INC(Ares,  k1)
1072                } else {
1073                    comp = (ARES > 2.0) ? ARES : 2.0 ;
1074                    ARES_INC = comp;
1075                    MAXDEC(AARG1,comp);
1076                    MAXDEC(AARG2,comp);
1077                    AARG1_INC_O;
1078                    AARG2_INC_O;
1079
1080                    inconv(k,Ares,Targ1,Aarg2);
1081                    inconv(k,Ares,Targ2,Aarg1);
1082
1083                    HOV_INC(Ares,  k)
1084                    HOV_INC(Aarg1, k)
1085                    HOV_INC(Aarg2, k)
1086                    HOS_OV_INC(Targ1, k)
1087                    HOS_OV_INC(Targ2, k)
1088                }
1089                break;
1090
1091                /*--------------------------------------------------------------------------*/
1092                /* olvo 991122: new op_code with recomputation */
1093            case eq_min_prod:   /* decrement a product of             eq_min_prod */
1094                /* two adoubles (*) */
1095                res  = get_locint_r();
1096                arg2 = get_locint_r();
1097                arg1 = get_locint_r();
1098
1099
1100                ASSIGN_A(Ares,  rpp_A[res])
1101                ASSIGN_A(Aarg2, rpp_A[arg2])
1102                ASSIGN_A(Aarg1, rpp_A[arg1])
1103                ASSIGN_T(Targ1, rpp_T[arg1])
1104                ASSIGN_T(Targ2, rpp_T[arg2])
1105
1106                /* RECOMPUTATION */
1107                ASSIGN_T( Tres,  rpp_T[res])
1108                inconv1(k,Targ1,Targ2,Tres);
1109
1110                FOR_0_LE_l_LT_p
1111                if (0 == ARES) {
1112                    HOV_INC(Aarg1, k1)
1113                    HOV_INC(Aarg2, k1)
1114                    HOV_INC(Ares,  k1)
1115                } else {
1116                    comp = (ARES > 2.0) ? ARES : 2.0 ;
1117                    ARES_INC = comp;
1118                    MAXDEC(AARG1,comp);
1119                    MAXDEC(AARG2,comp);
1120                    AARG1_INC_O;
1121                    AARG2_INC_O;
1122
1123                    deconv1(k,Ares,Targ1,Aarg2);
1124                    deconv1(k,Ares,Targ2,Aarg1);
1125
1126                    HOV_INC(Ares,  k)
1127                    HOV_INC(Aarg1, k)
1128                    HOV_INC(Aarg2, k)
1129                    HOS_OV_INC(Targ1, k)
1130                    HOS_OV_INC(Targ2, k)
1131                }
1132                break;
1133
1134                /*--------------------------------------------------------------------------*/
1135            case mult_d_a:         /* Multiply an adouble by a double    mult_d_a */
1136                /* (*) */
1137                res   = get_locint_r();
1138                arg   = get_locint_r();
1139                coval = get_val_r();
1140
1141                ASSIGN_A(Ares, rpp_A[res])
1142                ASSIGN_A(Aarg, rpp_A[arg])
1143
1144                FOR_0_LE_l_LT_p
1145                if (0 == ARES) {
1146                    HOV_INC(Ares, k1)
1147                    HOV_INC(Aarg, k1)
1148                } else {
1149                    aTmp = ARES;
1150                    ARES_INC = 0.0;
1151                    MAXDEC(AARG,aTmp);
1152                    AARG_INC_O;
1153                    FOR_0_LE_i_LT_k
1154                    { aTmp = ARES;
1155                      ARES_INC = 0.0;
1156                      AARG_INC += coval * aTmp;
1157                    }
1158                }
1159
1160                GET_TAYL(res,k,p)
1161                break;
1162
1163                /*--------------------------------------------------------------------------*/
1164            case div_a_a:           /* Divide an adouble by an adouble    div_a_a */
1165                /* (/) */
1166                res  = get_locint_r();
1167                arg2 = get_locint_r();
1168                arg1 = get_locint_r();
1169
1170                ASSIGN_A(Ares,  rpp_A[res])
1171                ASSIGN_A(Aarg2, rpp_A[arg2])
1172                ASSIGN_A(Aarg1, rpp_A[arg1])
1173                ASSIGN_T(Tres,  rpp_T[res])
1174                ASSIGN_T(Targ2, rpp_T[arg2])
1175
1176                /* olvo 980922 allows reflexive operation */
1177                if (arg2 == res) {
1178                    FOR_0_LE_l_LT_pk
1179                    rp_Ttemp2[l] = Tres[l];
1180                    Tres = rp_Ttemp2;
1181                    GET_TAYL(res,k,p)
1182                }
1183
1184                VEC_COMPUTED_INIT
1185                FOR_0_LE_l_LT_p
1186                { if (0 == ARES) {
1187                  HOV_INC(Ares,  k1)
1188                      HOV_INC(Aarg1, k1)
1189                      HOV_INC(Aarg2, k1)
1190                  } else {
1191                      aTmp = ARES;
1192                      ARES_INC = 0.0;
1193                      MAXDEC(AARG1,3.0);
1194                      MAXDEC(AARG1,aTmp);
1195                      MAXDEC(AARG2,3.0);
1196                      MAXDEC(AARG2,aTmp);
1197                      AARG1_INC_O;
1198                      AARG2_INC_O;
1199
1200                      VEC_COMPUTED_CHECK
1201                      recipr(k,1.0,Targ2,rp_Ttemp);
1202                      conv0(k ,rp_Ttemp,
1203                           Tres, rp_Atemp2);
1204                      VEC_COMPUTED_END
1205                      copyAndZeroset(k,Ares,rp_Atemp);
1206                      inconv(k, rp_Atemp,
1207                             rp_Ttemp, Aarg1);
1208                      deconv(k, rp_Atemp,
1209                             rp_Atemp2, Aarg2);
1210
1211                      HOV_INC(Ares,  k)
1212                      HOV_INC(Aarg1, k)
1213                      HOV_INC(Aarg2, k)
1214                      HOS_OV_INC(Tres, k)
1215                      HOS_OV_INC(Targ2, k)
1216                  }
1217            }
1218
1219                if (res != arg2)
1220                    GET_TAYL(res,k,p)
1221                    break;
1222
1223                /*--------------------------------------------------------------------------*/
1224            case div_d_a:             /* Division double - adouble (/)    div_d_a */
1225                res   = get_locint_r();
1226                arg   = get_locint_r();
1227                coval = get_val_r();
1228
1229                ASSIGN_A(Ares, rpp_A[res])
1230                ASSIGN_A(Aarg, rpp_A[arg])
1231                ASSIGN_T(Tres, rpp_T[res])
1232                ASSIGN_T(Targ, rpp_T[arg])
1233
1234                /* olvo 980922 allows reflexive operation */
1235                if (arg == res) {
1236                    FOR_0_LE_l_LT_pk
1237                    rp_Ttemp2[l] = Tres[l];
1238                    Tres = rp_Ttemp2;
1239                    GET_TAYL(arg,k,p)
1240                }
1241
1242                VEC_COMPUTED_INIT
1243                FOR_0_LE_l_LT_p
1244                { if (0 == ARES) {
1245                  HOV_INC(Ares, k1)
1246                      HOV_INC(Aarg, k1)
1247                  } else {
1248                      aTmp = ARES;
1249                      ARES_INC = 0.0;
1250                      MAXDEC(AARG,aTmp);
1251                      MAXDEC(AARG,3.0);
1252                      AARG_INC_O;
1253
1254                      VEC_COMPUTED_CHECK
1255                      recipr(k,1.0,Targ,rp_Ttemp);
1256                      conv0(k, rp_Ttemp,
1257                           Tres, rp_Atemp);
1258                      VEC_COMPUTED_END
1259                      deconvZeroR(k,Ares,rp_Atemp,Aarg);
1260
1261                      HOV_INC(Ares, k)
1262                      HOV_INC(Aarg, k)
1263                      HOS_OV_INC(Tres, k)
1264                      HOS_OV_INC(Targ, k)
1265                  }
1266            }
1267
1268                GET_TAYL(res,k,p)
1269                break;
1270
1271
1272                /****************************************************************************/
1273                /*                                                         SIGN  OPERATIONS */
1274
1275                /*--------------------------------------------------------------------------*/
1276            case pos_sign_a:                                        /* pos_sign_a */
1277                res   = get_locint_r();
1278                arg   = get_locint_r();
1279
1280                ASSIGN_A(Ares, rpp_A[res])
1281                ASSIGN_A(Aarg, rpp_A[arg])
1282
1283                FOR_0_LE_l_LT_p
1284                if  (0 == ARES) {
1285                    HOV_INC(Ares, k1)
1286                    HOV_INC(Aarg, k1)
1287                } else {
1288                    aTmp = ARES;
1289                    ARES_INC = 0.0;
1290                    MAXDEC(AARG,aTmp);
1291                    AARG_INC_O;
1292                    FOR_0_LE_i_LT_k
1293                    { aTmp = ARES;
1294                      ARES_INC = 0.0;
1295                      AARG_INC += aTmp;
1296                    }
1297                }
1298
1299                GET_TAYL(res,k,p)
1300                break;
1301
1302                /*--------------------------------------------------------------------------*/
1303            case neg_sign_a:                                        /* neg_sign_a */
1304                res   = get_locint_r();
1305                arg   = get_locint_r();
1306
1307                ASSIGN_A(Ares, rpp_A[res])
1308                ASSIGN_A(Aarg, rpp_A[arg])
1309
1310                FOR_0_LE_l_LT_p
1311                if  (0 == ARES) {
1312                    HOV_INC(Ares, k1)
1313                    HOV_INC(Aarg, k1)
1314                } else {
1315                    aTmp = ARES;
1316                    ARES_INC = 0.0;
1317                    MAXDEC(AARG,aTmp);
1318                    AARG_INC_O;
1319                    FOR_0_LE_i_LT_k
1320                    { aTmp = ARES;
1321                      ARES_INC = 0.0;
1322                      AARG_INC -= aTmp;
1323                    }
1324                }
1325
1326                GET_TAYL(res,k,p)
1327                break;
1328
1329
1330                /****************************************************************************/
1331                /*                                                         UNARY OPERATIONS */
1332
1333                /*--------------------------------------------------------------------------*/
1334            case exp_op:                          /* exponent operation    exp_op */
1335                res = get_locint_r();
1336                arg = get_locint_r();
1337
1338                ASSIGN_A(Ares, rpp_A[res])
1339                ASSIGN_A(Aarg, rpp_A[arg])
1340                ASSIGN_T(Tres, rpp_T[res])
1341                ASSIGN_T(Targ, rpp_T[arg])
1342
1343                FOR_0_LE_l_LT_p
1344                { if (0 == ARES) {
1345                  HOV_INC(Aarg, k1)
1346                      HOV_INC(Ares, k1)
1347                  } else {
1348                      aTmp = ARES;
1349                      ARES_INC = 0.0;
1350                      MAXDEC(AARG,aTmp);
1351                      MAXDEC(AARG,4.0);
1352                      AARG_INC_O;
1353
1354                      inconv0(k,Ares,Tres,Aarg);
1355
1356                      HOV_INC(Ares, k)
1357                      HOV_INC(Aarg, k)
1358                      HOS_OV_INC(Tres, k)
1359                  }
1360            }
1361
1362                GET_TAYL(res,k,p)
1363                break;
1364
1365                /*--------------------------------------------------------------------------*/
1366            case sin_op:                              /* sine operation    sin_op */
1367                res  = get_locint_r();
1368                arg2 = get_locint_r();
1369                arg1 = get_locint_r();
1370
1371                ASSIGN_A(Ares,  rpp_A[res])
1372                ASSIGN_A(Aarg1, rpp_A[arg1])
1373                ASSIGN_T(Targ2, rpp_T[arg2])
1374
1375                FOR_0_LE_l_LT_p
1376                { if (0 == ARES) {
1377                  HOV_INC(Aarg1, k1)
1378                      HOV_INC(Ares,  k1)
1379                  } else {
1380                      aTmp = ARES;
1381                      ARES_INC = 0.0;
1382                      MAXDEC(AARG1,aTmp);
1383                      MAXDEC(AARG1,4.0);
1384                      AARG1_INC_O;
1385
1386                      inconv0(k,Ares,Targ2,Aarg1);
1387
1388                      HOV_INC(Ares,  k)
1389                      HOV_INC(Aarg1, k)
1390                      HOS_OV_INC(Targ2, k)
1391                  }
1392            }
1393
1394                GET_TAYL(res,k,p)
1395                GET_TAYL(arg2,k,p) /* olvo 980710 covalue */
1396                /* NOTE: rpp_A[arg2] should be 0 already */
1397                break;
1398
1399                /*--------------------------------------------------------------------------*/
1400            case cos_op:                            /* cosine operation    cos_op */
1401                res  = get_locint_r();
1402                arg2 = get_locint_r();
1403                arg1 = get_locint_r();
1404
1405                ASSIGN_A(Ares,  rpp_A[res])
1406                ASSIGN_A(Aarg1, rpp_A[arg1])
1407                ASSIGN_T(Targ2, rpp_T[arg2])
1408
1409                FOR_0_LE_l_LT_p
1410                { if (0 == ARES) {
1411                  HOV_INC(Aarg1, k1)
1412                      HOV_INC(Ares,  k1)
1413                  } else {
1414                      aTmp = ARES;
1415                      ARES_INC = 0.0;
1416                      MAXDEC(AARG1,aTmp);
1417                      MAXDEC(AARG1,4.0);
1418                      AARG1_INC_O;
1419
1420                      deconv0(k,Ares,Targ2,Aarg1);
1421
1422                      HOV_INC(Ares,  k)
1423                      HOV_INC(Aarg1, k)
1424                      HOS_OV_INC(Targ2, k)
1425                  }
1426            }
1427
1428                GET_TAYL(res,k,p)
1429                GET_TAYL(arg2,k,p) /* olvo 980710 covalue */
1430                /* NOTE: rpp_A[arg2] should be 0 already */
1431                break;
1432                /*xxx*/
1433                /*--------------------------------------------------------------------------*/
1434            case atan_op:                                             /* atan_op  */
1435            case asin_op:                                             /* asin_op  */
1436            case acos_op:                                             /* acos_op  */
1437            case asinh_op:                                            /* asinh_op */
1438            case acosh_op:                                            /* acosh_op */
1439            case atanh_op:                                            /* atanh_op */
1440            case erf_op:                                              /* erf_op   */
1441                res  = get_locint_r();
1442                arg2 = get_locint_r();
1443                arg1 = get_locint_r();
1444
1445                GET_TAYL(res,k,p)
1446
1447                ASSIGN_A(Ares,  rpp_A[res])
1448                ASSIGN_A(Aarg1, rpp_A[arg1])
1449                ASSIGN_T(Targ2, rpp_T[arg2])
1450
1451                FOR_0_LE_l_LT_p
1452                { if (0 == ARES) {
1453                  HOV_INC(Aarg1, k1)
1454                      HOV_INC(Ares,  k1)
1455                  } else {
1456                      aTmp = ARES;
1457                      ARES_INC = 0.0;
1458                      MAXDEC(AARG1,aTmp);
1459                      MAXDEC(AARG1,4.0);
1460                      AARG1_INC_O;
1461
1462                      inconv0(k,Ares,Targ2,Aarg1);
1463
1464                      HOV_INC(Aarg1, k)
1465                      HOV_INC(Ares,  k)
1466                      HOS_OV_INC(Targ2, k)
1467                  }
1468            }
1469                break;
1470
1471                /*--------------------------------------------------------------------------*/
1472            case log_op:                                                /* log_op */
1473                res = get_locint_r();
1474                arg = get_locint_r();
1475
1476                GET_TAYL(res,k,p)
1477
1478                ASSIGN_A(Ares, rpp_A[res])
1479                ASSIGN_A(Aarg, rpp_A[arg])
1480                ASSIGN_T(Targ, rpp_T[arg])
1481
1482                VEC_COMPUTED_INIT
1483                FOR_0_LE_l_LT_p
1484                { if (0 == ARES) {
1485                  HOV_INC(Aarg, k1)
1486                      HOV_INC(Ares, k1)
1487                  } else {
1488                      aTmp = ARES;
1489                      ARES_INC = 0.0;
1490                      MAXDEC(AARG,aTmp);
1491                      MAXDEC(AARG,4.0);
1492                      AARG_INC_O;
1493
1494                      VEC_COMPUTED_CHECK
1495                      recipr(k,1.0,Targ,rp_Ttemp);
1496                      VEC_COMPUTED_END
1497                      inconv0(k,Ares,rp_Ttemp,Aarg);
1498
1499                      HOV_INC(Ares, k)
1500                      HOV_INC(Aarg, k)
1501                      HOS_OV_INC(Targ2, k)
1502                  }
1503            }
1504                break;
1505
1506                /*--------------------------------------------------------------------------*/
1507            case pow_op:                                                /* pow_op */
1508                res   = get_locint_r();
1509                arg   = get_locint_r();
1510                coval = get_val_r();
1511
1512                ASSIGN_T(Targ, rpp_T[arg])
1513                ASSIGN_T(Tres, rpp_T[res])
1514                ASSIGN_A(Ares, rpp_A[res])
1515                ASSIGN_A(Aarg, rpp_A[arg])
1516
1517                /* olvo 980921 allows reflexive operation */
1518                if (arg == res) {
1519                    FOR_0_LE_l_LT_pk
1520                    rp_Ttemp2[l] = Tres[l];
1521                    Tres = rp_Ttemp2;
1522                    GET_TAYL(arg,k,p)
1523                }
1524
1525                VEC_COMPUTED_INIT
1526                FOR_0_LE_l_LT_p
1527                if (0 == ARES) {
1528                    HOV_INC(Aarg, k1)
1529                    HOV_INC(Ares, k1)
1530                } else {
1531                    aTmp = ARES;
1532                    ARES_INC = 0.0;
1533                    MAXDEC(AARG,aTmp);
1534                    MAXDEC(AARG,4.0);
1535                    AARG_INC_O;
1536
1537                    VEC_COMPUTED_CHECK
1538                    if (fabs(Targ[0]) > ADOLC_EPS) {
1539                        divide(k,Tres,Targ,rp_Ttemp);
1540                        for (i=0;i<k;i++) {
1541                            rp_Ttemp[i] *= coval;
1542                            /*                 printf(" EPS i %d %f\n",i,rp_Ttemp[i]); */
1543                        }
1544                        inconv0(k,Ares,rp_Ttemp,Aarg);
1545                    } else {
1546                        if (coval <= 0.0) {
1547                            FOR_0_LE_i_LT_k
1548                            {
1549                                Aarg[i] = make_nan();
1550                                Ares[i] = 0;
1551                            }
1552                        } else {
1553                            /* coval not a whole number */
1554                            if (coval - floor(coval) != 0) {
1555                                i = 0;
1556                                FOR_0_LE_i_LT_k
1557                                {
1558                                    if (coval - i > 1) {
1559                                    Aarg[i] = 0;
1560                                        Ares[i] = 0;
1561                                    }
1562                                    if ((coval - i < 1) && (coval - i > 0)) {
1563                                    Aarg[i] = make_inf();
1564                                        Ares[i] = 0;
1565                                    }
1566                                    if (coval - i < 0) {
1567                                    Aarg[i] = make_nan();
1568                                        Ares[i] = 0;
1569                                    }
1570                                }
1571                            } else {
1572                                if (coval == 1) {
1573                                    FOR_0_LE_i_LT_k
1574                                    { /* ! no tempory */
1575                                        AARG_INC += ARES;
1576                                        ARES_INC = 0.0;
1577                                    }
1578                                } else {
1579                                    /* coval is an int > 1 */
1580                                    /* the following is not efficient but at least it works */
1581                                    /* it reformulates x^n into x* ... *x n times */
1582
1583                                    copyAndZeroset(k,Ares,rp_Atemp);
1584                                    inconv(k,rp_Atemp,Targ,Aarg);
1585                                    inconv(k,rp_Atemp,Targ,Aarg);
1586                                    if (coval == 3) {
1587                                        conv(k,Aarg,Targ,rp_Atemp);
1588                                        FOR_0_LE_i_LT_k
1589                                        Aarg[i] = 2.0 * rp_Atemp[i];
1590                                   }
1591                                }
1592                            }
1593                        }
1594                    }
1595                    VEC_COMPUTED_END
1596
1597                    HOV_INC(Ares, k)
1598                    HOV_INC(Aarg, k)
1599                    HOS_OV_INC(Tres, k)
1600                    HOS_OV_INC(Targ, k)
1601                }
1602
1603                GET_TAYL(res,k,p)
1604                break;
1605
1606                /*--------------------------------------------------------------------------*/
1607            case sqrt_op:                                              /* sqrt_op */
1608                res = get_locint_r();
1609                arg = get_locint_r();
1610
1611                ASSIGN_A(Ares, rpp_A[res])
1612                ASSIGN_A(Aarg, rpp_A[arg])
1613                ASSIGN_T(Tres, rpp_T[res])
1614
1615                VEC_COMPUTED_INIT
1616                FOR_0_LE_l_LT_p
1617                if (0 == ARES) {
1618                    HOV_INC(Aarg, k1)
1619                    HOV_INC(Ares, k1)
1620                } else {
1621                    aTmp = ARES;
1622                    ARES_INC = 0.0;
1623                    MAXDEC(AARG,aTmp);
1624                    MAXDEC(AARG,4.0);
1625                    AARG_INC_O;
1626
1627                    VEC_COMPUTED_CHECK
1628                    recipr(k,0.5,Tres,rp_Ttemp);
1629                    VEC_COMPUTED_END
1630                    inconv0(k,Ares,rp_Ttemp,Aarg);
1631
1632                    HOV_INC(Ares, k)
1633                    HOV_INC(Aarg, k)
1634                    HOS_OV_INC(Tres,k)
1635                }
1636
1637                GET_TAYL(res,k,p)
1638                break;
1639
1640                /*--------------------------------------------------------------------------*/
1641            case gen_quad:                                            /* gen_quad */
1642                res   = get_locint_r();
1643                arg2  = get_locint_r();
1644                arg1  = get_locint_r();
1645                coval = get_val_r();
1646                coval = get_val_r();
1647
1648                ASSIGN_A(Ares,  rpp_A[res])
1649                ASSIGN_A(Aarg1, rpp_A[arg1])
1650                ASSIGN_T(Targ2, rpp_T[arg2])
1651
1652                FOR_0_LE_l_LT_p
1653                if (0 == ARES) {
1654                    HOV_INC(Aarg1, k1)
1655                    HOV_INC(Ares,  k1)
1656                } else {
1657                    aTmp = ARES;
1658                    ARES_INC = 0.0;
1659                    MAXDEC(AARG1,aTmp);
1660                    MAXDEC(AARG1,4.0);
1661                    AARG1_INC_O;
1662
1663                    inconv0(k,Ares,Targ2,Aarg1);
1664
1665                    HOV_INC(Aarg1, k)
1666                    HOV_INC(Ares,  k)
1667                    HOS_OV_INC(Targ2,  k)
1668                }
1669
1670                GET_TAYL(res,k,p)
1671                break;
1672
1673                /*--------------------------------------------------------------------------*/
1674            case min_op:                                                /* min_op */
1675
1676#ifdef _HOS_OV_
1677
1678                fprintf(DIAG_OUT," operation min_op not implemented for hos_ov");
1679                break;
1680#endif
1681                res   = get_locint_r();
1682                arg2  = get_locint_r();
1683                arg1  = get_locint_r();
1684                coval = get_val_r();
1685
1686                GET_TAYL(res,k,p)
1687
1688                ASSIGN_A(Aarg1, rpp_A[arg1])
1689                ASSIGN_A(Aarg2, rpp_A[arg2])
1690                ASSIGN_A(Ares,  rpp_A[res])
1691                ASSIGN_T(Targ1, rpp_T[arg1])
1692                ASSIGN_T(Targ2, rpp_T[arg2])
1693                ASSIGN_A(AP1,   NULL)
1694                ASSIGN_A(AP2,   Ares)
1695
1696                if (Targ1[0] > Targ2[0]) {
1697                    FOR_0_LE_l_LT_p
1698                    { if ((coval) && (*AP2))
1699                      MINDEC(ret_c,2);
1700                      HOV_INC(AP2,k1)
1701                    }
1702                    AP1 = Aarg2;
1703                arg = 0;
1704            } else
1705                if (Targ1[0] < Targ2[0]) {
1706                        FOR_0_LE_l_LT_p
1707                        { if ((!coval) && (*AP2))
1708                          MINDEC(ret_c,2);
1709                          HOV_INC(AP2,k1)
1710                        }
1711                        AP1 = Aarg1;
1712                    arg = 0;
1713                } else /* both are equal */ /* must be changed for hos_ov, but how? */
1714                    /* seems to influence the return value */
1715                    for (i=1;i<k;i++) {
1716                            if (Targ1[i] > Targ2[i]) {
1717                                FOR_0_LE_l_LT_p
1718                                { if (*AP2)
1719                                  MINDEC(ret_c,1);
1720                                  HOV_INC(AP2,k1)
1721                                }
1722                                AP1 = Aarg2;
1723                            arg = i+1;
1724                        } else
1725                            if (Targ1[i] < Targ2[i]) {
1726                                    FOR_0_LE_l_LT_p
1727                                    { if (*AP2)
1728                                      MINDEC(ret_c,1);
1729                                      HOV_INC(AP2,k1)
1730                                    }
1731                                    AP1 = Aarg1;
1732                                arg = i+1;
1733                            }
1734                        if (AP1 != NULL)
1735                                break;
1736                        }
1737
1738                if (AP1 != NULL)
1739                    FOR_0_LE_l_LT_p
1740                    { if (0 == ARES) {
1741                      HOV_INC(AP1, k1)
1742                          HOV_INC(Ares,k1);
1743                      } else {
1744                          aTmp = ARES;
1745                          ARES_INC = 0.0;
1746                          if (arg)  /* we are at the tie */
1747                              *AP1 = 5.0;
1748                          else
1749                              MAXDEC(*AP1,aTmp);
1750                          AP1++;
1751                          for (i=0;i<k;i++) {
1752                              aTmp = ARES;
1753                              ARES_INC = 0.0;
1754                              *AP1++ += aTmp;
1755                          }
1756                      }
1757            }
1758                else /* both are identical */
1759                {
1760                    FOR_0_LE_l_LT_p
1761                    { if (0 == ARES) {
1762                      HOV_INC(Aarg1,k1)
1763                          HOV_INC(Aarg2,k1)
1764                          HOV_INC(Ares, k1)
1765                      } else {
1766                          aTmp = ARES;
1767                          ARES_INC = 0.0;
1768                          MAXDEC(AARG1,aTmp);  /*assume sthg like fmin(x,x) */
1769                          MAXDEC(AARG2,aTmp);
1770                          AARG1_INC_O;
1771                          AARG2_INC_O;
1772                          for (i=0;i<k;i++) {
1773                              aTmp = ARES;
1774                              ARES_INC = 0.0;
1775                              AARG1_INC += aTmp/2;
1776                              AARG2_INC += aTmp/2;
1777                          }
1778                      }
1779                    }
1780                    if (arg1 != arg2)
1781                        MINDEC(ret_c,1);
1782                }
1783                break;
1784
1785
1786                /*--------------------------------------------------------------------------*/
1787            case abs_val:                                              /* abs_val */
1788                res   = get_locint_r();
1789                arg   = get_locint_r();
1790                coval = get_val_r();
1791                /* must be changed for hos_ov, but how? */
1792                /* seems to influence the return value  */
1793                GET_TAYL(res,k,p)
1794
1795                ASSIGN_A(Ares, rpp_A[res])
1796                ASSIGN_A(Aarg, rpp_A[arg])
1797                ASSIGN_T(Targ, rpp_T[arg])
1798
1799                FOR_0_LE_l_LT_q
1800                {
1801                    x[l] = 0.0;
1802                    jj[l] = 0;
1803                    for (i=0;i<k;i++)
1804                    if ( (x[l] == 0.0) && (Targ[i] != 0.0) ) {
1805                        jj[l] = i;
1806                            if (Targ[i] < 0.0)
1807                                x[l] = -1.0;
1808                            else
1809                                x[l] = 1.0;
1810                        }
1811                    HOS_OV_INC(Targ,k)
1812            }
1813                ASSIGN_T(Targ, rpp_T[arg])
1814                FOR_0_LE_l_LT_p
1815                { if (0 == ARES) {
1816                  HOV_INC(Aarg, k1)
1817                      HOV_INC(Ares, k1)
1818                  } else {
1819                      if (Targ[0] == 0.0) {
1820                          ARES_INC = 0.0;
1821                          AARG_INC = 5.0;
1822                      } else {
1823                          aTmp = ARES;
1824                          ARES_INC = 0.0;
1825                          MAXDEC(AARG,aTmp);
1826                          AARG_INC_O;
1827                      }
1828                      if(Targ[0] == 0.0)
1829                          MINDEC(ret_c,1);
1830                      for (i=0;i<jj[l];i++)
1831                          ARES_INC = 0.0;
1832                      Aarg += jj[l];
1833                      for (i=jj[l];i<k;i++) {
1834                          aTmp = ARES;
1835                          ARES_INC = 0.0;
1836                          if ( (coval) && (x[l]<0) && (aTmp) )
1837                              MINDEC(ret_c,2);
1838                          if ( (!coval) && (x[l]>0) && (aTmp))
1839                              MINDEC(ret_c,2);
1840                          AARG_INC += x[l] * aTmp;
1841                      }
1842                  }
1843                  HOS_OV_INC(Targ,k)
1844            }
1845                break;
1846
1847                /*--------------------------------------------------------------------------*/
1848            case ceil_op:                                              /* ceil_op */
1849                res   = get_locint_r();
1850                arg   = get_locint_r();
1851                coval = get_val_r();
1852
1853                GET_TAYL(res,k,p)
1854
1855                coval = (coval != ceil(*rpp_T[arg]) );
1856
1857                ASSIGN_A(Ares, rpp_A[res])
1858
1859                FOR_0_LE_l_LT_p
1860                if (0 == ARES) {
1861                    HOV_INC(Aarg,  k1)
1862                    HOV_INC(Ares,  k1)
1863                } else {
1864                    ARES_INC = 0.0;
1865                    AARG_INC = 5.0;
1866                    FOR_0_LE_i_LT_k
1867                    { if ((coval) && (ARES))
1868                      MINDEC(ret_c,2);
1869                      ARES_INC = 0.0;
1870                    }
1871                    HOV_INC(Aarg, k)
1872                    }
1873                break;
1874
1875                /*--------------------------------------------------------------------------*/
1876            case floor_op:                                            /* floor_op */
1877                res   = get_locint_r();
1878                arg   = get_locint_r();
1879                coval = get_val_r();
1880
1881                GET_TAYL(res,k,p)
1882
1883                coval = ( coval != floor(*rpp_T[arg]) );
1884
1885                ASSIGN_A(Ares, rpp_A[res])
1886                ASSIGN_A(Aarg, rpp_A[arg])
1887
1888                FOR_0_LE_l_LT_p
1889                if (0 == ARES) {
1890                    HOV_INC(Aarg, k1)
1891                    HOV_INC(Ares, k1)
1892                } else {
1893                    ARES = 0.0;
1894                    AARG_INC = 5.0;
1895                    FOR_0_LE_i_LT_k
1896                    { if ( (coval) && (ARES) )
1897                      MINDEC(ret_c,2);
1898                      ARES_INC = 0.0;
1899                    }
1900                    HOV_INC(Aarg, k)
1901                    }
1902                break;
1903
1904
1905                /****************************************************************************/
1906                /*                                                             CONDITIONALS */
1907
1908                /*--------------------------------------------------------------------------*/
1909            case cond_assign:                                      /* cond_assign */
1910                res   = get_locint_r();
1911                arg2  = get_locint_r();
1912                arg1  = get_locint_r();
1913                arg   = get_locint_r();
1914                coval = get_val_r();
1915
1916                GET_TAYL(res,k,p)
1917
1918                ASSIGN_A(Aarg1, rpp_A[arg1])
1919                ASSIGN_A(Ares,  rpp_A[res])
1920                ASSIGN_A(Aarg2, rpp_A[arg2])
1921                ASSIGN_T(Targ,  rpp_T[arg])
1922
1923                /* olvo 980925 changed code a little bit */
1924                if (TARG > 0.0) {
1925                    if (res != arg1)
1926                        FOR_0_LE_l_LT_p
1927                        { if (0 == ARES) {
1928                          HOV_INC(Ares,  k1)
1929                              HOV_INC(Aarg1, k1)
1930                          } else {
1931                              if (coval <= 0.0)
1932                                  MINDEC(ret_c,2);
1933                              MAXDEC(AARG1,ARES);
1934                              ARES_INC = 0.0;
1935                              AARG1_INC_O;
1936                              FOR_0_LE_i_LT_k
1937                              { AARG1_INC += ARES;
1938                                ARES_INC = 0;
1939                              }
1940                          }
1941                    }
1942                    else
1943                        FOR_0_LE_l_LT_p {
1944                            if ((coval <= 0.0) && (ARES))
1945                            MINDEC(ret_c,2);
1946                            HOV_INC(Ares,  k1)
1947                        }
1948                    } else /* TARG <= 0.0 */
1949            {
1950                if (res != arg2)
1951                        FOR_0_LE_l_LT_p
1952                        { if (0 == ARES) {
1953                          HOV_INC(Ares,  k1)
1954                              HOV_INC(Aarg2, k1)
1955                          } else {
1956                              if (TARG == 0.0) /* we are at the tie */
1957                              { MINDEC(ret_c,0);
1958                                  AARG1 = 5.0;
1959                                  AARG2_INC = 5.0;
1960                              } else {
1961                                  if (coval <= 0.0)
1962                                      MINDEC(ret_c,2);
1963                                  MAXDEC(AARG2,ARES);
1964                                  AARG2_INC_O;
1965                              }
1966                              ARES_INC = 0.0;
1967
1968                              FOR_0_LE_i_LT_k
1969                              { AARG2_INC += ARES;
1970                                ARES_INC = 0;
1971                              }
1972                          }
1973                      HOV_INC(Aarg1, k1)
1974                    } else
1975                        FOR_0_LE_l_LT_p {
1976                            if (ARES) {
1977                            if (TARG == 0.0) /* we are at the tie */
1978                                { MINDEC(ret_c,0);
1979                                    AARG1 = 5.0;
1980                                    AARG2 = 5.0;
1981                                } else
1982                                    if (coval <= 0.0)
1983                                        MINDEC(ret_c,2);
1984                            }
1985                        HOV_INC(Ares,  k1)
1986                        HOV_INC(Aarg1, k1)
1987                        HOV_INC(Aarg2, k1)
1988                    }
1989                }
1990                break;
1991
1992                /*--------------------------------------------------------------------------*/
1993            case cond_assign_s:                                  /* cond_assign_s */
1994                res   = get_locint_r();
1995                arg1  = get_locint_r();
1996                arg   = get_locint_r();
1997                coval = get_val_r();
1998
1999                GET_TAYL(res,k,p)
2000
2001                ASSIGN_A(Aarg1, rpp_A[arg1])
2002                ASSIGN_A(Ares,  rpp_A[res])
2003                ASSIGN_T(Targ,  rpp_T[arg])
2004
2005                /* olvo 980925 changed code a little bit */
2006                if (TARG == 0.0) /* we are at the tie */
2007                { FOR_0_LE_l_LT_p
2008                    { if  (ARES)
2009                      AARG1 = 5.0;
2010                      HOV_INC(Aarg1, k1)
2011                      HOV_INC(Ares,  k1)
2012                    }
2013                    MINDEC(ret_c,0);
2014                } else
2015                    if (TARG > 0.0) {
2016                        if (res != arg1)
2017                            FOR_0_LE_l_LT_p
2018                            { if  (0 == ARES) {
2019                              HOV_INC(Ares,  k1)
2020                                  HOV_INC(Aarg1, k1)
2021                              } else {
2022                                  if (coval <= 0.0)
2023                                      MINDEC(ret_c,2);
2024                                  MAXDEC(AARG1,ARES);
2025                                  ARES_INC = 0.0;
2026                                  AARG1_INC_O;
2027                                  FOR_0_LE_i_LT_k
2028                                  { (AARG1_INC) += ARES;
2029                                    ARES_INC = 0;
2030                                  }
2031                              }
2032                        }
2033                        else
2034                            FOR_0_LE_l_LT_p {
2035                                if ((coval <= 0.0) && (ARES))
2036                                MINDEC(ret_c,2);
2037                                HOV_INC(Ares,  k1)
2038                            }
2039                        }
2040            break;
2041
2042                /****************************************************************************/
2043                /*                                                          REMAINING STUFF */
2044
2045                /*--------------------------------------------------------------------------*/
2046            case take_stock_op:                                  /* take_stock_op */
2047                res = get_locint_r();
2048                size = get_locint_r();
2049                d = get_val_v_r(size);
2050
2051                res += size;
2052                for (ls=size;ls>0;ls--) {
2053                    res--;
2054
2055                    ASSIGN_A( Ares, rpp_A[res])
2056
2057                    FOR_0_LE_l_LT_pk1
2058                    ARES_INC = 0.0;
2059                }
2060                break;
2061
2062                /*--------------------------------------------------------------------------*/
2063            case death_not:                                          /* death_not */
2064                arg2 = get_locint_r();
2065                arg1 = get_locint_r();
2066
2067                for (j=arg1;j<=arg2;j++) {
2068                    ASSIGN_A(Aarg1, rpp_A[j])
2069
2070                    FOR_0_LE_l_LT_p
2071                    for (i=0; i<k1; i++)
2072                        AARG1_INC = 0.0;
2073                }
2074               
2075                for (j=arg1;j<=arg2;j++)
2076                    GET_TAYL(j,k,p)
2077
2078                break;
2079
2080                /*--------------------------------------------------------------------------*/
2081            default:                                                   /* default */
2082                /*             Die here, we screwed up     */
2083
2084                fprintf(DIAG_OUT,"ADOL-C fatal error in " GENERATED_FILENAME " ("
2085                        __FILE__
2086                        ") : no such operation %d\n", operation);
2087                exit(-1);
2088                break;
2089        } /* endswitch */
2090
2091        /* Get the next operation */   
2092        operation=get_op_r();
2093#if defined(ADOLC_DEBUG)
2094        ++countPerOperation[operation];
2095#endif /* ADOLC_DEBUG */
2096    }
2097
2098#if defined(ADOLC_DEBUG)
2099    printf("\nTape contains:\n");
2100    for (v = 0; v < 256; ++v)
2101        if (countPerOperation[v] > 0)
2102            printf("operation %3d: %6d time(s) - %6d taylors read (%10.2f per operation)\n", v, countPerOperation[v], taylorPerOperation[v], (double)taylorPerOperation[v] / (double)countPerOperation[v]);
2103    printf("\n");
2104#endif /* ADOLC_DEBUG */
2105
2106    /* clean up */
2107    free((char*)*rpp_T);
2108    free((char*) rpp_T);
2109    free(*rpp_A);
2110    free(rpp_A);
2111    free(rp_Ttemp);
2112    free(rp_Ttemp2);
2113    free(rp_Atemp);
2114    free(rp_Atemp2);
2115
2116    free((char*) jj);
2117    free((char*) x);
2118
2119    end_sweep();
2120
2121    return ret_c;
2122}
2123
2124/****************************************************************************/
2125/*                                                               THAT'S ALL */
2126
2127END_C_DECLS
Note: See TracBrowser for help on using the repository browser.