source: trunk/ADOL-C/src/uni5_for.c @ 360

Last change on this file since 360 was 360, checked in by kulshres, 7 years ago

implement traced logical operators for condassign and advector

this will help to handle code branching as long as condassign or
advector based generalized conditional assignment is used.

These operators are only used if both operands are badoubles. If the
old behaviour is required then at least one of the two operands must
be a double.

For if-then-else branches the old behaviour should continue and we can do
nothing.

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

  • Property svn:keywords set to Author Date Id Revision
File size: 153.0 KB
Line 
1
2/*----------------------------------------------------------------------------
3 ADOL-C -- Automatic Differentiation by Overloading in C++
4 File:     uni5_for.c
5
6
7 Revision: $Id: uni5_for.c 360 2012-10-19 09:22:43Z kulshres $
8
9 Contents: Contains the routines :
10           zos_forward (zero-order-scalar forward mode):      define _ZOS_
11           fos_forward (first-order-scalar forward mode):     define _FOS_
12           hos_forward (higher-order-scalar forward mode):    define _HOS_
13           fov_forward (first-order-vector forward mode):     define _FOV_
14           hov_forward (higher-order-vector forward mode):    define _HOV_
15           hov_wk_forward (higher-order-vector forward mode): define _HOV_WK_
16           int_forward_safe:                                  define _INT_FOR_ and _NTIGHT__
17
18           Uses the preprocessor to compile the 7 different object files
19           with/without "keep" parameter:                     define _KEEP_
20
21 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz,
22               Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel,
23               Benjamin Letschert, Kshitij Kulshreshtha
24               
25
26 This file is part of ADOL-C. This software is provided as open source.
27 Any use, reproduction, or distribution of the software constitutes
28 recipient's acceptance of the terms of the accompanying license file.
29
30----------------------------------------------------------------------------*/
31
32#include <adolc/interfaces.h>
33#include <adolc/adalloc.h>
34#include <adolc/taping.h>
35#include "taping_p.h"
36#include "oplate.h"
37#include <adolc/externfcts.h>
38#include "externfcts_p.h"
39
40#include <math.h>
41
42#if defined(ADOLC_DEBUG)
43#include <string.h>
44#endif /* ADOLC_DEBUG */
45
46
47/****************************************************************************/
48/*                                                                   MACROS */
49#undef _ADOLC_VECTOR_
50#undef _HIGHER_ORDER_
51
52
53
54/*--------------------------------------------------------------------------*/
55#if defined(_ZOS_)
56#  define GENERATED_FILENAME "zos_forward"
57
58/*--------------------------------------------------------------------------*/
59#else
60#if defined(_FOS_)
61#define GENERATED_FILENAME "fos_forward"
62
63#define ARGUMENT(indexi,l,i) argument[indexi]
64#define TAYLORS(indexd,l,i)   taylors[indexd]
65
66/*--------------------------------------------------------------------------*/
67#else
68#if defined(_FOV_)
69#define GENERATED_FILENAME "fov_forward"
70
71#define _ADOLC_VECTOR_
72
73#if defined(_CHUNKED_)
74#define ARGUMENT(indexi,l,i) argument[indexi][l+offset]
75#define TAYLORS(indexd,l,i)   taylors[indexd][l+offset]
76#else
77#define ARGUMENT(indexi,l,i) argument[indexi][l]
78#define TAYLORS(indexd,l,i)   taylors[indexd][l]
79#endif
80
81/*--------------------------------------------------------------------------*/
82#else
83#if defined(_HOS_)
84#define GENERATED_FILENAME "hos_forward"
85
86#define _HIGHER_ORDER_
87
88#define ARGUMENT(indexi,l,i) argument[indexi][i]
89#define TAYLORS(indexd,l,i)   taylors[indexd][i]
90
91/*--------------------------------------------------------------------------*/
92#else
93#if defined(_HOV_)
94#define GENERATED_FILENAME "hov_forward"
95
96#define _ADOLC_VECTOR_
97#define _HIGHER_ORDER_
98
99#define ARGUMENT(indexi,l,i) argument[indexi][l][i]
100#define TAYLORS(indexd,l,i)   taylors[indexd][l][i]
101
102/*--------------------------------------------------------------------------*/
103#else
104#if defined(_HOV_WK_)
105#define GENERATED_FILENAME "hov_wk_forward"
106
107#define _ADOLC_VECTOR_
108#define _HIGHER_ORDER_
109
110#define ARGUMENT(indexi,l,i) argument[indexi][l][i]
111#define TAYLORS(indexd,l,i)   taylors[indexd][l][i]
112
113/*--------------------------------------------------------------------------*/
114#else
115#if defined(_INT_FOR_)
116#if defined(_TIGHT_)
117#define GENERATED_FILENAME "int_forward_t"
118#endif
119#if defined(_NTIGHT_)
120#define GENERATED_FILENAME "int_forward_s"
121#endif
122#define ARGUMENT(indexi,l,i) argument[indexi][l]
123#define TAYLORS(indexd,l,i)   taylors[indexd][l]
124/*--------------------------------------------------------------------------*/
125#else
126#if defined(_INDO_)
127void copy_index_domain(int res, int arg, locint **ind_dom);
128void merge_2_index_domains(int res, int arg, locint **ind_dom);
129void combine_2_index_domains(int res, int arg1, int arg2, locint **ind_dom);
130void merge_3_index_domains(int res, int arg1, int arg2, locint **ind_dom);
131
132#define NUMNNZ 20
133#define FMIN_ADOLC(x,y)  ((y<x)?y:x)
134
135#if defined(_INDOPRO_)
136#if defined(_TIGHT_)
137#define GENERATED_FILENAME "indopro_forward_t"
138#endif
139#if defined(_NTIGHT_)
140#define GENERATED_FILENAME "indopro_forward_s"
141#endif
142#endif
143#if defined(_NONLIND_)
144
145typedef struct IndexElement {
146    locint  entry;
147    struct IndexElement* left;
148    struct IndexElement* right;
149}
150IndexElement;
151
152typedef struct IndexElement_sod {
153    locint  entry;
154    struct IndexElement_sod* left;
155}
156IndexElement_sod;
157
158void free_tree(IndexElement* nlf, int num);
159void traverse_crs(IndexElement* fod2, IndexElement_sod* sod, int num);
160void traverse_unary(IndexElement* fod, IndexElement* nonl_dom,  IndexElement* fodi, int num, int maxopind);
161
162#if defined(_TIGHT_)
163#define GENERATED_FILENAME "nonl_ind_forward_t"
164#elif defined(_NTIGHT_)
165#define GENERATED_FILENAME "nonl_ind_forward_s"
166#endif
167#endif
168
169
170/*--------------------------------------------------------------------------*/
171#else
172#error Error ! Define [_ZOS_ | _FOS_ |\
173   _HOS_ | _FOV_ | _HOV_ | _HOV_WK_  | _INT_FOR_SAFE_ | _INT_FOR_TIGHT_ | _INDOPRO_ | _NONLIND_ ] [{_KEEP_}]
174#endif
175#endif
176#endif
177#endif
178#endif
179#endif
180#endif
181#endif
182
183/*--------------------------------------------------------------------------*/
184/*                                                               KEEP stuff */
185#if defined(_KEEP_)
186
187#if defined(_HOV_WK_) /* keep in this vector mode */
188#define IF_KEEP_TAYLOR_CLOSE \
189if (keep){\
190  fprintf(DIAG_OUT,"Succeeding reverse sweep will fail!\n");\
191  taylor_close(0);\
192}
193#define IF_KEEP_WRITE_TAYLOR(res,keep,k,p) \
194    { \
195        UPDATE_TAYLORWRITTEN(keep * k * p) \
196        if (keep) \
197        { \
198            ADOLC_WRITE_SCAYLOR(dp_T0[res]); \
199            if (keep > 1) \
200            write_taylors(res,(keep-1),k,p); \
201        } \
202    }
203#else
204#if defined(_ADOLC_VECTOR_) /* otherwise no keep */
205#define IF_KEEP_TAYLOR_CLOSE
206#define IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
207#else /* _ZOS_, _FOS_, _HOS_ */
208#define IF_KEEP_TAYLOR_CLOSE \
209if (keep){\
210  fprintf(DIAG_OUT,"Otherwise succeeding reverse sweep will fail!\n");\
211  taylor_close(0);\
212}
213#if defined(_ZOS_)
214#define IF_KEEP_WRITE_TAYLOR(res,keep,k,p) \
215    { \
216        UPDATE_TAYLORWRITTEN(keep) \
217        if (keep) \
218            ADOLC_WRITE_SCAYLOR(dp_T0[res]); \
219    }
220#else
221#if defined(_FOS_)
222#define IF_KEEP_WRITE_TAYLOR(res,keep,k,p) \
223    { \
224        UPDATE_TAYLORWRITTEN(keep) \
225        if (keep) \
226        { \
227            ADOLC_WRITE_SCAYLOR(dp_T0[res]); \
228            if (keep > 1) \
229                ADOLC_WRITE_SCAYLOR(dp_T[res]); \
230        } \
231    }
232#else
233#if defined(_HOS_)
234#define IF_KEEP_WRITE_TAYLOR(res,keep,k,p) \
235    { \
236        UPDATE_TAYLORWRITTEN(keep) \
237        if (keep) \
238        { \
239            ADOLC_WRITE_SCAYLOR(dp_T0[res]); \
240            if (keep > 1) \
241                write_taylor(res,keep-1); \
242        } \
243    }
244#endif
245#endif
246#endif
247#endif
248#endif
249
250#else  /* no _KEEP_ */
251#define IF_KEEP_TAYLOR_CLOSE
252#define IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
253#endif
254
255/*--------------------------------------------------------------------------*/
256/*                                                      access to variables */
257#if !defined(_ZOS_)
258#if defined(_FOS_)
259#define TRES         *Tres
260#define TARG         *Targ
261#define TARG1        *Targ1
262#define TARG2        *Targ2
263#define TQO          *Tqo
264
265#define TRES_INC     *Tres
266#define TARG_INC     *Targ
267#define TARG1_INC    *Targ1
268#define TARG2_INC    *Targ2
269#define TQO_INC      *Tqo
270
271#define TRES_DEC     *Tres
272#define TARG_DEC     *Targ
273#define TARG1_DEC    *Targ1
274#define TARG2_DEC    *Targ2
275#define TQO_DEC      *Tqo
276
277#define TRES_FOINC   *Tres
278#define TARG_FOINC   *Targ
279#define TARG1_FOINC  *Targ1
280#define TARG2_FOINC  *Targ2
281#define TQO_FOINC    *Tqo
282
283#define TRES_FODEC   *Tres
284#define DEC_TRES_FO
285#define TARG_FODEC   *Targ
286#define TARG1_FODEC  *Targ1
287#define TARG2_FODEC  *Targ2
288#define TQO_FODEC    *Tqo
289
290#define ASSIGN_T(a,b)  a = &b;
291
292#else
293#if defined(_INT_FOR_)
294#define TRES         *Tres
295#define TARG         *Targ
296#define TARG1        *Targ1
297#define TARG2        *Targ2
298#define TQO          *Tqo
299
300#define TRES_INC     *Tres++
301#define TARG_INC     *Targ++
302#define TARG1_INC    *Targ1++
303#define TARG2_INC    *Targ2++
304#define TQO_INC      *Tqo++
305
306#define TRES_DEC     *Tres--
307#define TARG_DEC     *Targ--
308#define TARG1_DEC    *Targ1--
309#define TARG2_DEC    *Targ2--
310#define TQO_DEC      *Tqo--
311
312#define TRES_FOINC   *Tres++
313#define TARG_FOINC   *Targ++
314#define TARG1_FOINC  *Targ1++
315#define TARG2_FOINC  *Targ2++
316#define TQO_FOINC    *Tqo++
317
318#define TRES_FODEC   *Tres--
319#define TARG_FODEC   *Targ--
320#define TARG1_FODEC  *Targ1--
321#define TARG2_FODEC  *Targ2--
322#define TQO_FODEC    *Tqo--
323
324
325#define ASSIGN_T(a,b)  a = b;
326
327#else  /* _HOS_, _FOV_, _HOV_, _HOV_WK */
328#define TRES         *Tres
329#define TARG         *Targ
330#define TARG1        *Targ1
331#define TARG2        *Targ2
332#define TQO          *Tqo
333
334#define TRES_INC     *Tres++
335#define TARG_INC     *Targ++
336#define TARG1_INC    *Targ1++
337#define TARG2_INC    *Targ2++
338#define TQO_INC      *Tqo++
339
340#define TRES_DEC     *Tres--
341#define TARG_DEC     *Targ--
342#define TARG1_DEC    *Targ1--
343#define TARG2_DEC    *Targ2--
344#define TQO_DEC      *Tqo--
345
346#if defined(_FOV_)
347#define TRES_FOINC   *Tres++
348#define TARG_FOINC   *Targ++
349#define TARG1_FOINC  *Targ1++
350#define TARG2_FOINC  *Targ2++
351#define TQO_FOINC    *Tqo++
352
353#define TRES_FODEC   *Tres
354#define DEC_TRES_FO  Tres--;
355#define TARG_FODEC   *Targ--
356#define TARG1_FODEC  *Targ1--
357#define TARG2_FODEC  *Targ2--
358#define TQO_FODEC    *Tqo--
359#else /* _HOS_, _HOV_, _HOV_WK */
360#define TRES_FOINC   *Tres
361#define TARG_FOINC   *Targ
362#define TARG1_FOINC  *Targ1
363#define TARG2_FOINC  *Targ2
364#define TQO_FOINC    *Tqo
365
366#define TRES_FODEC   *Tres
367#define DEC_TRES_FO
368#define TARG_FODEC   *Targ
369#define TARG1_FODEC  *Targ1
370#define TARG2_FODEC  *Targ2
371#define TQO_FODEC    *Tqo
372#endif
373#endif
374
375#define ASSIGN_T(a,b)  a = b;
376#endif
377#endif
378
379
380/*--------------------------------------------------------------------------*/
381/*                                                               loop stuff */
382#if defined(_ADOLC_VECTOR_)
383#define FOR_0_LE_l_LT_p for (l=0; l<p; l++)
384#define FOR_p_GT_l_GE_0 for (l=p-1; l>=0; l--)
385#else
386#if defined(_INT_FOR_)
387#define FOR_0_LE_l_LT_p for (l=0; l<p; l++)
388#define FOR_p_GT_l_GE_0 for (l=p-1; l>=0; l--)
389#else
390#define FOR_0_LE_l_LT_p
391#define FOR_p_GT_l_GE_0
392#endif
393#endif
394
395#if defined(_HIGHER_ORDER_)
396#define FOR_0_LE_i_LT_k for (i=0; i<k; i++)
397#define FOR_k_GT_i_GE_0 for (i=k-1; i>=0; i--)
398#else
399#define FOR_0_LE_i_LT_k
400#define FOR_k_GT_i_GE_0
401#endif
402
403#if defined(_HOV_)
404#define FOR_0_LE_l_LT_pk for (l=0; l<pk; l++)
405#define INC_pk_1(T)      T += pk-1;
406#define VEC_INC(T,inc)   T += inc;
407#define HOV_INC(T,inc)   T += inc;
408#else
409#if defined(_HOV_WK_)
410#define FOR_0_LE_l_LT_pk for (l=0; l<pk; l++)
411#define INC_pk_1(T)      T += pk-1;
412#define VEC_INC(T,inc)   T += inc;
413#define HOV_INC(T,inc)   T += inc;
414#else
415#if defined(_FOV_)
416#define FOR_0_LE_l_LT_pk for (l=0; l<p; l++)
417#define INC_pk_1(T)      T += p-1;
418#define VEC_INC(T,inc)   T++;
419#define HOV_INC(T,inc)
420#else
421#if defined(_HOS_)
422#define FOR_0_LE_l_LT_pk for (l=0; l<k; l++)
423#define INC_pk_1(T)      T += k-1;
424#define VEC_INC(T,inc)
425#define HOV_INC(T,inc)
426#else
427#if defined(_INT_FOR_)
428#define FOR_0_LE_l_LT_pk for (l=0; l<p; l++)
429#define INC_pk_1(T)      T += p-1;
430#define VEC_INC(T,inc)   T++;
431#else
432#define FOR_0_LE_l_LT_pk
433#define INC_pk_1(T)
434#define VEC_INC(T,inc)
435#define HOV_INC(T,inc)
436#endif
437#endif
438#endif
439#endif
440#endif
441
442/*--------------------------------------------------------------------------*/
443/*                                                        higher order case */
444#if defined(_HIGHER_ORDER_)
445#define BREAK_FOR_I break;
446#else
447#define BREAK_FOR_I ;
448#endif
449
450/* END Macros */
451
452BEGIN_C_DECLS
453
454#if defined(_ZOS_)
455/****************************************************************************/
456/* Zero Order Scalar version of the forward mode.                           */
457/****************************************************************************/
458#if defined(_KEEP_)
459int  zos_forward(
460#else
461int  zos_forward_nk(
462#endif
463    short  tnum,              /* tape id */
464    int    depcheck,          /* consistency chk on # of deps */
465    int    indcheck,          /* consistency chk on # of indeps */
466#if defined(_KEEP_)
467    int    keep,              /* flag for reverse sweep */
468#endif
469    const double *basepoint,  /* independant variable values */
470    double       *valuepoint) /* dependent variable values */
471
472#else
473#if defined(_FOS_)
474/****************************************************************************/
475/* First Order Scalar version of the forward mode.                          */
476/****************************************************************************/
477#if defined(_KEEP_)
478int  fos_forward(
479#else
480int  fos_forward_nk(
481#endif
482    short  tnum,        /* tape id */
483    int    depcheck,    /* consistency chk on # of deps */
484    int    indcheck,    /* consistency chk on # of indeps */
485#if defined(_KEEP_)
486    int    keep,        /* flag for reverse sweep */
487#endif
488    const double *basepoint,  /* independent variable values */
489    double *argument,   /* Taylor coefficients (input) */
490    double *valuepoint, /* Taylor coefficients (output) */
491    double *taylors)    /* matrix of coefficient vectors */
492/* the order of the indices in argument and taylors is [var][taylor] */
493
494#else
495#if defined(_INT_FOR_)
496#if defined(_TIGHT_)
497/****************************************************************************/
498/* First Order Vector version of the forward mode for bit patterns, tight   */
499/****************************************************************************/
500int int_forward_tight(
501    short               tnum,     /* tape id                              */
502    int                 depcheck, /* consistency chk on # of dependents   */
503    int                 indcheck, /* consistency chk on # of independents */
504    int                 p,        /* # of taylor series, bit pattern      */
505    const double       *basepoint,  /* independent variable values   (in)*/
506    unsigned long int **argument,  /* Taylor coeff.                 (in)*/
507    double             *valuepoint, /* dependent variable values    (out)*/
508    unsigned long int **taylors)   /* matrix of coefficient vectors(out)*/
509
510/* int_forward_tight( tag, m, n, p, x[n], X[n][p], y[m], Y[m][p]),
511
512     nBV = number of Boolean Vectors to be packed
513                      (see Chapter Dependence Analysis, ADOL-C Documentation)
514     bits_per_long = 8*sizeof(unsigned long int)
515     p = nBV / bits_per_long + ( (nBV % bits_per_long) != 0 )
516
517     The order of the indices in argument and taylors is [var][taylor]
518
519     For the full Jacobian matrix set
520     p = indep / bits_per_long + ((indep % bits_per_long) != 0)
521     and pass a bit pattern version of the identity matrix as an argument   */
522
523
524#endif
525#if defined (_NTIGHT_)
526/****************************************************************************/
527/* First Order Vector version of the forward mode, bit pattern, safe        */
528/****************************************************************************/
529int int_forward_safe(
530    short             tnum,     /* tape id                              */
531    int               depcheck, /* consistency chk on # of dependents   */
532    int               indcheck, /* consistency chk on # of independents */
533    int               p,        /* # of taylor series, bit pattern      */
534    unsigned long int **argument, /* Taylor coeff.                  (in)*/
535    unsigned long int **taylors)  /* matrix of coefficient vectors (out)*/
536
537/* int_forward_safe( tag, m, n, p, X[n][p], Y[m][p]),
538
539nBV = number of Boolean Vectors to be packed
540(see Chapter Dependence Analysis, ADOL-C Documentation)
541bits_per_long = 8*sizeof(unsigned long int)
542p = nBV / bits_per_long + ( (nBV % bits_per_long) != 0 )
543
544The order of the indices in argument and taylors is [var][taylor]
545
546For the full Jacobian matrix set
547p = indep / bits_per_long + ((indep % bits_per_long) != 0)
548and pass a bit pattern version of the identity matrix as an argument    */
549#endif
550#else
551#if defined(_INDOPRO_)
552#if defined(_TIGHT_)
553/****************************************************************************/
554/* First Order Vector version of the forward mode for bit patterns, tight   */
555/****************************************************************************/
556int indopro_forward_tight(
557    short             tnum,        /* tape id                              */
558    int               depcheck,    /* consistency chk on # of dependents   */
559    int               indcheck,    /* consistency chk on # of independents */
560    const double     *basepoint,  /* independent variable values   (in)   */
561    unsigned int    **crs)        /* returned row index storage (out)     */
562
563/* indopro_forward_tight( tag, m, n, x[n], *crs[m]),
564
565  */
566
567
568#endif
569#if defined (_NTIGHT_)
570/****************************************************************************/
571/* First Order Vector version of the forward mode, bit pattern, safe        */
572/****************************************************************************/
573int indopro_forward_safe(
574    short             tnum,        /* tape id                              */
575    int               depcheck,    /* consistency chk on # of dependents   */
576    int               indcheck,    /* consistency chk on # of independents */
577    const double     *basepoint,   /* independent variable values   (in)   */
578    unsigned int    **crs)         /* returned row index storage (out)     */
579
580/* indopro_forward_safe( tag, m, n, x[n], *crs[m]),
581
582  */
583#endif
584#else
585#if defined(_NONLIND_)
586#if defined(_TIGHT_)
587/****************************************************************************/
588/* First Order Vector version of the forward mode for bit patterns, tight   */
589/****************************************************************************/
590int nonl_ind_forward_tight(
591    short             tnum,        /* tape id                              */
592    int               depcheck,    /* consistency chk on # of dependents   */
593    int               indcheck,    /* consistency chk on # of independents */
594    const double     *basepoint,  /* independent variable values   (in)   */
595    unsigned int     **crs)        /* returned row index storage (out)     */
596
597#endif
598#if defined (_NTIGHT_)
599/****************************************************************************/
600/* First Order Vector version of the forward mode, bit pattern, safe        */
601/****************************************************************************/
602int nonl_ind_forward_safe(
603    short             tnum,        /* tape id                              */
604    int               depcheck,    /* consistency chk on # of dependents   */
605    int               indcheck,    /* consistency chk on # of independents */
606    const double      *basepoint,  /* independent variable values   (in)   */
607    unsigned int    **crs)        /* returned row index storage (out)     */
608
609/* indopro_forward_safe( tag, m, n, x[n], *crs[m]),
610
611  */
612#endif
613#else
614#if defined(_FOV_)
615#if defined(_CHUNKED_)
616/****************************************************************************/
617/* First Order Vector version of the forward mode with p-offset in          */
618/* **argument and **taylors                                                 */
619/****************************************************************************/
620int  fov_offset_forward(
621    short  tnum,        /* tape id */
622    int    depcheck,    /* consistency chk on # of deps */
623    int    indcheck,    /* consistency chk on # of indeps */
624    int    p,           /* # of taylor series */
625    int    offset,      /* offset for assignments */
626    double *basepoint,  /* independent variable values */
627    double **argument,  /* Taylor coefficients (input) */
628    double *valuepoint, /* Taylor coefficients (output) */
629    double **taylors)   /* matrix of coifficient vectors */
630/* the order of the indices in argument and taylors is [var][taylor] */
631#else
632/****************************************************************************/
633/* First Order Vector version of the forward mode.                          */
634/****************************************************************************/
635int  fov_forward(
636    short         tnum,        /* tape id */
637    int           depcheck,    /* consistency chk on # of deps */
638    int           indcheck,    /* consistency chk on # of indeps */
639    int           p,           /* # of taylor series */
640    const double *basepoint,   /* independent variable values */
641    double      **argument,    /* Taylor coefficients (input) */
642    double       *valuepoint,  /* Taylor coefficients (output) */
643    double      **taylors)     /* matrix of coifficient vectors */
644/* the order of the indices in argument and taylors is [var][taylor] */
645#endif
646
647#else
648#if defined(_HOS_)
649/****************************************************************************/
650/* Higher Order Scalar version of the forward mode.                         */
651/****************************************************************************/
652#if defined(_KEEP_)
653int  hos_forward(
654#else
655int  hos_forward_nk(
656#endif
657    short  tnum,        /* tape id */
658    int    depcheck,    /* consistency chk on # of dependents */
659    int    indcheck,    /* consistency chk on # of independents */
660    int    gdegree,     /* highest derivative degree */
661#if defined(_KEEP_)
662    int    keep,        /* flag for reverse sweep */
663#endif
664    double *basepoint,  /* independent variable values */
665    double **argument,  /* independant variable values */
666    double *valuepoint, /* Taylor coefficients (output) */
667    double **taylors)   /* matrix of coifficient vectors */
668
669
670#else
671/****************************************************************************/
672/* Higher Order Vector version of the forward mode.                         */
673/****************************************************************************/
674#if defined(_KEEP_)
675int  hov_wk_forward(
676#else
677int  hov_forward(
678#endif
679    short  tnum,        /* tape id */
680    int    depcheck,    /* consistency chk on # of deps */
681    int    indcheck,    /* consistency chk on # of indeps */
682    int    gdegree,     /* highest derivative degree */
683#if defined(_KEEP_)
684    int    keep,        /* flag for reverse sweep */
685#endif
686    int    p,           /* # of taylor series */
687    double *basepoint,  /* independent variable values */
688    double ***argument, /* Taylor coefficients (input) */
689    double *valuepoint, /* Taylor coefficients (output) */
690    double ***taylors)  /* matrix of coifficient vectors */
691/* the order of the indices in argument and taylors is [var][taylor][deriv] */
692
693#endif
694#endif
695#endif
696#endif
697#endif
698#endif
699#endif
700{
701    /****************************************************************************/
702    /*                                                            ALL VARIABLES */
703
704    unsigned char operation;   /* operation code */
705    int ret_c =3;              /* return value */
706
707    locint size = 0;
708    locint res  = 0;
709    locint arg  = 0;
710    locint arg1 = 0;
711    locint arg2 = 0;
712
713    double coval = 0, *d = 0;
714
715    int indexi = 0,  indexd = 0;
716
717    /* loop indices */
718#if !defined (_ZOS_)
719#if !defined (_INT_FOR_)
720    int i;
721#if !defined (_INDO_)
722    int ii;
723#endif
724#endif
725#endif
726#if defined (_HIGHER_ORDER_)
727    int j, l=0;
728#endif
729    int ls;
730#if defined(_ADOLC_VECTOR_)
731#if !defined (_HIGHER_ORDER_)
732    int l=0;
733#endif
734#endif
735#if defined (_INT_FOR_)
736    int l=0;
737#endif
738#if defined (_INDO_)
739    int l=0;
740    int max_ind_dom;
741#if defined(_INDOPRO_)
742    /* index domains */
743    locint** ind_dom;
744#endif
745#if defined(_NONLIND_)
746    /* nonlinear interaction domains */
747    locint* arg_index;
748    IndexElement* fod;
749    IndexElement* nonl_dom;
750    IndexElement_sod* sod;
751    IndexElement_sod* temp;
752    IndexElement_sod* temp1;
753    int maxopind;
754    int opind;
755    int ii;
756#endif
757#endif
758
759    /* other necessary variables */
760#if !defined (_ZOS_)
761#if !defined (_INDO_)
762#if !defined (_INT_FOR_)
763    double r0=0.0, x, y, divs;
764    int even;
765#endif
766#endif
767#endif
768
769#if defined(_INT_FOR_)
770#ifdef _TIGHT_
771    double  *dp_T0;
772    double y, divs;
773#endif /* _TIGHT_ */
774
775    /* Taylor stuff */
776    unsigned long int  **up_T;
777
778    unsigned long int         *Tres, *Targ, *Targ1, *Targ2;
779#ifdef _TIGHT_
780    unsigned long int         *Tqo;
781    unsigned long int         *Targ1OP, *Targ2OP;
782#endif
783
784#define T0res  T0temp
785#else
786#if defined(_INDO_)
787#ifdef _TIGHT_
788    double  *dp_T0;
789    double  T0temp;
790    double divs;
791#endif /* _TIGHT_ */
792#define T0res  T0temp
793#define T0arg  T0temp
794
795#else
796    double *dp_T0;
797#if !defined(_ZOS_)
798#if  defined(_FOS_)
799    double  *dp_T;
800# define T_TEMP Ttemp;
801# else
802    double *dp_Ttemp, **dpp_T;
803#endif
804    double         *Tres, *Targ, *Targ1, *Targ2, *Tqo;
805
806#if defined (_HIGHER_ORDER_)
807    double         *TresOP, *TresOP2, *zOP;
808    double *dp_z;
809#endif
810   double         *TargOP, *Targ1OP, *Targ2OP;
811   double         T0temp;
812#endif
813#define T0res  T0temp
814#define T0arg  T0temp
815#endif
816#endif
817
818#if defined(_HIGHER_ORDER_)
819    int k = gdegree;
820#endif
821
822#if defined(_KEEP_)
823    int taylbuf=0;
824#endif
825
826#if defined(_HOV_)
827    int pk = k*p;
828#else
829#if defined(_HOV_WK_)
830    int pk = k*p;
831#endif
832#endif
833
834    /* extern diff. function variables */
835#if defined(_EXTERN_)
836#  undef (_EXTERN_)
837#endif
838    /* ZOS_FORWARD */
839#if defined(_ZOS_)
840#   define _EXTERN_ 1
841#   define ADOLC_EXT_FCT_POINTER zos_forward
842#   define ADOLC_EXT_FCT_COMPLETE \
843    zos_forward(n, edfct->dp_x, m, edfct->dp_y)
844#   define ADOLC_EXT_LOOP
845#   define ADOLC_EXT_SUBSCRIPT
846#endif
847    /* FOS_FORWARD */
848#if defined(_FOS_)
849#   define _EXTERN_ 1
850#   define ADOLC_EXT_FCT_POINTER fos_forward
851#   define ADOLC_EXT_FCT_COMPLETE \
852    fos_forward(n, edfct->dp_x, edfct->dp_X, m, edfct->dp_y, edfct->dp_Y)
853#   define ADOLC_EXT_POINTER_X edfct->dp_X
854#   define ADOLC_EXT_POINTER_Y edfct->dp_Y
855#   define ADOLC_EXT_LOOP
856#   define ADOLC_EXT_SUBSCRIPT
857#endif
858    /* FOV_FORWARD */
859#if defined(_FOV_)
860#   define _EXTERN_ 1
861#   define ADOLC_EXT_FCT_POINTER fov_forward
862#   define ADOLC_EXT_FCT_COMPLETE \
863    fov_forward(n, edfct->dp_x, p, edfct->dpp_X, m, edfct->dp_y, edfct->dpp_Y)
864#   define ADOLC_EXT_POINTER_X edfct->dpp_X
865#   define ADOLC_EXT_POINTER_Y edfct->dpp_Y
866#   define ADOLC_EXT_LOOP for (loop2 = 0; loop2 < p; ++loop2)
867#   define ADOLC_EXT_SUBSCRIPT [loop2]
868#endif
869
870#if defined(_EXTERN_)
871    locint n, m;
872    ext_diff_fct *edfct;
873    int loop;
874#   if defined(_FOV_)
875        int loop2;
876#   endif
877    int ext_retc;
878#endif
879
880    ADOLC_OPENMP_THREAD_NUMBER;
881
882#if defined(ADOLC_DEBUG)
883    /****************************************************************************/
884    /*                                                           DEBUG MESSAGES */
885    fprintf(DIAG_OUT,"Call of %s(..) with tag: %d, n: %d, m %d,\n",
886            GENERATED_FILENAME, tnum, indcheck, depcheck);
887#if defined(_KEEP_)
888    fprintf(DIAG_OUT,"                    keep: %d\n", keep);
889#endif
890#if defined(_HIGHER_ORDER_)
891    fprintf(DIAG_OUT,"                    degree: %d\n",gdegree);
892#endif
893#if defined(_ADOLC_VECTOR_)
894    fprintf(DIAG_OUT,"                    p: %d\n\n",p);
895#endif
896
897#endif
898
899    /****************************************************************************/
900    /*                                                                    INITs */
901
902   /* Set up stuff for the tape */
903    ADOLC_OPENMP_GET_THREAD_NUMBER;
904
905    /* Initialize the Forward Sweep */
906
907    init_for_sweep(tnum);
908
909      if ( (depcheck != ADOLC_CURRENT_TAPE_INFOS.stats[NUM_DEPENDENTS]) ||
910            (indcheck != ADOLC_CURRENT_TAPE_INFOS.stats[NUM_INDEPENDENTS]) ) {
911        fprintf(DIAG_OUT,"ADOL-C error: forward sweep on tape %d  aborted!\n"
912                "Number of dependent(%u) and/or independent(%u) variables passed"
913                " to forward is\ninconsistent with number "
914                "recorded on tape (%d, %d) \n", tnum,
915                depcheck, indcheck,
916                ADOLC_CURRENT_TAPE_INFOS.stats[NUM_DEPENDENTS],
917                ADOLC_CURRENT_TAPE_INFOS.stats[NUM_INDEPENDENTS]);
918        exit (-1);
919    }
920
921    /****************************************************************************/
922    /*                                                        MEMORY ALLOCATION */
923    /* olvo 980626 has to be revised for common blocks */
924
925    /*--------------------------------------------------------------------------*/
926#if !defined(_NTIGHT_)
927    dp_T0 = myalloc1(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES]);
928    ADOLC_CURRENT_TAPE_INFOS.dp_T0 = dp_T0;
929#endif /* !_NTIGHT_ */
930#if defined(_ZOS_)                                                   /* ZOS */
931
932#if defined(_KEEP_)
933    if (keep>1) {
934        fprintf(DIAG_OUT,"\n ADOL-C error: zero order scalar forward cannot save"
935                " more\nthan zero order taylor coefficients!\n");
936        exit (-1);
937    }
938#endif
939#if defined(_KEEP_)
940    if (keep) {
941      taylbuf = ADOLC_CURRENT_TAPE_INFOS.stats[TAY_BUFFER_SIZE];
942
943        taylor_begin(taylbuf,&dp_T0,keep-1);
944    }
945#endif
946
947    /*--------------------------------------------------------------------------*/
948#else                                                                /* FOS */
949#if defined(_FOS_)
950#if defined(_KEEP_)
951    if (keep>2) {
952        fprintf(DIAG_OUT,"\n ADOL-C error: first order scalar forward cannot save"
953                " more  \nthan first order taylor coefficients!\n");
954        exit (-1);
955    }
956#endif
957    dp_T = myalloc1(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES]);
958# define TAYLOR_BUFFER dp_T
959#if defined(_KEEP_)
960    if (keep) {
961        taylbuf = ADOLC_CURRENT_TAPE_INFOS.stats[TAY_BUFFER_SIZE];
962        taylor_begin(taylbuf,&dp_T,keep-1);
963    }
964#endif
965
966    /*--------------------------------------------------------------------------*/
967#else                                                                /* INF_FOR */
968#if defined(_INT_FOR_)
969        up_T     = myalloc2_ulong(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES],p);
970#define TAYLOR_BUFFER up_T
971
972    /*--------------------------------------------------------------------------*/
973#else                                                                /* INDOPRO */
974#if defined(_INDO_)
975#if defined(_INDOPRO_)
976    /* index domains */
977    ind_dom = (locint **)  malloc(sizeof(locint*) * ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES]);
978    max_ind_dom = ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES];
979
980    for(i=0;i<max_ind_dom;i++)
981    {
982        ind_dom[i] = (locint *)  malloc(sizeof(locint) * (NUMNNZ+2));
983        ind_dom[i][0] = 0;
984        ind_dom[i][1] = NUMNNZ;
985    }
986#endif
987#if defined(_NONLIND_)
988        maxopind=ADOLC_CURRENT_TAPE_INFOS.stats[NUM_OPERATIONS]+ADOLC_CURRENT_TAPE_INFOS.stats[NUM_EQ_PROD];
989        fod = (struct IndexElement*) malloc(sizeof(struct IndexElement) * maxopind);
990        nonl_dom = (struct IndexElement*) malloc(sizeof(struct IndexElement) * indcheck);
991        arg_index = (locint *)  malloc(sizeof(locint) * (ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES]));
992        opind = 0;
993        sod = (struct IndexElement_sod*) malloc(sizeof(struct IndexElement_sod) * indcheck);
994
995        for(i=0;i<indcheck;i++)
996          {
997            sod[i].entry = 0;
998            sod[i].left = NULL;
999            nonl_dom[i].entry = indcheck+2;
1000            nonl_dom[i].left = NULL;
1001            nonl_dom[i].right = NULL;
1002          }
1003#endif
1004
1005    /*--------------------------------------------------------------------------*/
1006#else                                                                /* FOV */
1007#if defined(_FOV_)
1008    dpp_T = myalloc2(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES],p);
1009# define TAYLOR_BUFFER dpp_T
1010    dp_Ttemp = myalloc1(p);
1011# define T_TEMP dp_Ttemp;
1012
1013    /*--------------------------------------------------------------------------*/
1014#else                                                                /* HOS */
1015#if defined(_HOS_)
1016    dpp_T = myalloc2(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES],k);
1017# define TAYLOR_BUFFER dpp_T
1018    dp_z  = myalloc1(k);
1019    dp_Ttemp = myalloc1(k);
1020# define T_TEMP dp_Ttemp;
1021#if defined(_KEEP_)
1022    if (keep) {
1023        taylbuf = ADOLC_CURRENT_TAPE_INFOS.stats[TAY_BUFFER_SIZE];
1024        taylor_begin(taylbuf,dpp_T,keep-1);
1025    }
1026#endif
1027
1028    /*--------------------------------------------------------------------------*/
1029#else                                                     /* HOV and HOV_WK */
1030    dpp_T = myalloc2(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES],p*k);
1031# define TAYLOR_BUFFER dpp_T
1032    dp_z  = myalloc1(k);
1033    dp_Ttemp = myalloc1(p*k);
1034# define T_TEMP dp_Ttemp;
1035#if defined(_KEEP_)
1036    if (keep) {
1037        taylbuf = ADOLC_CURRENT_TAPE_INFOS.stats[TAY_BUFFER_SIZE];
1038        taylor_begin(taylbuf,dpp_T,keep-1);
1039    }
1040#endif
1041#endif
1042#endif
1043#endif
1044#endif
1045#endif
1046#endif
1047    /****************************************************************************/
1048    /*                                                            FORWARD SWEEP */
1049
1050#if defined(ADOLC_DEBUG)
1051/* #include <string.h> */
1052    int v = 0;
1053    unsigned int countPerOperation[256], taylorPerOperation[256];
1054    memset(countPerOperation, 0, 1024);
1055    memset(taylorPerOperation, 0, 1024);
1056#   define UPDATE_TAYLORWRITTEN(X) taylorPerOperation[operation] += X;
1057#else
1058#   define UPDATE_TAYLORWRITTEN(X)
1059#endif /* ADOLC_DEBUG */
1060
1061    operation=get_op_f();
1062#if defined(ADOLC_DEBUG)
1063    ++countPerOperation[operation];
1064#endif /* ADOLC_DEBUG */
1065
1066    while (operation !=end_of_tape) {
1067
1068
1069      switch (operation) {
1070
1071
1072                /****************************************************************************/
1073                /*                                                                  MARKERS */
1074
1075                /*--------------------------------------------------------------------------*/
1076            case end_of_op:                                          /* end_of_op */
1077                get_op_block_f();
1078                operation=get_op_f();
1079                /* Skip next operation, it's another end_of_op */
1080                break;
1081
1082                /*--------------------------------------------------------------------------*/
1083            case end_of_int:                                        /* end_of_int */
1084                get_loc_block_f();
1085                break;
1086
1087                /*--------------------------------------------------------------------------*/
1088            case end_of_val:                                        /* end_of_val */
1089               get_val_block_f();
1090                break;
1091                /*--------------------------------------------------------------------------*/
1092            case start_of_tape:                                  /* start_of_tape */
1093            case end_of_tape:                                      /* end_of_tape */
1094                break;
1095
1096
1097                /****************************************************************************/
1098                /*                                                               COMPARISON */
1099
1100                /*--------------------------------------------------------------------------*/
1101            case eq_zero:                                              /* eq_zero */
1102                arg = get_locint_f();
1103
1104#if !defined(_NTIGHT_)
1105                if (dp_T0[arg] != 0) {
1106                    if (ADOLC_GLOBAL_TAPE_VARS.branchSwitchWarning)
1107                        fprintf(DIAG_OUT,
1108                                "ADOL-C Warning: Branch switch detected in comparison "
1109                                "(operator eq_zero).\n"
1110                                "Forward sweep aborted! Retaping recommended!\n");
1111                    ret_c = -1;
1112                    operation = end_of_tape;
1113                    continue;
1114                }
1115                ret_c = 0;
1116#endif /* !_NTIGHT_ */
1117                break;
1118
1119                /*--------------------------------------------------------------------------*/
1120            case neq_zero:                                            /* neq_zero */
1121                arg = get_locint_f();
1122
1123#if !defined(_NTIGHT_)
1124                if (dp_T0[arg] == 0) {
1125                    if (ADOLC_GLOBAL_TAPE_VARS.branchSwitchWarning)
1126                        fprintf(DIAG_OUT,
1127                                "ADOL-C Warning: Branch switch detected in comparison "
1128                                "(operator neq_zero).\n"
1129                                "Forward sweep aborted! Retaping recommended!\n");
1130                    ret_c = -1;
1131                    operation = end_of_tape;
1132                    continue;
1133                }
1134#endif /* !_NTIGHT_ */
1135                break;
1136
1137                /*--------------------------------------------------------------------------*/
1138            case le_zero:                                              /* le_zero */
1139                arg = get_locint_f();
1140
1141#if !defined(_NTIGHT_)
1142                if (dp_T0[arg] > 0) {
1143                    if (ADOLC_GLOBAL_TAPE_VARS.branchSwitchWarning)
1144                        fprintf(DIAG_OUT,
1145                                "ADOL-C Warning: Branch switch detected in comparison "
1146                                "(operator le_zero).\n"
1147                                "Forward sweep aborted! Retaping recommended!\n");
1148                    ret_c = -1;
1149                    operation = end_of_tape;
1150                    continue;
1151                }
1152                if (dp_T0[arg] == 0)
1153                    ret_c = 0;
1154#endif /* !_NTIGHT_ */
1155                break;
1156
1157                /*--------------------------------------------------------------------------*/
1158            case gt_zero:                                              /* gt_zero */
1159                arg = get_locint_f();
1160
1161#if !defined(_NTIGHT_)
1162                if (dp_T0[arg] <= 0) {
1163                    if (ADOLC_GLOBAL_TAPE_VARS.branchSwitchWarning)
1164                        fprintf(DIAG_OUT,
1165                                "ADOL-C Warning: Branch switch detected in comparison "
1166                                "(operator gt_zero).\n"
1167                                "Forward sweep aborted! Retaping recommended!\n");
1168                    ret_c = -1;
1169                    operation = end_of_tape;
1170                    continue;
1171                }
1172#endif /* !_NTIGHT_ */
1173                break;
1174
1175                /*--------------------------------------------------------------------------*/
1176            case ge_zero:                                              /* ge_zero */
1177                arg = get_locint_f();
1178
1179#if !defined(_NTIGHT_)
1180                if (dp_T0[arg] < 0) {
1181                    if (ADOLC_GLOBAL_TAPE_VARS.branchSwitchWarning)
1182                        fprintf(DIAG_OUT,
1183                                "ADOL-C Warning: Branch switch detected in comparison "
1184                                "(operator ge_zero).\n"
1185                                "Forward sweep aborted! Retaping recommended!\n");
1186                    ret_c = -1;
1187                    operation = end_of_tape;
1188                    continue;
1189                }
1190                if (dp_T0[arg] == 0)
1191                    ret_c = 0;
1192#endif /* !_NTIGHT_ */
1193                break;
1194
1195                /*--------------------------------------------------------------------------*/
1196            case lt_zero:                                              /* lt_zero */
1197                arg = get_locint_f();
1198
1199#if !defined(_NTIGHT_)
1200                if (dp_T0[arg] >= 0) {
1201                    if (ADOLC_GLOBAL_TAPE_VARS.branchSwitchWarning)
1202                        fprintf(DIAG_OUT,
1203                                "ADOL-C Warning: Branch switch detected in comparison "
1204                                "(operator lt_zero).\n"
1205                                "Forward sweep aborted! Retaping recommended!\n");
1206                    ret_c = -1;
1207                    operation = end_of_tape;
1208                    continue;
1209                }
1210#endif /* !_NTIGHT_ */
1211                break;
1212
1213
1214                /****************************************************************************/
1215                /*                                                              ASSIGNMENTS */
1216
1217                /*--------------------------------------------------------------------------*/
1218            case assign_a:           /* assign an adouble variable an    assign_a */
1219                /* adouble value. (=) */
1220                arg = get_locint_f();
1221                res = get_locint_f();
1222
1223                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1224
1225
1226#if !defined(_NTIGHT_)
1227                dp_T0[res] = dp_T0[arg];
1228#endif /* !_NTIGHT_ */
1229
1230#if defined(_INDO_)
1231#if defined(_INDOPRO_)
1232                copy_index_domain(res, arg, ind_dom);
1233#endif
1234#if defined(_NONLIND_)
1235                arg_index[res] = arg_index[arg];               
1236#endif           
1237#else
1238#if !defined(_ZOS_) /* BREAK_ZOS */
1239                ASSIGN_T(Targ,TAYLOR_BUFFER[arg])
1240                ASSIGN_T(Tres,TAYLOR_BUFFER[res])
1241
1242                FOR_0_LE_l_LT_pk
1243                TRES_INC = TARG_INC;
1244#endif
1245#endif /* ALL_TOGETHER_AGAIN */
1246                break;
1247
1248                /*--------------------------------------------------------------------------*/
1249            case assign_d:            /* assign an adouble variable a    assign_d */
1250                /* double value. (=) */
1251                res   = get_locint_f();
1252                coval = get_val_f();
1253
1254                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1255
1256#if !defined(_NTIGHT_)
1257                dp_T0[res] = coval;
1258#endif /* !_NTIGHT_ */
1259
1260#if defined(_INDO_)
1261#if defined(_INDOPRO_)
1262                ind_dom[res][0]=0;
1263#endif
1264#if defined(_NONLIND_)
1265                fod[opind].entry = maxopind+2;
1266                fod[opind].left = NULL;
1267                fod[opind].right = NULL;
1268                arg_index[res] = opind++;               
1269#endif
1270#else
1271#if !defined(_ZOS_) /* BREAK_ZOS */
1272                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1273
1274                FOR_0_LE_l_LT_pk
1275                TRES_INC = 0;
1276#endif
1277#endif /* ALL_TOGETHER_AGAIN */
1278                break;
1279
1280                /*--------------------------------------------------------------------------*/
1281            case assign_d_zero:  /* assign an adouble variable a    assign_d_zero */
1282                /* double value. (0) (=) */
1283                res   = get_locint_f();
1284
1285                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1286
1287#if !defined(_NTIGHT_)
1288                dp_T0[res] = 0.0;
1289#endif /* !_NTIGHT_ */
1290
1291#if defined(_INDO_)
1292#if defined(_INDOPRO_)
1293                ind_dom[res][0]=0;
1294#endif
1295#if defined(_NONLIND_)
1296                fod[opind].entry = maxopind+2;
1297                fod[opind].left = NULL;
1298                fod[opind].right = NULL;
1299                arg_index[res] = opind++;               
1300#endif
1301#else
1302#if !defined(_ZOS_) /* BREAK_ZOS */
1303                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1304
1305                FOR_0_LE_l_LT_pk
1306                TRES_INC = 0;
1307#endif
1308#endif /* ALL_TOGETHER_AGAIN */
1309                break;
1310
1311                /*--------------------------------------------------------------------------*/
1312            case assign_d_one:    /* assign an adouble variable a    assign_d_one */
1313                /* double value. (1) (=) */
1314                res   = get_locint_f();
1315
1316                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1317
1318#if !defined(_NTIGHT_)
1319                dp_T0[res] = 1.0;
1320#endif /* !_NTIGHT_ */
1321
1322#if defined(_INDO_)
1323#if defined(_INDOPRO_)
1324                ind_dom[res][0]=0;
1325#endif
1326#if defined(_NONLIND_)
1327                fod[opind].entry = maxopind+2;
1328                fod[opind].left = NULL;
1329                fod[opind].right = NULL;
1330                arg_index[res] = opind++;               
1331#endif
1332#else
1333#if !defined(_ZOS_) /* BREAK_ZOS */
1334                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1335
1336                FOR_0_LE_l_LT_pk
1337                TRES_INC = 0;
1338
1339#endif
1340#endif /* ALL_TOGETHER_AGAIN */
1341                break;
1342
1343                /*--------------------------------------------------------------------------*/
1344            case assign_ind:       /* assign an adouble variable an    assign_ind */
1345                /* independent double value (<<=) */
1346                res = get_locint_f();
1347
1348                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1349
1350#if !defined(_NTIGHT_)
1351                dp_T0[res] = basepoint[indexi];
1352#endif /* !_NTIGHT_ */
1353
1354#if defined(_INDO_)
1355#if defined(_INDOPRO_)
1356                ind_dom[res][0] = 1;
1357                ind_dom[res][2] = indexi;
1358#endif         
1359#if defined(_NONLIND_)
1360                fod[opind].entry = indexi;
1361                fod[opind].left = NULL;
1362                fod[opind].right = NULL;
1363                arg_index[res] = opind++;               
1364#endif         
1365#else
1366#if !defined(_ZOS_) /* BREAK_ZOS */
1367                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1368
1369#ifdef _INT_FOR_
1370                FOR_0_LE_l_LT_p
1371                TRES_INC = ARGUMENT(indexi,l,i);
1372#else
1373                FOR_0_LE_l_LT_p
1374                FOR_0_LE_i_LT_k
1375                TRES_INC = ARGUMENT(indexi,l,i);
1376#endif
1377#endif
1378#endif /* ALL_TOGETHER_AGAIN */
1379                ++indexi;
1380                break;
1381
1382                /*--------------------------------------------------------------------------*/
1383            case assign_dep:           /* assign a float variable a    assign_dep */
1384                /* dependent adouble value. (>>=) */
1385                res = get_locint_f();
1386
1387#if !defined(_INDO_)
1388#if !defined(_NTIGHT_)
1389                if ( valuepoint != NULL )
1390                  valuepoint[indexd] = dp_T0[res];
1391#endif /* !_NTIGHT_ */
1392#endif
1393
1394#if defined(_INDO_)
1395#if defined(_INDOPRO_)
1396          if (ind_dom[res][0] != 0) {
1397            crs[indexd] = (unsigned int*) malloc(sizeof(unsigned int) * (ind_dom[res][0]+1));
1398            crs[indexd][0] = ind_dom[res][0];
1399            for(l=1;l<=crs[indexd][0];l++) {
1400              crs[indexd][l] = ind_dom[res][l+1];
1401            }
1402          }
1403          else {
1404            crs[indexd] = (unsigned int*) malloc(sizeof(unsigned int));
1405            crs[indexd][0] =0;
1406          }
1407#endif
1408#else
1409#if !defined(_ZOS_) /* BREAK_ZOS */
1410                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1411
1412#ifdef _INT_FOR_
1413                if (taylors != 0 )  /* ??? question: why here? */
1414                    FOR_0_LE_l_LT_p
1415                    TAYLORS(indexd,l,i) = TRES_INC;
1416#else
1417                if (taylors != 0 )  /* ??? question: why here? */
1418                    FOR_0_LE_l_LT_p
1419                    FOR_0_LE_i_LT_k
1420                    TAYLORS(indexd,l,i) = TRES_INC;
1421#endif
1422#endif
1423#endif /* ALL_TOGETHER_AGAIN */
1424                indexd++;
1425                break;
1426
1427
1428                /****************************************************************************/
1429                /*                                                   OPERATION + ASSIGNMENT */
1430
1431                /*--------------------------------------------------------------------------*/
1432            case eq_plus_d:            /* Add a floating point to an    eq_plus_d */
1433                /* adouble. (+=) */
1434                res   = get_locint_f();
1435                coval = get_val_f();
1436
1437                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1438
1439#if !defined(_NTIGHT_)
1440                dp_T0[res] += coval;
1441#endif /* !_NTIGHT_ */
1442                break;
1443
1444                /*--------------------------------------------------------------------------*/
1445            case eq_plus_a:             /* Add an adouble to another    eq_plus_a */
1446                /* adouble. (+=) */
1447                arg = get_locint_f();
1448                res = get_locint_f();
1449
1450                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1451
1452#if !defined(_NTIGHT_)
1453                dp_T0[res] += dp_T0[arg];
1454#endif /* !_NTIGHT_ */
1455
1456#if defined(_INDO_)
1457#if defined(_INDOPRO_)
1458                merge_2_index_domains(res, arg, ind_dom);
1459#endif
1460#if defined(_NONLIND_)
1461                fod[opind].entry = maxopind+2;
1462                fod[opind].left = &fod[arg_index[res]];
1463                fod[opind].right = &fod[arg_index[arg]];
1464                arg_index[res] = opind++;               
1465#endif
1466#else
1467#if !defined(_ZOS_) /* BREAK_ZOS */
1468                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1469                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
1470
1471#ifdef _INT_FOR_
1472                FOR_0_LE_l_LT_pk
1473                TRES_INC |= TARG_INC;
1474#else
1475                FOR_0_LE_l_LT_pk
1476                TRES_INC += TARG_INC;
1477#endif
1478#endif
1479#endif /* ALL_TOGETHER_AGAIN */
1480                break;
1481
1482                /*--------------------------------------------------------------------------*/
1483            case eq_min_d:       /* Subtract a floating point from an    eq_min_d */
1484                /* adouble. (-=) */
1485                res = get_locint_f();
1486                coval = get_val_f();
1487
1488                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1489
1490#if !defined(_NTIGHT_)
1491                dp_T0[res] -= coval;
1492#endif /* !_NTIGHT_ */
1493                break;
1494
1495                /*--------------------------------------------------------------------------*/
1496            case eq_min_a:        /* Subtract an adouble from another    eq_min_a */
1497                /* adouble. (-=) */
1498                arg = get_locint_f();
1499                res = get_locint_f();
1500
1501                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1502
1503#if !defined(_NTIGHT_)
1504                dp_T0[res] -= dp_T0[arg];
1505#endif /* !_NTIGHT_ */
1506
1507#if defined(_INDO_)
1508#if defined(_INDOPRO_)
1509                merge_2_index_domains(res, arg, ind_dom);
1510#endif
1511#if defined(_NONLIND_)
1512                fod[opind].entry = maxopind+2;
1513                fod[opind].left = &fod[arg_index[res]];
1514                fod[opind].right = &fod[arg_index[arg]];
1515                arg_index[res] = opind++;               
1516#endif
1517#else
1518#if !defined(_ZOS_) /* BREAK_ZOS */
1519                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1520                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
1521
1522#ifdef _INT_FOR_
1523                FOR_0_LE_l_LT_pk
1524                TRES_INC |= TARG_INC;
1525#else
1526                FOR_0_LE_l_LT_pk
1527                TRES_INC -= TARG_INC;
1528#endif
1529#endif
1530#endif /* ALL_TOGETHER_AGAIN */
1531                break;
1532
1533                /*--------------------------------------------------------------------------*/
1534            case eq_mult_d:              /* Multiply an adouble by a    eq_mult_d */
1535                /* flaoting point. (*=) */
1536                res   = get_locint_f();
1537                coval = get_val_f();
1538
1539                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1540
1541#if !defined(_NTIGHT_)
1542                dp_T0[res] *= coval;
1543#endif /* !_NTIGHT_ */
1544
1545#if !defined(_INDO_)
1546#if !defined(_ZOS_) /* BREAK_ZOS */
1547#if !defined( _INT_FOR_)
1548
1549                FOR_0_LE_l_LT_pk
1550                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1551
1552                FOR_0_LE_l_LT_pk
1553                TRES_INC *= coval;
1554#endif
1555#endif
1556#endif /* ALL_TOGETHER_AGAIN */
1557                break;
1558
1559                /*--------------------------------------------------------------------------*/
1560            case eq_mult_a:       /* Multiply one adouble by another    eq_mult_a */
1561                /* (*=) */
1562                arg = get_locint_f();
1563                res = get_locint_f();
1564
1565                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1566
1567#if defined(_INDO_)
1568#if defined(_INDOPRO_)
1569                merge_2_index_domains(res, arg, ind_dom);
1570#endif
1571#if defined(_NONLIND_)
1572                fod[opind].entry = maxopind+2;
1573                fod[opind].left = &fod[arg_index[res]];
1574                fod[opind].right = &fod[arg_index[arg]];
1575                traverse_unary(&fod[arg_index[res]], nonl_dom, &fod[arg_index[arg]], indcheck+1,maxopind+2);
1576                traverse_unary(&fod[arg_index[arg]], nonl_dom, &fod[arg_index[res]], indcheck+1,maxopind+2);
1577                arg_index[res] = opind++;               
1578#endif
1579#else
1580#if !defined(_ZOS_) /* BREAK_ZOS */
1581                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1582                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
1583
1584                INC_pk_1(Tres)
1585                INC_pk_1(Targ)
1586
1587#ifdef _INT_FOR_
1588                FOR_p_GT_l_GE_0
1589                TRES_FODEC |= TARG_DEC;
1590#else
1591                FOR_p_GT_l_GE_0
1592                FOR_k_GT_i_GE_0
1593                { TRES_FODEC = dp_T0[res]*TARG_DEC +
1594                               TRES*dp_T0[arg];
1595                  DEC_TRES_FO
1596#ifdef _HIGHER_ORDER_
1597                  TresOP = Tres-i;
1598                  TargOP = Targ;
1599
1600                  for (j=0;j<i;j++)
1601                  *Tres += (*TresOP++) * (*TargOP--);
1602                  Tres--;
1603#endif /* _HIGHER_ORDER_ */
1604                }
1605#endif
1606#endif
1607#endif /* ALL_TOGETHER_AGAIN */
1608#if !defined(_NTIGHT_)
1609               dp_T0[res] *= dp_T0[arg];
1610#endif /* !_NTIGHT_ */
1611                break;
1612
1613                /*--------------------------------------------------------------------------*/
1614            case incr_a:                        /* Increment an adouble    incr_a */
1615                res   = get_locint_f();
1616
1617                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1618
1619#if !defined(_NTIGHT_)
1620                dp_T0[res]++;
1621#endif /* !_NTIGHT_ */
1622                break;
1623
1624                /*--------------------------------------------------------------------------*/
1625            case decr_a:                        /* Increment an adouble    decr_a */
1626                res   = get_locint_f();
1627
1628                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1629
1630#if !defined(_NTIGHT_)
1631                dp_T0[res]--;
1632#endif /* !_NTIGHT_ */
1633                break;
1634
1635
1636                /****************************************************************************/
1637                /*                                                        BINARY OPERATIONS */
1638
1639                /*--------------------------------------------------------------------------*/
1640            case plus_a_a:                 /* : Add two adoubles. (+)    plus a_a */
1641                arg1 = get_locint_f();
1642                arg2 = get_locint_f();
1643                res  = get_locint_f();
1644
1645                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1646
1647#if !defined(_NTIGHT_)
1648                dp_T0[res] = dp_T0[arg1] + dp_T0[arg2];
1649#endif /* !_NTIGHT_ */
1650
1651#if defined(_INDO_)
1652#if defined(_INDOPRO_)
1653                combine_2_index_domains(res, arg1, arg2, ind_dom);
1654#endif     
1655#if defined(_NONLIND_)
1656                fod[opind].entry = maxopind+2;
1657                fod[opind].left = &fod[arg_index[arg1]];
1658                fod[opind].right = &fod[arg_index[arg2]];
1659                arg_index[res] = opind++;               
1660#endif
1661#else
1662#if !defined(_ZOS_) /* BREAK_ZOS */
1663                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
1664                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
1665                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
1666
1667#ifdef _INT_FOR_
1668                FOR_0_LE_l_LT_pk
1669                TRES_INC = TARG1_INC | TARG2_INC;
1670#else
1671                FOR_0_LE_l_LT_pk
1672                TRES_INC = TARG1_INC + TARG2_INC;
1673#endif
1674#endif
1675#endif /* ALL_TOGETHER_AGAIN */
1676                break;
1677
1678                /*--------------------------------------------------------------------------*/
1679            case plus_d_a:             /* Add an adouble and a double    plus_d_a */
1680                /* (+) */
1681                arg   = get_locint_f();
1682                res   = get_locint_f();
1683                coval = get_val_f();
1684
1685                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1686
1687#if !defined(_NTIGHT_)
1688                dp_T0[res] = dp_T0[arg] + coval;
1689#endif /* !_NTIGHT_ */
1690
1691#if defined(_INDO_)
1692#if defined(_INDOPRO_)
1693                copy_index_domain(res, arg, ind_dom);
1694#endif
1695#if defined(_NONLIND_)
1696                arg_index[res] = arg_index[arg];
1697#endif               
1698#else
1699#if !defined(_ZOS_) /* BREAK_ZOS */
1700                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1701                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
1702
1703                FOR_0_LE_l_LT_pk
1704                TRES_INC = TARG_INC;
1705#endif
1706#endif /* ALL_TOGETHER_AGAIN */
1707                break;
1708
1709                /*--------------------------------------------------------------------------*/
1710            case min_a_a:              /* Subtraction of two adoubles     min_a_a */
1711                /* (-) */
1712                arg1 = get_locint_f();
1713                arg2 = get_locint_f();
1714                res  = get_locint_f();
1715
1716                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1717
1718#if !defined(_NTIGHT_)
1719                dp_T0[res] = dp_T0[arg1] -
1720                                               dp_T0[arg2];
1721#endif /* !_NTIGHT_ */
1722
1723
1724#if defined(_INDO_)   
1725#if defined(_INDOPRO_)
1726                combine_2_index_domains(res, arg1, arg2, ind_dom);
1727#endif
1728#if defined(_NONLIND_)
1729                fod[opind].entry = maxopind+2;
1730                fod[opind].left = &fod[arg_index[arg1]];
1731                fod[opind].right = &fod[arg_index[arg2]];
1732                arg_index[res] = opind++;               
1733#endif
1734#else
1735#if !defined(_ZOS_) /* BREAK_ZOS */
1736                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
1737                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
1738                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
1739
1740#ifdef _INT_FOR_
1741                FOR_0_LE_l_LT_pk
1742                TRES_INC = TARG1_INC | TARG2_INC;
1743#else
1744                 FOR_0_LE_l_LT_pk
1745                TRES_INC = TARG1_INC - TARG2_INC;
1746#endif
1747#endif
1748#endif /* ALL_TOGETHER_AGAIN */
1749                break;
1750
1751                /*--------------------------------------------------------------------------*/
1752            case min_d_a:                /* Subtract an adouble from a    min_d_a */
1753                /* double (-) */
1754                arg =get_locint_f();
1755                res = get_locint_f();
1756                coval = get_val_f();
1757
1758                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1759
1760#if !defined(_NTIGHT_)
1761                dp_T0[res] = coval - dp_T0[arg];
1762#endif /* !_NTIGHT_ */
1763
1764#if defined(_INDO_)
1765#if defined(_INDOPRO_)
1766                copy_index_domain(res, arg, ind_dom);
1767#endif
1768#if defined(_NONLIND_)
1769                arg_index[res] = arg_index[arg];               
1770#endif
1771#else
1772#if !defined(_ZOS_) /* BREAK_ZOS */
1773                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1774                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
1775
1776#ifdef _INT_FOR_
1777                FOR_0_LE_l_LT_pk
1778                TRES_INC = TARG_INC;
1779#else
1780                FOR_0_LE_l_LT_pk
1781                TRES_INC = -TARG_INC;
1782#endif
1783#endif
1784#endif /* ALL_TOGETHER_AGAIN */
1785                break;
1786
1787                /*--------------------------------------------------------------------------*/
1788            case mult_a_a:               /* Multiply two adoubles (*)    mult_a_a */
1789                arg1 = get_locint_f();
1790                arg2 = get_locint_f();
1791                res  = get_locint_f();
1792
1793                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1794
1795#if defined(_INDO_)
1796#if defined(_INDOPRO_)
1797                combine_2_index_domains(res, arg1, arg2, ind_dom);
1798#endif
1799#if defined(_NONLIND_)
1800                fod[opind].entry = maxopind+2;
1801                fod[opind].left = &fod[arg_index[arg1]];
1802                fod[opind].right = &fod[arg_index[arg2]];
1803                traverse_unary(&fod[arg_index[arg1]], nonl_dom, &fod[arg_index[arg2]], indcheck+1,maxopind+2);
1804                traverse_unary(&fod[arg_index[arg2]], nonl_dom, &fod[arg_index[arg1]], indcheck+1,maxopind+2);
1805                arg_index[res] = opind++;               
1806#endif
1807#else
1808#if !defined(_ZOS_) /* BREAK_ZOS */
1809                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
1810                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
1811                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
1812
1813#ifdef _INT_FOR_
1814                FOR_0_LE_l_LT_p
1815                TRES_FOINC = TARG2_INC | TARG1_INC;
1816#else
1817                /* olvo 980915 now in reverse order to allow x = x*x etc. */
1818                INC_pk_1(Tres)
1819                INC_pk_1(Targ1)
1820                INC_pk_1(Targ2)
1821
1822                FOR_p_GT_l_GE_0
1823                FOR_k_GT_i_GE_0
1824                { TRES_FODEC = dp_T0[arg1]*TARG2_DEC +
1825                               TARG1_DEC*dp_T0[arg2];
1826                  DEC_TRES_FO
1827#if defined(_HIGHER_ORDER_)
1828                  Targ1OP = Targ1-i+1;
1829                  Targ2OP = Targ2;
1830
1831                  for (j=0;j<i;j++) {
1832                  *Tres += (*Targ1OP++) * (*Targ2OP--);
1833                  }
1834                  Tres--;
1835#endif /* _HIGHER_ORDER_ */
1836            }
1837#endif
1838#endif
1839#endif /* ALL_TOGETHER_AGAIN */
1840#if !defined(_NTIGHT_)
1841                dp_T0[res] = dp_T0[arg1] *
1842                                               dp_T0[arg2];
1843#endif /* !_NTIGHT_ */
1844                break;
1845
1846                /*--------------------------------------------------------------------------*/
1847                /* olvo 991122: new op_code with recomputation */
1848            case eq_plus_prod:   /* increment a product of           eq_plus_prod */
1849                /* two adoubles (*) */
1850                arg1 = get_locint_f();
1851                arg2 = get_locint_f();
1852                res  = get_locint_f();
1853
1854#if defined(_INDO_)
1855#if defined(_INDOPRO_)
1856                merge_3_index_domains(res, arg1, arg2, ind_dom);
1857#endif
1858#if defined(_NONLIND_)
1859                // operation: v = v+u*w
1860                // first step: z = u*w, index domains
1861                fod[opind].entry = maxopind+2;
1862                fod[opind].left = &fod[arg_index[arg1]];
1863                fod[opind].right = &fod[arg_index[arg2]];
1864                // first step: z = u*w,
1865                traverse_unary(&fod[arg_index[arg1]], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
1866                traverse_unary(&fod[arg_index[arg2]], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
1867                opind++;
1868                // second step: v = v+z, index domains
1869                fod[opind].entry = maxopind+2;
1870                fod[opind].left = &fod[arg_index[res]];
1871                fod[opind].right = &fod[opind-1];
1872                // second step: v = v+z,
1873                arg_index[res] = opind++;               
1874#endif
1875#else
1876#if !defined(_ZOS_) /* BREAK_ZOS */
1877                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
1878                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
1879                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
1880
1881#ifdef _INT_FOR_
1882                FOR_0_LE_l_LT_p
1883                TRES_FOINC |= TARG2_INC | TARG1_INC;
1884#else
1885                /* olvo 980915 now in reverse order to allow x = x*x etc. */
1886                INC_pk_1(Tres)
1887                INC_pk_1(Targ1)
1888                INC_pk_1(Targ2)
1889
1890                FOR_p_GT_l_GE_0
1891                FOR_k_GT_i_GE_0
1892                { TRES_FODEC += dp_T0[arg1]*TARG2_DEC +
1893                                TARG1_DEC*dp_T0[arg2];
1894                  DEC_TRES_FO
1895#if defined(_HIGHER_ORDER_)
1896                  Targ1OP = Targ1-i+1;
1897                  Targ2OP = Targ2;
1898
1899                  for (j=0;j<i;j++)
1900                  *Tres += (*Targ1OP++) * (*Targ2OP--);
1901                  Tres--;
1902#endif /* _HIGHER_ORDER_ */
1903                }
1904#endif
1905#endif
1906#endif /* ALL_TOGETHER_AGAIN */
1907#if !defined(_NTIGHT_)
1908                dp_T0[res] += dp_T0[arg1] *
1909                                                    dp_T0[arg2];
1910#endif /* !_NTIGHT_ */
1911                break;
1912
1913                /*--------------------------------------------------------------------------*/
1914                /* olvo 991122: new op_code with recomputation */
1915            case eq_min_prod:    /* decrement a product of            eq_min_prod */
1916                /* two adoubles (*) */
1917                arg1 = get_locint_f();
1918                arg2 = get_locint_f();
1919                res  = get_locint_f();
1920
1921#if defined(_INDO_)
1922#if defined(_INDOPRO_)
1923                merge_3_index_domains(res, arg1, arg2, ind_dom);
1924#endif
1925#if defined(_NONLIND_)
1926                // operation: v = v-u*w
1927                // first step: z = u*w, index domains
1928                fod[opind].entry = maxopind+2;
1929                fod[opind].left = &fod[arg_index[arg1]];
1930                fod[opind].right = &fod[arg_index[arg2]];
1931                // first step: z = u*w,
1932                traverse_unary(&fod[arg_index[arg1]], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
1933                traverse_unary(&fod[arg_index[arg2]], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
1934                opind++;
1935                // second step: v = v-z, index domains
1936                fod[opind].entry = maxopind+2;
1937                fod[opind].left = &fod[arg_index[res]];
1938                fod[opind].right = &fod[opind-1];
1939                // second step: v = v-z,
1940                arg_index[res] = opind++;       
1941#endif
1942#else
1943#if !defined(_ZOS_) /* BREAK_ZOS */
1944                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
1945                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
1946                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
1947
1948#ifdef _INT_FOR_
1949                FOR_0_LE_l_LT_p
1950                TRES_FOINC |= TARG2_INC | TARG1_INC;
1951#else
1952                /* olvo 980915 now in reverse order to allow x = x*x etc. */
1953                INC_pk_1(Tres)
1954                INC_pk_1(Targ1)
1955                INC_pk_1(Targ2)
1956
1957                FOR_p_GT_l_GE_0
1958                FOR_k_GT_i_GE_0
1959                { TRES_FODEC -= dp_T0[arg1]*TARG2_DEC +
1960                                TARG1_DEC*dp_T0[arg2];
1961                  DEC_TRES_FO
1962#if defined(_HIGHER_ORDER_)
1963                  Targ1OP = Targ1-i+1;
1964                  Targ2OP = Targ2;
1965
1966                  for (j=0;j<i;j++)
1967                  *Tres -= (*Targ1OP++) * (*Targ2OP--);
1968                  Tres--;
1969#endif /* _HIGHER_ORDER_ */
1970                }
1971#endif
1972#endif
1973#endif /* ALL_TOGETHER_AGAIN */
1974
1975#if !defined(_NTIGHT_)
1976                dp_T0[res] -= dp_T0[arg1] *
1977                                                    dp_T0[arg2];
1978#endif /* !_NTIGHT_ */
1979                break;
1980
1981                /*--------------------------------------------------------------------------*/
1982            case mult_d_a:         /* Multiply an adouble by a double    mult_d_a */
1983                /* (*) */
1984                arg   = get_locint_f();
1985                res   = get_locint_f();
1986                coval = get_val_f();
1987
1988                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1989
1990#if !defined(_NTIGHT_)
1991                dp_T0[res] = dp_T0[arg] * coval;
1992#endif /* !_NTIGHT_ */
1993
1994#if defined(_INDO_)
1995#if defined(_INDOPRO_)
1996                copy_index_domain(res, arg, ind_dom);
1997#endif
1998#if defined(_NONLIND_)
1999                arg_index[res] = arg_index[arg];               
2000#endif
2001#else
2002#if !defined(_ZOS_) /* BREAK_ZOS */
2003                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
2004                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
2005
2006#ifdef _INT_FOR_
2007                FOR_0_LE_l_LT_pk
2008                TRES_INC = TARG_INC;
2009#else
2010                FOR_0_LE_l_LT_pk
2011                TRES_INC = TARG_INC * coval;
2012#endif
2013#endif
2014#endif /* ALL_TOGETHER_AGAIN */
2015                break;
2016
2017                /*--------------------------------------------------------------------------*/
2018            case div_a_a:           /* Divide an adouble by an adouble    div_a_a */
2019                /* (/) */
2020                arg1 = get_locint_f();
2021                arg2 = get_locint_f();
2022                res  = get_locint_f();
2023
2024                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2025
2026#if !defined(_NTIGHT_)
2027#if !defined(_ZOS_) /* BREAK_ZOS */
2028                divs = 1.0 / dp_T0[arg2];
2029#endif /* ALL_TOGETHER_AGAIN */
2030
2031                dp_T0[res] = dp_T0[arg1] /
2032                                               dp_T0[arg2];
2033#endif /* !_NTIGHT_ */
2034
2035#if defined(_INDO_)
2036#if defined(_INDOPRO_)
2037                combine_2_index_domains(res, arg1, arg2, ind_dom);
2038#endif
2039#if defined(_NONLIND_)
2040                fod[opind].entry = maxopind+2;
2041                fod[opind].left = &fod[arg_index[arg1]];
2042                fod[opind].right = &fod[arg_index[arg2]];
2043                traverse_unary(&fod[arg_index[arg1]], nonl_dom, &fod[arg_index[arg2]], indcheck+1,maxopind+2);
2044                traverse_unary(&fod[arg_index[arg2]], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2045                arg_index[res] = opind++;               
2046#endif
2047#else
2048#if !defined(_ZOS_) /* BREAK_ZOS */
2049                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2050                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2051                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2052
2053#ifdef _INT_FOR_
2054                FOR_0_LE_l_LT_p
2055                TRES_FOINC = TARG1_INC | TARG2_FOINC;
2056#else
2057                FOR_0_LE_l_LT_p
2058                FOR_0_LE_i_LT_k
2059                { /* olvo 980922 changed order to allow x = y/x */
2060#if defined(_HIGHER_ORDER_)
2061                    zOP      = dp_z+i;
2062                    (*zOP--) = -(*Targ2) * divs;
2063#endif /* _HIGHER_ORDER_ */
2064
2065                    TRES_FOINC = TARG1_INC * divs + dp_T0[res] *
2066                                 (-TARG2_INC * divs);
2067
2068#if defined(_HIGHER_ORDER_)
2069                    TresOP = Tres-i;
2070
2071                    for (j=0;j<i;j++)
2072                    *Tres += (*TresOP++) * (*zOP--);
2073                    Tres++;
2074#endif /* _HIGHER_ORDER_ */
2075                }
2076#endif
2077#endif
2078#endif /* ALL_TOGETHER_AGAIN */
2079                break;
2080
2081            /*--------------------------------------------------------------------------*/
2082        case div_d_a:             /* Division double - adouble (/)    div_d_a */
2083            arg   = get_locint_f();
2084                res   = get_locint_f();
2085                coval = get_val_f();
2086
2087                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2088
2089                /* olvo 980922 necessary for reverse */
2090                if (arg == res) {
2091                    IF_KEEP_WRITE_TAYLOR(arg,keep,k,p)
2092                }
2093
2094#if !defined(_NTIGHT_)
2095#if !defined(_ZOS_) /* BREAK_ZOS */
2096                divs = 1.0 / dp_T0[arg];
2097#endif /* ALL_TOGETHER_AGAIN */
2098
2099                dp_T0[res] = coval / dp_T0[arg];
2100#endif /* !_NTIGHT_ */
2101
2102#if defined(_INDO_)
2103#if defined(_INDOPRO_)
2104                copy_index_domain(res, arg, ind_dom);
2105#endif
2106#if defined(_NONLIND_)
2107                fod[opind].entry = maxopind+2;
2108                fod[opind].left = &fod[arg_index[arg]];
2109                fod[opind].right = NULL;
2110                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2111                arg_index[res] = opind++;               
2112#endif
2113#else
2114#if !defined(_ZOS_) /* BREAK_ZOS */
2115                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
2116                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
2117
2118#ifdef _INT_FOR_
2119                FOR_0_LE_l_LT_p
2120                TRES_FOINC = TARG_FOINC;
2121#else
2122                FOR_0_LE_l_LT_p
2123                FOR_0_LE_i_LT_k
2124                { /* olvo 980922 changed order to allow x = d/x */
2125#if defined(_HIGHER_ORDER_)
2126                    zOP      = dp_z+i;
2127                    (*zOP--) = -(*Targ) * divs;
2128#endif /* _HIGHER_ORDER_ */
2129
2130                    TRES_FOINC = dp_T0[res] * (-TARG_INC * divs);
2131
2132#if defined(_HIGHER_ORDER_)
2133                    TresOP = Tres-i;
2134
2135                    for (j=0;j<i;j++)
2136                    *Tres += (*TresOP++) * (*zOP--);
2137                    Tres++;
2138#endif /* _HIGHER_ORDER_ */
2139                }
2140#endif
2141#endif
2142#endif /* ALL_TOGETHER_AGAIN */
2143                break;
2144
2145
2146            /****************************************************************************/
2147            /*                                                         SIGN  OPERATIONS */
2148
2149            /*--------------------------------------------------------------------------*/
2150        case pos_sign_a:                                        /* pos_sign_a */
2151            arg   = get_locint_f();
2152                res   = get_locint_f();
2153
2154                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2155
2156#if !defined(_NTIGHT_)
2157                dp_T0[res] = dp_T0[arg];
2158#endif /* !_NTIGHT_ */
2159
2160#if defined(_INDO_)
2161#if defined(_INDOPRO_)
2162                copy_index_domain(res, arg, ind_dom);
2163#endif
2164#if defined(_NONLIND_)
2165                arg_index[res] = arg_index[arg];                           
2166#endif 
2167#else
2168#if !defined(_ZOS_) /* BREAK_ZOS */
2169                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
2170                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
2171
2172                FOR_0_LE_l_LT_pk
2173                TRES_INC = TARG_INC;
2174#endif
2175#endif /* ALL_TOGETHER_AGAIN */
2176                break;
2177
2178                /*--------------------------------------------------------------------------*/
2179            case neg_sign_a:                                        /* neg_sign_a */
2180                arg   = get_locint_f();
2181                res   = get_locint_f();
2182
2183                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2184
2185#if !defined(_NTIGHT_)
2186                dp_T0[res] = -dp_T0[arg];
2187#endif /* !_NTIGHT_ */
2188
2189#if defined(_INDO_)
2190#if defined(_INDOPRO_)
2191                copy_index_domain(res, arg, ind_dom);
2192#endif
2193#if defined(_NONLIND_)
2194                arg_index[res] = arg_index[arg];                           
2195#endif
2196#else
2197#if !defined(_ZOS_) /* BREAK_ZOS */
2198                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
2199                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
2200
2201#ifdef _INT_FOR_
2202                FOR_0_LE_l_LT_pk
2203                TRES_INC = TARG_INC;
2204#else
2205                FOR_0_LE_l_LT_pk
2206                TRES_INC = -TARG_INC;
2207#endif
2208#endif
2209#endif /* ALL_TOGETHER_AGAIN */
2210                break;
2211
2212
2213                /****************************************************************************/
2214                /*                                                         UNARY OPERATIONS */
2215
2216                /*--------------------------------------------------------------------------*/
2217            case exp_op:                          /* exponent operation    exp_op */
2218                arg = get_locint_f();
2219                res = get_locint_f();
2220
2221                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2222
2223#if !defined(_NTIGHT_)
2224                dp_T0[res] = exp(dp_T0[arg]);
2225#endif /* !_NTIGHT_ */
2226
2227                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2228
2229#if defined(_INDO_)
2230#if defined(_INDOPRO_)
2231                copy_index_domain(res, arg, ind_dom);
2232#endif
2233#if defined(_NONLIND_)
2234                fod[opind].entry = maxopind+2;
2235                fod[opind].left = &fod[arg_index[arg]];
2236                fod[opind].right = NULL;
2237                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2238                arg_index[res] = opind++;               
2239#endif
2240#else
2241#if !defined(_ZOS_) /* BREAK_ZOS */
2242                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
2243                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
2244
2245#ifdef _INT_FOR_
2246                FOR_0_LE_l_LT_p
2247                TRES_FOINC = TARG_FOINC;
2248#else
2249                FOR_0_LE_l_LT_p
2250                FOR_0_LE_i_LT_k
2251                { /* olvo 980915 changed order to allow x = exp(x) */
2252#if defined(_HIGHER_ORDER_)
2253                    zOP      = dp_z+i;
2254                    (*zOP--) = (i+1) * (*Targ);
2255#endif /* _HIGHER_ORDER_ */
2256
2257                    TRES_FOINC = dp_T0[res] * TARG_INC;
2258
2259#if defined(_HIGHER_ORDER_)
2260                    TresOP = Tres-i;
2261
2262                    *Tres *= (i+1);
2263                    for (j=0;j<i;j++)
2264                    *Tres += (*TresOP++) * (*zOP--);
2265                    *Tres++ /= (i+1); /* important only for i>0 */
2266#endif /* _HIGHER_ORDER_ */
2267                }
2268
2269#endif
2270#endif
2271#endif /* ALL_TOGETHER_AGAIN */
2272                break;
2273
2274            /*--------------------------------------------------------------------------*/
2275        case sin_op:                              /* sine operation    sin_op */
2276                arg1 = get_locint_f();
2277                arg2 = get_locint_f();
2278                res  = get_locint_f();
2279
2280                IF_KEEP_WRITE_TAYLOR(arg2,keep,k,p) /* olvo 980710 covalue */
2281                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2282
2283#if !defined(_NTIGHT_)
2284                /* Note: always arg2 != arg1 */
2285                dp_T0[arg2] = cos(dp_T0[arg1]);
2286                dp_T0[res]  = sin(dp_T0[arg1]);
2287#endif /* !_NTIGHT_ */
2288
2289                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2290
2291#if defined(_INDO_)
2292#if defined(_INDOPRO_)
2293                copy_index_domain(res, arg1, ind_dom);
2294#endif
2295#if defined(_NONLIND_)
2296                fod[opind].entry = maxopind+2;
2297                fod[opind].left = &fod[arg_index[arg1]];
2298                fod[opind].right = NULL;
2299                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2300                arg_index[res] = opind++;               
2301#endif
2302#else
2303#if !defined(_ZOS_) /* BREAK_ZOS */
2304                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2305                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2306                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2307
2308#ifdef _INT_FOR_
2309                FOR_0_LE_l_LT_p
2310                { /* olvo 980923 changed order to allow x = sin(x) */
2311                    TARG2_FOINC =  TARG1;
2312                    TRES_FOINC  =  TARG1_FOINC;
2313            }
2314#else
2315                FOR_0_LE_l_LT_p
2316                FOR_0_LE_i_LT_k
2317                { /* olvo 980921 changed order to allow x = sin(x) */
2318#if defined(_HIGHER_ORDER_)
2319                    zOP      = dp_z+i;
2320                    (*zOP--) = (i+1) * (*Targ1);
2321#endif /* _HIGHER_ORDER_ */
2322
2323                    /* Note: always arg2 != arg1 */
2324                    TARG2_FOINC = -dp_T0[res]  * TARG1;
2325                    TRES_FOINC  =  dp_T0[arg2] * TARG1_INC;
2326
2327#if defined(_HIGHER_ORDER_)
2328                    TresOP  = Tres-i;
2329                    Targ2OP = Targ2-i;
2330
2331                    *Tres  *= (i+1);
2332                    *Targ2 *= (i+1);
2333                    for (j=0;j<i;j++) {
2334                    *Tres  += (*Targ2OP++) * (*zOP);
2335                        *Targ2 -= (*TresOP++)  * (*zOP--);
2336                    }
2337                    *Targ2++ /= (i+1);
2338                    *Tres++  /= (i+1);
2339#endif /* _HIGHER_ORDER_ */
2340            }
2341#endif
2342#endif
2343#endif /* ALL_TOGETHER_AGAIN */
2344                break;
2345
2346                /*--------------------------------------------------------------------------*/
2347            case cos_op:                            /* cosine operation    cos_op */
2348                arg1 = get_locint_f();
2349                arg2 = get_locint_f();
2350                res  = get_locint_f();
2351
2352                IF_KEEP_WRITE_TAYLOR(arg2,keep,k,p) /* olvo 980710 covalue */
2353                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2354
2355#if !defined(_NTIGHT_)
2356                /* Note: always arg2 != arg1 */
2357                dp_T0[arg2] = sin(dp_T0[arg1]);
2358                dp_T0[res]  = cos(dp_T0[arg1]);
2359#endif /* !_NTIGHT_ */
2360
2361                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2362
2363#if defined(_INDO_)
2364#if defined(_INDOPRO_)
2365                copy_index_domain(res, arg1, ind_dom);
2366#endif
2367#if defined(_NONLIND_)
2368                fod[opind].entry = maxopind+2;
2369                fod[opind].left = &fod[arg_index[arg1]];
2370                fod[opind].right = NULL;
2371                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2372                arg_index[res] = opind++;               
2373#endif
2374
2375#else
2376#if !defined(_ZOS_) /* BREAK_ZOS */
2377                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2378                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2379                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2380
2381#ifdef _INT_FOR_
2382                FOR_0_LE_l_LT_p
2383                { /* olvo 980923 changed order to allow x = cos(x) */
2384                    TARG2_FOINC = TARG1;
2385                    TRES_FOINC  = TARG1_FOINC;
2386            }
2387#else
2388                FOR_0_LE_l_LT_p
2389                FOR_0_LE_i_LT_k
2390                { /* olvo 980921 changed order to allow x = cos(x) */
2391#if defined(_HIGHER_ORDER_)
2392                    zOP      = dp_z+i;
2393                    (*zOP--) = (i+1) * (*Targ1);
2394#endif /* _HIGHER_ORDER_ */
2395
2396                    /* Note: always arg2 != arg1 */
2397                    TARG2_FOINC =  dp_T0[res]  * TARG1;
2398                    TRES_FOINC  = -dp_T0[arg2] * TARG1_INC;
2399
2400#if defined(_HIGHER_ORDER_)
2401                    TresOP  = Tres-i;
2402                    Targ2OP = Targ2-i;
2403
2404                    *Tres  *= (i+1);
2405                    *Targ2 *= (i+1);
2406                    for (j=0;j<i;j++) {
2407                    *Tres  -= (*Targ2OP++) * (*zOP);
2408                        *Targ2 += (*TresOP++)  * (*zOP--);
2409                    }
2410                    *Targ2++ /= (i+1);
2411                    *Tres++  /= (i+1);
2412#endif /* _HIGHER_ORDER_ */
2413            }
2414#endif
2415#endif
2416#endif /* ALL_TOGETHER_AGAIN */
2417                break;
2418
2419                /*--------------------------------------------------------------------------*/
2420            case atan_op:                                              /* atan_op */
2421                arg1 = get_locint_f();
2422                arg2 = get_locint_f();
2423                res  = get_locint_f();
2424
2425                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2426
2427#if !defined(_NTIGHT_)
2428                dp_T0[res]=atan(dp_T0[arg1]);
2429#endif /* !_NTIGHT_ */
2430
2431                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2432
2433#if defined(_INDO_)
2434#if defined(_INDOPRO_)
2435                copy_index_domain(res, arg1, ind_dom);
2436#endif
2437#if defined(_NONLIND_)
2438                fod[opind].entry = maxopind+2;
2439                fod[opind].left = &fod[arg_index[arg1]];
2440                fod[opind].right = NULL;
2441                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2442                arg_index[res] = opind++;               
2443#endif
2444#else
2445#if !defined(_ZOS_) /* BREAK_ZOS */
2446                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2447                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2448#ifdef _INT_FOR_
2449                FOR_0_LE_l_LT_p
2450                TRES_FOINC = TARG1_FOINC;
2451#else
2452                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2453
2454                FOR_0_LE_l_LT_p
2455                { FOR_0_LE_i_LT_k
2456                  { /* olvo 980921 changed order to allow x = atan(x) */
2457#if defined(_HIGHER_ORDER_)
2458                      zOP      = dp_z+i;
2459                      (*zOP--) = (i+1) * (*Targ1);
2460#endif /* _HIGHER_ORDER_ */
2461
2462                      TRES_FOINC = dp_T0[arg2] * TARG1_INC;
2463
2464#if defined(_HIGHER_ORDER_)
2465                      Targ2OP = Targ2;
2466
2467                      *Tres *= (i+1);
2468                      for (j=0;j<i;j++)
2469                      *Tres  += (*Targ2OP++) * (*zOP--);
2470                      *Tres++ /= (i+1);
2471#endif /* _HIGHER_ORDER_ */
2472                  }
2473                  HOV_INC(Targ2, k)
2474                }
2475#endif
2476#endif
2477#endif /* ALL_TOGETHER_AGAIN */
2478                break;
2479
2480            /*--------------------------------------------------------------------------*/
2481        case asin_op:                                              /* asin_op */
2482            arg1 = get_locint_f();
2483                arg2 = get_locint_f();
2484                res  = get_locint_f();
2485
2486                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2487
2488#if !defined(_NTIGHT_)
2489                dp_T0[res] = asin(dp_T0[arg1]);
2490#endif /* !_NTIGHT_ */
2491
2492                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2493
2494#if defined(_INDO_)
2495#if defined(_INDOPRO_)
2496                copy_index_domain(res, arg1, ind_dom);
2497#endif
2498#if defined(_NONLIND_)
2499                fod[opind].entry = maxopind+2;
2500                fod[opind].left = &fod[arg_index[arg1]];
2501                fod[opind].right = NULL;
2502                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2503                arg_index[res] = opind++;               
2504#endif
2505#else
2506#if !defined(_ZOS_) /* BREAK_ZOS */
2507                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2508                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2509#ifdef _INT_FOR_
2510                FOR_0_LE_l_LT_p
2511                TRES_FOINC = TARG1_FOINC;
2512#else
2513                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2514
2515                if (dp_T0[arg1] == 1.0)
2516                    FOR_0_LE_l_LT_p
2517                    { FOR_0_LE_i_LT_k
2518                      if (TARG1 > 0.0) {
2519                      r0 = make_nan();
2520                          VEC_INC(Targ1, k-i)
2521                          BREAK_FOR_I
2522                      } else
2523                          if (TARG1 < 0.0) {
2524                          r0 = make_inf();
2525                              VEC_INC(Targ1, k-i)
2526                              BREAK_FOR_I
2527                          } else {
2528                              r0 = 0.0;
2529                              Targ1++;
2530                          }
2531                  TRES = r0;
2532                  VEC_INC(Tres, k)
2533            } else
2534                    if (dp_T0[arg1] == -1.0)
2535                        FOR_0_LE_l_LT_p
2536                        { FOR_0_LE_i_LT_k
2537                          if (TARG1 > 0.0) {
2538                          r0 = make_inf();
2539                              VEC_INC(Targ1, k-i)
2540                              BREAK_FOR_I
2541                          } else
2542                              if (TARG1 < 0.0) {
2543                              r0 = make_nan();
2544                                  VEC_INC(Targ1, k-i)
2545                                  BREAK_FOR_I
2546                              } else {
2547                                  r0 = 0.0;
2548                                  Targ1++;
2549                              }
2550                  TRES = r0;
2551                  VEC_INC(Tres, k)
2552                } else
2553                        FOR_0_LE_l_LT_p {
2554                            FOR_0_LE_i_LT_k
2555                            { /* olvo 980921 changed order to allow x = asin(x) */
2556#if defined(_HIGHER_ORDER_)
2557                                zOP      = dp_z+i;
2558                                (*zOP--) = (i+1) * (*Targ1);
2559#endif /* _HIGHER_ORDER_ */
2560
2561                                TRES_FOINC = dp_T0[arg2] * TARG1_INC;
2562
2563#if defined(_HIGHER_ORDER_)
2564                                Targ2OP = Targ2;
2565
2566                                *Tres *= (i+1);
2567                                for (j=0;j<i;j++)
2568                                *Tres += (*Targ2OP++) * (*zOP--);
2569                                *Tres++ /= (i+1);
2570#endif /* _HIGHER_ORDER_ */
2571                            }
2572                            HOV_INC(Targ2, k)
2573                        }
2574#endif
2575#endif
2576#endif /* ALL_TOGETHER_AGAIN */
2577                        break;
2578
2579            /*--------------------------------------------------------------------------*/
2580        case acos_op:                                              /* acos_op */
2581                arg1 = get_locint_f();
2582                arg2 = get_locint_f();
2583                res  = get_locint_f();
2584
2585                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2586
2587#if !defined(_NTIGHT_)
2588                dp_T0[res] = acos(dp_T0[arg1]);
2589#endif /* !_NTIGHT_ */
2590
2591                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2592
2593#if defined(_INDO_)
2594#if defined(_INDOPRO_)
2595                copy_index_domain(res, arg1, ind_dom);
2596#endif
2597#if defined(_NONLIND_)
2598                fod[opind].entry = maxopind+2;
2599                fod[opind].left = &fod[arg_index[arg1]];
2600                fod[opind].right = NULL;
2601                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2602                arg_index[res] = opind++;               
2603#endif
2604#else
2605#if !defined(_ZOS_) /* BREAK_ZOS */
2606                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2607                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2608#ifdef _INT_FOR_
2609                FOR_0_LE_l_LT_p
2610                TRES_FOINC = TARG1_FOINC;
2611#else
2612                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2613
2614                if (dp_T0[arg1] == 1.0)
2615                    FOR_0_LE_l_LT_p
2616                    { FOR_0_LE_i_LT_k
2617                      if (TARG1 > 0.0) {
2618                      r0 = make_nan();
2619                          VEC_INC(Targ1, k-i)
2620                          BREAK_FOR_I
2621                      } else
2622                          if (TARG1 < 0.0) {
2623                          r0 = -make_inf();
2624                              VEC_INC(Targ1, k-i)
2625                              BREAK_FOR_I
2626                          } else {
2627                              r0 = 0.0;
2628                              Targ1++;
2629                          }
2630                  TRES = r0;
2631                  VEC_INC(Tres, k)
2632            } else
2633                    if (dp_T0[arg1] == -1.0)
2634                        FOR_0_LE_l_LT_p
2635                        { FOR_0_LE_i_LT_k
2636                          if (TARG1 > 0.0) {
2637                          r0 = -make_inf();
2638                              VEC_INC(Targ1, k-i)
2639                              BREAK_FOR_I
2640                          } else
2641                              if (TARG1 < 0.0) {
2642                              r0 = make_nan();
2643                                  VEC_INC(Targ1, k-i)
2644                                  BREAK_FOR_I
2645                              } else {
2646                                  r0 = 0.0;
2647                                  Targ1++;
2648                              }
2649                  TRES = r0;
2650                  VEC_INC(Tres, k)
2651                } else
2652                        FOR_0_LE_l_LT_p {
2653                            FOR_0_LE_i_LT_k
2654                            { /* olvo 980921 changed order to allow x = acos(x) */
2655#if defined(_HIGHER_ORDER_)
2656                                zOP      = dp_z+i;
2657                                (*zOP--) = (i+1) * (*Targ1);
2658#endif /* _HIGHER_ORDER_ */
2659
2660                                TRES_FOINC = dp_T0[arg2] * TARG1_INC;
2661
2662#if defined(_HIGHER_ORDER_)
2663                                Targ2OP = Targ2;
2664
2665                                *Tres *= (i+1);
2666                                for (j=0;j<i;j++)
2667                                *Tres += (*Targ2OP++) * (*zOP--);
2668                                *Tres++ /= (i+1);
2669#endif /* _HIGHER_ORDER_ */
2670                            }
2671                            HOV_INC(Targ2, k)
2672                        }
2673#endif
2674#endif
2675#endif /* ALL_TOGETHER_AGAIN */
2676                        break;
2677
2678#ifdef ATRIG_ERF
2679
2680            /*--------------------------------------------------------------------------*/
2681        case asinh_op:                                            /* asinh_op */
2682                arg1 = get_locint_f();
2683                arg2 = get_locint_f();
2684                res  = get_locint_f();
2685
2686                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2687
2688#if !defined(_NTIGHT_)
2689                dp_T0[res] = asinh(dp_T0[arg1]);
2690#endif /* !_NTIGHT_ */
2691
2692                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2693
2694#if defined(_INDO_)
2695#if defined(_INDOPRO_)
2696                copy_index_domain(res, arg1, ind_dom);
2697#endif
2698#if defined(_NONLIND_)
2699                fod[opind].entry = maxopind+2;
2700                fod[opind].left = &fod[arg_index[arg1]];
2701                fod[opind].right = NULL;
2702                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2703                arg_index[res] = opind++;               
2704#endif
2705#else
2706#if !defined(_ZOS_) /* BREAK_ZOS */
2707                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2708                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2709#ifdef _INT_FOR_
2710                FOR_0_LE_l_LT_p
2711                TRES_FOINC = TARG1_FOINC;
2712#else
2713                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2714
2715                FOR_0_LE_l_LT_p
2716                { FOR_0_LE_i_LT_k
2717                  { /* olvo 980921 changed order to allow x = asinh(x) */
2718#if defined(_HIGHER_ORDER_)
2719                      zOP      = dp_z+i;
2720                      (*zOP--) = (i+1) * (*Targ1);
2721#endif /* _HIGHER_ORDER_ */
2722
2723                      TRES_FOINC = dp_T0[arg2] * TARG1_INC;
2724
2725#if defined(_HIGHER_ORDER_)
2726                      Targ2OP = Targ2;
2727
2728                      *Tres *= (i+1);
2729                      for (j=0;j<i;j++)
2730                      *Tres += (*Targ2OP++) * (*zOP--);
2731                      *Tres++ /= (i+1);
2732#endif /* _HIGHER_ORDER_ */
2733                  }
2734                  HOV_INC(Targ2, k)
2735                }
2736#endif
2737#endif
2738#endif /* ALL_TOGETHER_AGAIN */
2739                break;
2740
2741            /*--------------------------------------------------------------------------*/
2742        case acosh_op:                                           /* acosh_op */
2743                arg1 = get_locint_f();
2744                arg2 = get_locint_f();
2745                res  = get_locint_f();
2746
2747                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2748
2749#if !defined(_NTIGHT_)
2750                dp_T0[res] = acosh(dp_T0[arg1]);
2751#endif /* !_NTIGHT_ */
2752
2753                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2754
2755#if defined(_INDO_)
2756#if defined(_INDOPRO_)
2757                copy_index_domain(res, arg1, ind_dom);
2758#endif
2759#if defined(_NONLIND_)
2760                fod[opind].entry = maxopind+2;
2761                fod[opind].left = &fod[arg_index[arg1]];
2762                fod[opind].right = NULL;
2763                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2764                arg_index[res] = opind++;               
2765#endif
2766#else
2767#if !defined(_ZOS_) /* BREAK_ZOS */
2768                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2769                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2770#ifdef _INT_FOR_
2771                FOR_0_LE_l_LT_p
2772                TRES_FOINC = TARG1_FOINC;
2773#else
2774                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2775
2776                if (dp_T0[arg1] == 1.0)
2777                    FOR_0_LE_l_LT_p
2778                    { FOR_0_LE_i_LT_k
2779                      if (TARG1 > 0.0) {
2780                      r0 = make_inf();
2781                          VEC_INC(Targ1, k-i)
2782                          BREAK_FOR_I
2783                      } else
2784                          if (TARG1 < 0.0) {
2785                          r0 = make_nan();
2786                              VEC_INC(Targ1, k-i)
2787                              BREAK_FOR_I
2788                          } else {
2789                              r0 = 0.0;
2790                              Targ1++;
2791                          }
2792                  TRES_INC = r0;
2793#if defined(_HIGHER_ORDER_)
2794                  for (i=1;i<k;i++)
2795                  *Tres++ = make_nan();
2796#endif /* _HIGHER_ORDER_ */
2797                } else
2798                    FOR_0_LE_l_LT_p {
2799                        FOR_0_LE_i_LT_k
2800                        { /* olvo 980921 changed order to allow x = acosh(x) */
2801#if defined(_HIGHER_ORDER_)
2802                            zOP      = dp_z+i;
2803                            (*zOP--) = (i+1) * (*Targ1);
2804#endif /* _HIGHER_ORDER_ */
2805
2806                            TRES_FOINC = dp_T0[arg2] * TARG1_INC;
2807
2808#if defined(_HIGHER_ORDER_)
2809                            Targ2OP = Targ2;
2810
2811                            *Tres *= (i+1);
2812                            for (j=0;j<i;j++)
2813                                *Tres += (*Targ2OP++) * (*zOP--);
2814                                *Tres++ /= (i+1);
2815#endif /* _HIGHER_ORDER_ */
2816                            }
2817                            HOV_INC(Targ2, k)
2818                        }
2819#endif
2820#endif
2821#endif /* ALL_TOGETHER_AGAIN */
2822                        break;
2823
2824            /*--------------------------------------------------------------------------*/
2825        case atanh_op:                                            /* atanh_op */
2826                arg1 = get_locint_f();
2827                arg2 = get_locint_f();
2828                res  = get_locint_f();
2829
2830                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2831
2832#if !defined(_NTIGHT_)
2833                dp_T0[res] = atanh(dp_T0[arg1]);
2834#endif /* !_NTIGHT_ */
2835
2836                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2837
2838#if defined(_INDO_)
2839#if defined(_INDOPRO_)
2840                copy_index_domain(res, arg1, ind_dom);
2841#endif
2842#if defined(_NONLIND_)
2843                fod[opind].entry = maxopind+2;
2844                fod[opind].left = &fod[arg_index[arg1]];
2845                fod[opind].right = NULL;
2846                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2847                arg_index[res] = opind++;               
2848#endif
2849#else
2850#if !defined(_ZOS_) /* BREAK_ZOS */
2851                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2852                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2853#ifdef _INT_FOR_
2854                FOR_0_LE_l_LT_p
2855                TRES_FOINC = TARG1_FOINC;
2856#else
2857                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2858
2859                if (dp_T0[arg1] == 1.0)
2860                    FOR_0_LE_l_LT_p
2861                    { FOR_0_LE_i_LT_k
2862                      if (TARG1 > 0.0) {
2863                      r0 = make_nan();
2864                          VEC_INC(Targ1, k-i)
2865                          BREAK_FOR_I
2866                      } else
2867                          if (TARG1 < 0.0) {
2868                          r0 = make_inf();
2869                              VEC_INC(Targ1, k-i)
2870                              BREAK_FOR_I
2871                          } else {
2872                              r0 = 0.0;
2873                              Targ1++;
2874                          }
2875                  TRES_INC = r0;
2876#if defined(_HIGHER_ORDER_)
2877                  for (i=1;i<k;i++)
2878                  *Tres++ = make_nan();
2879#endif /* _HIGHER_ORDER_ */
2880                } else
2881                    if (dp_T0[arg1] == -1.0)
2882                            FOR_0_LE_l_LT_p
2883                            { FOR_0_LE_i_LT_k
2884                              if (TARG1 > 0.0) {
2885                              r0 = make_inf();
2886                                  VEC_INC(Targ1, k-i)
2887                                  BREAK_FOR_I
2888                              } else
2889                                  if (TARG1 < 0.0) {
2890                                  r0 = make_nan();
2891                                      VEC_INC(Targ1, k-i)
2892                                      BREAK_FOR_I
2893                                  } else {
2894                                      r0 = 0.0;
2895                                      Targ1++;
2896                                  }
2897                  TRES_INC = r0;
2898#if defined(_HIGHER_ORDER_)
2899                  for (i=1;i<k;i++)
2900                  *Tres++ = make_nan();
2901#endif /* _HIGHER_ORDER_ */
2902                        } else
2903                            FOR_0_LE_l_LT_p {
2904                                FOR_0_LE_i_LT_k
2905                                { /* olvo 980921 changed order to allow x = atanh(x) */
2906#if defined(_HIGHER_ORDER_)
2907                                    zOP      = dp_z+i;
2908                                    (*zOP--) = (i+1) * (*Targ1);
2909#endif /* _HIGHER_ORDER_ */
2910
2911                                    TRES_FOINC = dp_T0[arg2] * TARG1_INC;
2912
2913#if defined(_HIGHER_ORDER_)
2914                                    Targ2OP = Targ2;
2915
2916                                    *Tres *= (i+1);
2917                                    for (j=0;j<i;j++)
2918                                        *Tres += (*Targ2OP++) * (*zOP--);
2919                                        *Tres++ /= (i+1);
2920#endif /* _HIGHER_ORDER_ */
2921                                    }
2922                                    HOV_INC(Targ2, k)
2923                                }
2924#endif
2925#endif
2926#endif /* ALL_TOGETHER_AGAIN */
2927                                break;
2928
2929            /*--------------------------------------------------------------------------*/
2930        case erf_op:                                                /* erf_op */
2931                arg1 = get_locint_f();
2932                arg2 = get_locint_f();
2933                res  = get_locint_f();
2934
2935                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2936
2937#if !defined(_NTIGHT_)
2938                dp_T0[res] = erf(dp_T0[arg1]);
2939#endif /* !_NTIGHT_ */
2940
2941                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2942#if defined(_INDO_)
2943#if defined(_INDOPRO_)
2944                copy_index_domain(res, arg1, ind_dom);
2945#endif
2946#if defined(_NONLIND_)
2947                fod[opind].entry = maxopind+2;
2948                fod[opind].left = &fod[arg_index[arg1]];
2949                fod[opind].right = NULL;
2950                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2951                arg_index[res] = opind++;               
2952#endif       
2953#else
2954#if !defined(_ZOS_) /* BREAK_ZOS */
2955                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
2956                ASSIGN_T(Targ1,TAYLOR_BUFFER[arg1])
2957#ifdef _INT_FOR_
2958                FOR_0_LE_l_LT_p
2959                TRES_FOINC = TARG1_FOINC;
2960#else
2961                ASSIGN_T(Targ2,TAYLOR_BUFFER[arg2])
2962
2963                FOR_0_LE_l_LT_p
2964                { FOR_0_LE_i_LT_k
2965                  { /* olvo 980921 changed order to allow x = erf(x) */
2966#if defined(_HIGHER_ORDER_)
2967                      zOP      = dp_z+i;
2968                      (*zOP--) = (i+1) * (*Targ1);
2969#endif /* _HIGHER_ORDER_ */
2970
2971                      TRES_FOINC = dp_T0[arg2] * TARG1_INC;
2972
2973#if defined(_HIGHER_ORDER_)
2974                      Targ2OP = Targ2;
2975
2976                      *Tres *= (i+1);
2977                      for (j=0;j<i;j++)
2978                      *Tres += (*Targ2OP++) * (*zOP--);
2979                      *Tres++ /= (i+1);
2980#endif /* _HIGHER_ORDER_ */
2981                  }
2982                  HOV_INC(Targ2, k)
2983                }
2984#endif
2985#endif
2986#endif /* ALL_TOGETHER_AGAIN */
2987                break;
2988
2989#endif
2990
2991            /*--------------------------------------------------------------------------*/
2992        case log_op:                                                /* log_op */
2993                arg = get_locint_f();
2994                res = get_locint_f();
2995
2996                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2997
2998#if defined(_INDO_)
2999#if defined(_INDOPRO_)
3000                copy_index_domain(res, arg, ind_dom);
3001#endif
3002#if defined(_NONLIND_)
3003                fod[opind].entry = maxopind+2;
3004                fod[opind].left = &fod[arg_index[arg]];
3005                fod[opind].right = NULL;
3006                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
3007                arg_index[res] = opind++;               
3008#endif
3009#else
3010#if !defined(_ZOS_) /* BREAK_ZOS */
3011                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
3012                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
3013
3014#ifdef _INT_FOR_
3015                FOR_0_LE_l_LT_p
3016                TRES_FOINC = TARG_INC;
3017#else
3018                divs = 1.0 / dp_T0[arg];
3019                FOR_0_LE_l_LT_p
3020                { if (dp_T0[arg] == 0.0) {
3021                  TargOP = Targ;
3022                  FOR_0_LE_i_LT_k
3023                  { if (*TargOP++ < 0.0) {
3024                        divs = make_nan();
3025                            BREAK_FOR_I
3026                        }
3027                      }
3028                  }
3029
3030                  /* olvo 980921 changed order to allow x = log(x) */
3031                  FOR_0_LE_i_LT_k
3032                  { TRES_FOINC = TARG_INC * divs;
3033#if defined(_HIGHER_ORDER_)
3034                    TresOP = Tres - i;
3035                    zOP    = dp_z+i;
3036
3037                    (*zOP--) = *Tres;
3038                    (*Tres) *= i+1;
3039                    for (j=0;j<i;j++)
3040                    (*Tres) -= (*zOP--) * (*TresOP++) * (j+1);
3041                    *Tres++ /= i+1;
3042#endif /* _HIGHER_ORDER_ */
3043                  }
3044                }
3045#endif
3046#endif
3047#endif /* ALL_TOGETHER_AGAIN */
3048#if !defined(_NTIGHT_)
3049                dp_T0[res] = log(dp_T0[arg]);
3050#endif /* !_NTIGHT_ */
3051
3052                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
3053                break;
3054
3055                /*--------------------------------------------------------------------------*/
3056            case pow_op:                                                /* pow_op */
3057                arg   = get_locint_f();
3058                res   = get_locint_f();
3059
3060                coval = get_val_f();
3061
3062                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3063
3064                /* olvo 980921 necessary for reverse */
3065                if (arg == res) {
3066                    IF_KEEP_WRITE_TAYLOR(arg,keep,k,p)
3067                }
3068
3069#if !defined(_NTIGHT_)
3070
3071#ifndef _ZOS_ /* BREAK_ZOS */
3072#if !defined(_INT_FOR_)
3073                T0arg   = dp_T0[arg];
3074#endif
3075#endif /* ALL_TOGETHER_AGAIN */
3076
3077                dp_T0[res] =
3078                    pow(dp_T0[arg], coval);
3079#endif /* !_NTIGHT_ */
3080
3081                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
3082#if defined(_INDO_)
3083#if defined(_INDOPRO_)
3084                copy_index_domain(res, arg, ind_dom);
3085#endif
3086#if defined(_NONLIND_)
3087                fod[opind].entry = maxopind+2;
3088                fod[opind].left = &fod[arg_index[arg]];
3089                fod[opind].right = NULL;
3090                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
3091                arg_index[res] = opind++;               
3092#endif
3093#else
3094#ifndef _ZOS_ /* BREAK_ZOS */
3095                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
3096                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
3097
3098#ifdef _INT_FOR_
3099                FOR_0_LE_l_LT_p
3100                TRES_FOINC = TARG_INC;
3101#else
3102                if (T0arg == 0.0) {
3103                    if (coval <= 0.0)
3104                        FOR_0_LE_l_LT_pk
3105                        TRES_INC = make_nan();
3106                    else {
3107                        /* coval not a whole number */
3108                        if (coval - floor(coval) != 0) {
3109                            FOR_0_LE_l_LT_p
3110                            {
3111                                i = 0;
3112                                FOR_0_LE_i_LT_k
3113                                {
3114                                    if (coval - i > 1)
3115                                    TRES_INC = 0;
3116                                    if ((coval - i < 1) && (coval - i > 0))
3117                                        TRES_INC = make_inf();
3118                                        if (coval - i < 0)
3119                                            TRES_INC = make_nan();
3120                                        }
3121                                    }
3122                                } else {
3123                        if (coval == 1) {
3124                                FOR_0_LE_l_LT_pk
3125                                TRES_INC = TARG_INC;
3126                            } else
3127                                /* coval is an int > 1 */
3128                                /* the following is not efficient but at least it works */
3129                                /* it reformulates x^n into x* ... *x n times */
3130                            {
3131                                INC_pk_1(Targ)
3132                                INC_pk_1(Tres)
3133
3134                                FOR_p_GT_l_GE_0
3135                                {
3136                                    FOR_k_GT_i_GE_0
3137                                    {
3138                                        *Tres = 0;
3139                                        DEC_TRES_FO
3140#if defined(_HIGHER_ORDER_)
3141                                        if (i == k-1) {
3142                                        zOP = dp_z+k-1;
3143                                        for(j=k-1;j>=0;j--) {
3144                                                (*zOP--) = (*Targ--);
3145                                            }
3146                                        }
3147                                        for (j=0;j<i;j++) {
3148                                        *Tres += dp_z[j] *
3149                                                     dp_z[i-j-1];
3150                                        }
3151                                        Tres--;
3152#endif /* _HIGHER_ORDER_ */
3153                                    }
3154                                }
3155                                for(ii=3;ii<=coval;ii++) {
3156                                    ASSIGN_T(Tres, TAYLOR_BUFFER[res])
3157                                    ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
3158                                    INC_pk_1(Targ)
3159                                    INC_pk_1(Tres)
3160
3161                                    FOR_p_GT_l_GE_0
3162                                    {
3163                                        FOR_k_GT_i_GE_0
3164                                        {
3165                                            *Tres = 0;
3166                                            DEC_TRES_FO
3167#if defined(_HIGHER_ORDER_)
3168                                            TresOP = Tres-i;
3169                                            for (j=0;j<i;j++)
3170                                            *Tres += TresOP[j] * dp_z[i-j-1];
3171                                            Tres--;
3172#endif /* _HIGHER_ORDER_ */
3173                                        }
3174                                    }
3175                                }
3176                        }
3177                    }
3178                }
3179            } else {
3180                r0 = 1.0 / T0arg;
3181                FOR_0_LE_l_LT_p
3182                FOR_0_LE_i_LT_k
3183                { /* olvo 980921 changed order to allow x = pow(x,n) */
3184#ifdef _HIGHER_ORDER_
3185                    zOP      = dp_z+i;
3186                    (*zOP--) = (*Targ) * r0;
3187#endif /* _HIGHER_ORDER_ */
3188
3189                    TRES_FOINC = dp_T0[res] *
3190                                 TARG_INC * coval * r0;
3191
3192#ifdef _HIGHER_ORDER_
3193                    TresOP = Tres-i;
3194
3195                    (*Tres) *= i+1;
3196                    y = coval*i -1;
3197                    for (j=0;j<i;j++) {
3198                        *Tres += (*TresOP++) * (*zOP--) * y;
3199                            y -= coval + 1;
3200                        }
3201                        *Tres++ /= (i+1);
3202#endif /* _HIGHER_ORDER_ */
3203                    }
3204                }
3205#endif
3206#endif
3207#endif /* ALL_TOGETHER_AGAIN */
3208                break;
3209
3210                /*--------------------------------------------------------------------------*/
3211            case sqrt_op:                                              /* sqrt_op */
3212                arg = get_locint_f();
3213                res = get_locint_f();
3214
3215                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3216
3217#if !defined(_NTIGHT_)
3218                dp_T0[res] = sqrt(dp_T0[arg]);
3219#endif /* !_NTIGHT_ */
3220
3221                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
3222
3223#if defined(_INDO_)
3224#if defined(_INDOPRO_)
3225                copy_index_domain(res, arg, ind_dom);
3226#endif
3227#if defined(_NONLIND_)
3228                fod[opind].entry = maxopind+2;
3229                fod[opind].left = &fod[arg_index[arg]];
3230                fod[opind].right = NULL;
3231                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
3232                arg_index[res] = opind++;               
3233#endif
3234#else
3235#if !defined(_ZOS_) /* BREAK_ZOS */
3236                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
3237                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
3238
3239#ifdef _INT_FOR_
3240                FOR_0_LE_l_LT_p
3241                TRES_FOINC = TARG_INC;
3242#else
3243                FOR_0_LE_l_LT_p
3244                { TargOP = Targ;
3245                  if (dp_T0[arg] == 0.0)
3246                  /* Note: <=> dp_T0[res] == 0.0 */
3247              { r0 = 0.0;
3248                  FOR_0_LE_i_LT_k
3249                  { if (TARG>0.0) {
3250                        r0 = make_inf();
3251                            VEC_INC(Targ, k-i)
3252                            BREAK_FOR_I
3253                        } else
3254                            if (TARG<0.0) {
3255                            r0 = make_nan();
3256                                VEC_INC(Targ, k-i)
3257                                BREAK_FOR_I
3258                            } else
3259                                Targ++;
3260                              }
3261                          }
3262                  else {
3263                      r0 = 0.5/dp_T0[res];
3264                  }
3265                  Targ = TargOP;
3266
3267                  even = 1;
3268                  FOR_0_LE_i_LT_k
3269                  { TRES_FOINC = r0 * TARG_INC;
3270#if defined(_HIGHER_ORDER_)
3271                    TresOP  = Tres-i;
3272                    TresOP2 = Tres-1;
3273
3274                    x = 0;
3275                    for (j=1;2*j-1<i;j++)
3276                    x += (*TresOP++) * (*TresOP2--);
3277                    x *= 2;
3278                    if (!even)
3279                        x += (*TresOP) * (*TresOP2); /* !!! */
3280                        even = !even;
3281                        *Tres++ -= r0*x;
3282#endif /* _HIGHER_ORDER_ */
3283                      }
3284                    }
3285#endif
3286#endif
3287#endif /* ALL_TOGETHER_AGAIN */
3288                    break;
3289
3290            /*--------------------------------------------------------------------------*/
3291        case gen_quad:                                            /* gen_quad */
3292            arg1 = get_locint_f();
3293                arg2 = get_locint_f();
3294                res  = get_locint_f();
3295
3296#if !defined(_NTIGHT_)
3297                if (get_val_f()!=dp_T0[arg1]) {
3298                    fprintf(DIAG_OUT,
3299                            "ADOL-C Warning: forward sweep aborted; tape invalid!\n");
3300                    IF_KEEP_TAYLOR_CLOSE
3301                    end_sweep();
3302                    return -2;
3303                }
3304#endif /* !_NTIGHT_ */
3305
3306                coval = get_val_f();
3307
3308                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3309
3310#if !defined(_NTIGHT_)
3311                dp_T0[res] = coval;
3312#endif /* !_NTIGHT_ */
3313
3314#if defined(_INDO_)
3315               fprintf(DIAG_OUT,
3316                    "ADOL-C Warning: forward sweep aborted; sparse mode not available for gen_quad!\n");
3317               end_sweep();
3318               return -2;
3319#else
3320#if !defined(_ZOS_) /* BREAK_ZOS */
3321                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
3322                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
3323#ifdef _INT_FOR_
3324                FOR_0_LE_l_LT_p
3325                TRES_FOINC = TARG1_FOINC;
3326#else
3327                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
3328
3329                FOR_0_LE_l_LT_p
3330                { FOR_0_LE_i_LT_k
3331                  { /* olvo 980922 changed order to allow x = gen_quad(x) */
3332#if defined(_HIGHER_ORDER_)
3333                      zOP      = dp_z+i;
3334                      (*zOP--) = (i+1) * (*Targ1);
3335#endif /* _HIGHER_ORDER_ */
3336
3337                      TRES_FOINC = dp_T0[arg2] * TARG1_INC;
3338
3339#if defined(_HIGHER_ORDER_)
3340                      Targ2OP = Targ2;
3341
3342                      *Tres *= (i+1);
3343                      for (j=0;j<i;j++)
3344                      *Tres += (*Targ2OP++) * (*zOP--);
3345                      *Tres++ /= (i+1);
3346#endif /* _HIGHER_ORDER_ */
3347                  }
3348                  HOV_INC(Targ2, k)
3349                }
3350#endif
3351#endif
3352#endif /* ALL_TOGETHER_AGAIN */
3353                break;
3354
3355            /*--------------------------------------------------------------------------*/
3356        case min_op:                                                /* min_op */
3357            arg1  = get_locint_f();
3358                arg2  = get_locint_f();
3359                res   = get_locint_f();
3360                coval = get_val_f();
3361
3362                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3363
3364#if !defined(_NTIGHT_)
3365                /* olvo 980923 changed order to allow x = min(x,y) etc. */
3366
3367                /* olvo/mitev 980721 return value (taken from below) */
3368                if (dp_T0[arg1] > dp_T0[arg2]) {
3369                    if (coval)
3370                        MINDEC(ret_c,2);
3371                } else
3372                    if (dp_T0[arg1] < dp_T0[arg2]) {
3373                        if (!coval)
3374                            MINDEC(ret_c,2);
3375                    } else
3376                        if (arg1 != arg2)
3377                            MINDEC(ret_c,1);
3378#endif /* !_NTIGHT_ */
3379
3380#if defined (_INDO_)
3381#if defined (_INDOPRO_)
3382#if defined (_TIGHT_)
3383                    if (dp_T0[arg1] < dp_T0[arg2])
3384                        copy_index_domain(res, arg1, ind_dom);
3385                    else {
3386                        if (dp_T0[arg1] > dp_T0[arg2])
3387                            copy_index_domain(res, arg2, ind_dom);
3388                        else
3389                            combine_2_index_domains(res, arg1, arg2, ind_dom);
3390                    }
3391#else
3392                    combine_2_index_domains(res, arg1, arg2, ind_dom);
3393#endif
3394#endif
3395#if defined(_NONLIND_)
3396#ifdef _TIGHT_
3397                    if (dp_T0[arg1] < dp_T0[arg2])
3398                      {
3399                        fod[opind].entry = maxopind+2;
3400                        fod[opind].left = &fod[arg_index[arg1]];
3401                        fod[opind].right = NULL;
3402                        arg_index[res] = opind++;               
3403                      }           
3404                    else {
3405                        if (dp_T0[arg1] > dp_T0[arg2])
3406                          {
3407                            fod[opind].entry = maxopind+2;
3408                            fod[opind].left = &fod[arg_index[arg2]];
3409                            fod[opind].right = NULL;
3410                            arg_index[res] = opind++;           
3411
3412                          }               
3413                        else
3414                          {
3415                            fod[opind].entry = maxopind+2;
3416                            fod[opind].left = &fod[arg_index[arg1]];
3417                            fod[opind].right = &fod[arg_index[arg2]];
3418                            arg_index[res] = opind++;           
3419                          }
3420                    }
3421#else
3422                    fod[opind].entry = maxopind+2;
3423                    fod[opind].left = &fod[arg_index[arg1]];
3424                    fod[opind].right = &fod[arg_index[arg2]];
3425                    arg_index[res] = opind++;           
3426                    arg_index[res] = opind++;           
3427#endif
3428#endif
3429#else
3430#if !defined(_ZOS_) /* BREAK_ZOS */
3431                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
3432                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
3433                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
3434
3435#ifdef _INT_FOR_
3436#ifdef _TIGHT_
3437                Tqo = NULL;
3438                if (dp_T0[arg1] > dp_T0[arg2])
3439                    Tqo = Targ2;
3440                else
3441                    if (dp_T0[arg1] < dp_T0[arg2])
3442                        Tqo = Targ1;
3443
3444                FOR_0_LE_l_LT_p
3445                { Targ = Tqo;
3446                  if (Targ == NULL) /* e.g. T0[arg1] == T0[arg2] */
3447                    { Targ1OP = Targ1;
3448                      Targ2OP = Targ2;
3449                      if (TARG1 > TARG2)
3450                          Targ = Targ2OP;
3451                      else
3452                          if (TARG1 < TARG2)
3453                              Targ = Targ1OP;
3454                      Targ1++;
3455                      Targ2++;
3456                      if (Targ == NULL) /* e.g. both are equal */
3457                          Targ = Targ1OP;
3458                  }
3459
3460                  TRES_INC = TARG_INC;
3461
3462                  if (Tqo)
3463                    Tqo++;
3464                }
3465
3466                dp_T0[res] = MIN_ADOLC(dp_T0[arg1], dp_T0[arg2]);
3467#endif /* _TIGHT_ */
3468#ifdef _NTIGHT_
3469                TRES_INC = TARG1_INC | TARG2_INC;
3470#endif /* _NTIGHT_ */
3471#else
3472                Tqo = NULL;
3473                if (dp_T0[arg1] > dp_T0[arg2])
3474                    Tqo = Targ2;
3475                else
3476                    if (dp_T0[arg1] < dp_T0[arg2])
3477                        Tqo = Targ1;
3478
3479                FOR_0_LE_l_LT_p
3480                { Targ = Tqo;
3481                  if (Targ == NULL) /* e.g. dp_T0[arg1] ==
3482                                                                                 dp_T0[arg2] */
3483              { Targ1OP = Targ1;
3484                  Targ2OP = Targ2;
3485                  FOR_0_LE_i_LT_k
3486                  { if (TARG1 > TARG2) {
3487                        Targ = Targ2OP;
3488                        VEC_INC(Targ1, k-i)
3489                            VEC_INC(Targ2, k-i)
3490                            BREAK_FOR_I
3491                        } else
3492                            if (TARG1 < TARG2) {
3493                            Targ = Targ1OP;
3494                            VEC_INC(Targ1, k-i)
3495                                VEC_INC(Targ2, k-i)
3496                                BREAK_FOR_I
3497                            }
3498                        Targ1++;
3499                        Targ2++;
3500                      }
3501                      if (Targ == NULL) /* e.g. both are equal */
3502                          Targ = Targ1OP;
3503                  }
3504
3505                  FOR_0_LE_i_LT_k
3506                  TRES_INC = TARG_INC;
3507
3508                  if (Tqo) {
3509                  VEC_INC(Tqo, k)
3510                  }
3511            }
3512#endif
3513#endif
3514#endif /* ALL_TOGETHER_AGAIN */
3515#if !defined(_NTIGHT_)
3516                dp_T0[res] =
3517                    MIN_ADOLC( dp_T0[arg1],
3518                               dp_T0[arg2] );
3519#endif /* !_NTIGHT_ */
3520                break;
3521
3522                /*--------------------------------------------------------------------------*/
3523            case abs_val:                                              /* abs_val */
3524                arg   = get_locint_f();
3525                res   = get_locint_f();
3526                coval = get_val_f();
3527
3528                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3529
3530#if !defined(_NTIGHT_)
3531                /* olvo 980923 changed order to allow x = min(x,y) etc. */
3532
3533                /* olvo/mitev 980721 ec n3l (taken from below) */
3534                if (dp_T0[arg] < 0.0) {
3535                    if (coval)
3536                        MINDEC(ret_c,2);
3537                } else
3538                    if (dp_T0[arg] > 0.0) {
3539                        if (!coval)
3540                            MINDEC(ret_c,2);
3541                    }
3542#endif /* !_NTIGHT_ */
3543
3544#if defined(_INDO_)
3545#if defined(_INDOPRO_)
3546                copy_index_domain(res, arg, ind_dom);
3547#endif
3548#if defined(_NONLIND_)
3549                arg_index[res] = arg_index[arg];               
3550#endif
3551#else
3552#if !defined(_ZOS_) /* BREAK_ZOS */
3553                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
3554                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
3555
3556#ifdef _INT_FOR_
3557#ifdef _TIGHT_
3558                y = 0.0;
3559                if (dp_T0[arg] != 0.0) {
3560                    if (dp_T0[arg] < 0.0)
3561                        y = -1.0;
3562                    else
3563                        y = 1.0;
3564                }
3565                FOR_0_LE_l_LT_p
3566                { if ((y == 0.0) && (TARG != 0.0))
3567                  MINDEC(ret_c,1);
3568
3569                  TRES_INC = TARG_INC;
3570                }
3571
3572                dp_T0[res] = fabs(dp_T0[arg]);
3573#endif /* _TIGHT_ */
3574#ifdef _NTIGHT_
3575                FOR_0_LE_l_LT_p
3576                TRES_INC = TARG_INC;
3577#endif /* _NTIGHT_ */
3578#else
3579                y = 0.0;
3580                if (dp_T0[arg] != 0.0) {
3581                    if (dp_T0[arg] < 0.0)
3582                        y = -1.0;
3583                    else
3584                        y = 1.0;
3585                }
3586
3587                FOR_0_LE_l_LT_p
3588                { x = y;
3589                  FOR_0_LE_i_LT_k
3590                  { if ((x == 0.0) && (TARG != 0.0)) {
3591                    MINDEC(ret_c,1);
3592                        if (TARG < 0.0)
3593                            x = -1.0;
3594                        else
3595                            x = 1.0;
3596                    }
3597                    TRES_INC = x * TARG_INC;
3598              }
3599            }
3600#endif
3601#endif
3602#endif /* ALL_TOGETHER_AGAIN */
3603#if !defined(_NTIGHT_)
3604                dp_T0[res] = fabs(dp_T0[arg]);
3605#endif /* !_NTIGHT_ */
3606                break;
3607
3608                /*--------------------------------------------------------------------------*/
3609            case ceil_op:                                              /* ceil_op */
3610                arg   = get_locint_f();
3611                res   = get_locint_f();
3612                coval = get_val_f();
3613
3614                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3615
3616#if !defined(_NTIGHT_)
3617                dp_T0[res]=ceil(dp_T0[arg]);
3618                /* olvo/mitev 980721 ec n2l (taken from below) */
3619                if (coval != dp_T0[res])
3620                    MINDEC(ret_c,2);
3621#endif /* !_NTIGHT_ */
3622
3623#if defined(_INDO_)
3624#if defined(_INDOPRO_)
3625#ifdef _TIGHT_
3626                ind_dom[res][0] = 0;
3627#else
3628                copy_index_domain(res, arg, ind_dom);
3629#endif /* _TIGHT_ */
3630#endif
3631#if defined(_NONLIND_)
3632                arg_index[res] = arg_index[arg];               
3633#endif
3634#else
3635#if !defined(_ZOS_) /* BREAK_ZOS */
3636                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
3637
3638                FOR_0_LE_l_LT_pk
3639                TRES_INC = 0.0;
3640#endif
3641#endif /* ALL_TOGETHER_AGAIN */
3642                break;
3643
3644                /*--------------------------------------------------------------------------*/
3645            case floor_op:                 /* Compute ceil of adouble    floor_op */
3646                arg   = get_locint_f();
3647                res   = get_locint_f();
3648                coval = get_val_f();
3649
3650                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3651
3652#if !defined(_NTIGHT_)
3653                dp_T0[res] = floor(dp_T0[arg]);
3654                /* olvo/mitev 980721 ec n2l (taken from below) */
3655                if (coval != dp_T0[res])
3656                    MINDEC(ret_c,2);
3657#endif /* !_NTIGHT_ */
3658
3659#if defined(_INDO_)
3660#if defined(_INDOPRO_)
3661#ifdef _TIGHT_
3662                ind_dom[res][0] = 0;
3663#else
3664                copy_index_domain(res, arg, ind_dom);
3665#endif /* _TIGHT_ */
3666#endif
3667#if defined(_NONLIND_)
3668                arg_index[res] = arg_index[arg];               
3669#endif
3670#else
3671#if !defined(_ZOS_) /* BREAK_ZOS */
3672                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
3673
3674                FOR_0_LE_l_LT_pk
3675                TRES_INC = 0.0;
3676#endif
3677#endif /* ALL_TOGETHER_AGAIN */
3678                break;
3679
3680
3681                /****************************************************************************/
3682                /*                                                             CONDITIONALS */
3683
3684                /*--------------------------------------------------------------------------*/
3685            case cond_assign:                                      /* cond_assign */
3686                arg   = get_locint_f();
3687                arg1  = get_locint_f();
3688                arg2  = get_locint_f();
3689                res   = get_locint_f();
3690                coval = get_val_f();
3691
3692                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3693
3694                /* olvo 980924 changed order to allow reflexive ops */
3695#if defined (_INDO_)
3696#if defined (_INDOPRO_)
3697#if defined (_TIGHT_)
3698                if (dp_T0[arg] > 0) {
3699                    if (coval <= 0.0)
3700                        MINDEC(ret_c,2);
3701                    dp_T0[res] = dp_T0[arg1];
3702
3703                    copy_index_domain(res, arg1, ind_dom);
3704
3705                } else {
3706                    if (coval > 0.0)
3707                        MINDEC(ret_c,2);
3708                    if (dp_T0[arg] == 0)
3709                        MINDEC(ret_c,0);
3710                    dp_T0[res] = dp_T0[arg2];
3711                    copy_index_domain(res, arg2, ind_dom);
3712                }
3713#else
3714                    combine_2_index_domains(res, arg1, arg2, ind_dom);
3715#endif
3716#endif
3717#if defined (_NONLIND_)
3718#ifdef _TIGHT_
3719                if (dp_T0[arg] > 0) {
3720                    if (coval <= 0.0)
3721                        MINDEC(ret_c,2);
3722                    dp_T0[res] = dp_T0[arg1];
3723
3724                    arg_index[res] = arg_index[arg1];           
3725                } else {
3726                    if (coval > 0.0)
3727                        MINDEC(ret_c,2);
3728                    if (dp_T0[arg] == 0)
3729                        MINDEC(ret_c,0);
3730                    dp_T0[res] = dp_T0[arg2];
3731
3732                    arg_index[res] = arg_index[arg2];           
3733                }
3734
3735#else
3736               arg_index[res] = opind++;               
3737#endif
3738#endif
3739#else
3740#if !defined(_ZOS_) /* BREAK_ZOS */
3741                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
3742                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
3743                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
3744#endif /* ALL_TOGETHER_AGAIN */
3745
3746#ifdef _INT_FOR_
3747#ifdef _TIGHT_
3748                coval = get_val_f();
3749
3750                if (dp_T0[arg] > 0)
3751                    FOR_0_LE_l_LT_pk
3752                    TRES_INC = TARG1_INC;
3753                else
3754                    FOR_0_LE_l_LT_pk
3755                    TRES_INC = TARG2_INC;
3756
3757                if (dp_T0[arg] > 0) {
3758                    if (coval <= 0.0)
3759                        MINDEC(ret_c,2);
3760                    dp_T0[res] = dp_T0[arg1];
3761                } else {
3762                    if (coval > 0.0)
3763                        MINDEC(ret_c,2);
3764                    if (dp_T0[arg] == 0)
3765                        MINDEC(ret_c,0);
3766                    dp_T0[res] = dp_T0[arg2];
3767                }
3768#endif /* _TIGHT_ */
3769#ifdef _NTIGHT_
3770                FOR_0_LE_l_LT_pk
3771                TRES_INC = TARG1_INC | TARG2_INC;
3772#endif /* _NTIGHT_ */
3773#else
3774#if !defined(_ZOS_) /* BREAK_ZOS */
3775                if (dp_T0[arg] > 0)
3776                    FOR_0_LE_l_LT_pk
3777                    TRES_INC = TARG1_INC;
3778                else
3779                    FOR_0_LE_l_LT_pk
3780                    TRES_INC = TARG2_INC;
3781#endif
3782
3783                if (dp_T0[arg] > 0) {
3784                    if (coval <= 0.0)
3785                        MINDEC(ret_c,2);
3786                    dp_T0[res] = dp_T0[arg1];
3787                } else {
3788                    if (coval > 0.0)
3789                        MINDEC(ret_c,2);
3790                    if (dp_T0[arg] == 0)
3791                        MINDEC(ret_c,0);
3792                    dp_T0[res] = dp_T0[arg2];
3793                }
3794#endif
3795#endif /* ALL_TOGETHER_AGAIN */
3796                break;
3797
3798                /*--------------------------------------------------------------------------*/
3799            case cond_assign_s:                                  /* cond_assign_s */
3800                arg   = get_locint_f();
3801                arg1  = get_locint_f();
3802                res   = get_locint_f();
3803                coval = get_val_f();
3804
3805                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3806
3807                /* olvo 980924 changed order to allow reflexive ops */
3808#if defined(_INDO_)
3809#if defined(_INDOPRO_)
3810                copy_index_domain(res, arg1, ind_dom);
3811#endif
3812#if defined(_NONLIND_)
3813                arg_index[res] = arg_index[arg1];               
3814#endif
3815#else
3816#if !defined(_ZOS_) /* BREAK_ZOS */
3817                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
3818                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
3819#endif /* ALL_TOGETHER_AGAIN */
3820
3821#ifdef _INT_FOR_
3822#ifdef _TIGHT_
3823                coval = get_val_f();
3824
3825                if (dp_T0[arg] > 0)
3826#endif /* _TIGHT_ */
3827                    FOR_0_LE_l_LT_pk
3828                    TRES_INC = TARG1_INC;
3829
3830#ifdef _TIGHT_
3831                if (dp_T0[arg] > 0) {
3832                    if (coval <= 0.0)
3833                        MINDEC(ret_c,2);
3834                    dp_T0[res] = dp_T0[arg1];
3835                } else
3836                    if (dp_T0[arg] == 0)
3837                        MINDEC(ret_c,0);
3838#endif /* _TIGHT_ */
3839#else
3840#if !defined(_ZOS_) /* BREAK_ZOS */
3841                if (dp_T0[arg] > 0)
3842                    FOR_0_LE_l_LT_pk
3843                    TRES_INC = TARG1_INC;
3844#endif
3845                if (dp_T0[arg] > 0) {
3846                    if (coval <= 0.0)
3847                        MINDEC(ret_c,2);
3848                    dp_T0[res] = dp_T0[arg1];
3849                } else
3850                    if (dp_T0[arg] == 0)
3851                        MINDEC(ret_c,0);
3852#endif
3853#endif /* ALL_TOGETHER_AGAIN */
3854                break;
3855
3856
3857                /*--------------------------------------------------------------------------*/
3858                /* NEW CONDITIONALS */
3859                /*--------------------------------------------------------------------------*/
3860            case neq_a_a:
3861            case eq_a_a:
3862            case le_a_a:
3863            case ge_a_a:
3864            case lt_a_a:
3865            case gt_a_a:
3866                coval = get_val_f();
3867                arg = get_locint_f();
3868                arg1 = get_locint_f();
3869                res = get_locint_f();
3870#if !defined(_NTIGHT_)
3871                {
3872                    revreal retval = -1;
3873                    const char* opname = "";
3874                    switch (operation) {
3875                    case neq_a_a:
3876                        retval = (revreal)(dp_T0[arg] != dp_T0[arg1]);
3877                        opname = "neq_a_a";
3878                        break;
3879                    case eq_a_a:
3880                        retval = (revreal)(dp_T0[arg] == dp_T0[arg1]);
3881                        opname = "eq_a_a";
3882                        break;
3883                    case ge_a_a:
3884                        retval = (revreal)(dp_T0[arg] >= dp_T0[arg1]);
3885                        opname = "ge_a_a";
3886                        break;
3887                    case le_a_a:
3888                        retval = (revreal)(dp_T0[arg] <= dp_T0[arg1]);
3889                        opname = "le_a_a";
3890                        break;
3891                    case gt_a_a:
3892                        retval = (revreal)(dp_T0[arg] > dp_T0[arg1]);
3893                        opname = "gt_a_a";
3894                        break;
3895                    case lt_a_a:
3896                        retval = (revreal)(dp_T0[arg] < dp_T0[arg1]);
3897                        opname = "lt_a_a";
3898                        break;
3899                    }
3900                    if (retval != coval && ADOLC_GLOBAL_TAPE_VARS.branchSwitchWarning)
3901                        fprintf(DIAG_OUT,
3902                                "ADOL-C Warning: Branch switch detected in comparison "
3903                                "(operator %s).\n"
3904                                "Results may be unpredictable! Retaping recommended!\n",opname);
3905                    dp_T0[res] = retval;
3906                }
3907#endif
3908#if defined(_INDO_)
3909#if defined(_INDOPRO_)
3910                ind_dom[res][0]=0;
3911#endif
3912#if defined(_NONLIND_)
3913                fod[opind].entry = maxopind+2;
3914                fod[opind].left = NULL;
3915                fod[opind].right = NULL;
3916                arg_index[res] = opind++;
3917#endif
3918#else
3919#if !defined(_ZOS_) /* BREAK_ZOS */
3920                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
3921
3922                FOR_0_LE_l_LT_pk
3923                TRES_INC = 0;
3924#endif
3925#endif /* ALL_TOGETHER_AGAIN */
3926
3927                break;
3928
3929                /*--------------------------------------------------------------------------*/
3930            case subscript:
3931                coval = get_val_f();
3932                arg = get_locint_f();
3933                res = get_locint_f();
3934                {
3935                    size_t cnt, idx, numvar = (size_t)trunc(fabs(coval));
3936                    locint vectorloc;
3937                    vectorloc = get_locint_f();
3938#if !defined(_NTIGHT_)
3939                    idx = (size_t)trunc(fabs(dp_T0[arg]));
3940                    if (idx >= numvar)
3941                        fprintf(DIAG_OUT, "ADOL-C warning: index out of bounds while subscripting n=%z, idx=%z\n", numvar, idx);
3942                    arg1 = vectorloc+idx;
3943                    IF_KEEP_WRITE_TAYLOR(res,keep,k,p);
3944                    dp_T0[res] = dp_T0[arg1];
3945#if defined(_INDO_)
3946                    //copy_index_domain(res, arg1, ind_dom);
3947#else
3948#if !defined(_ZOS_) /* BREAK_ZOS */
3949                    ASSIGN_T(Targ1,TAYLOR_BUFFER[arg1])
3950                    ASSIGN_T(Tres,TAYLOR_BUFFER[res])
3951
3952                    FOR_0_LE_l_LT_pk
3953                    TRES_INC = TARG1_INC;
3954#endif
3955#endif
3956#else
3957                    fprintf(DIAG_OUT, "ADOL-C error: active subscripting does not work in safe mode, please use tight mode\n");
3958                    exit(-2);
3959#endif /* ALL_TOGETHER_AGAIN */
3960                }
3961                break;
3962
3963            case subscript_ref:
3964                coval = get_val_f();
3965                arg = get_locint_f();
3966                res = get_locint_f();
3967                {
3968                    size_t cnt, idx, numvar = (size_t)trunc(fabs(coval));
3969                    locint vectorloc;
3970                    vectorloc = get_locint_f();
3971#if !defined(_NTIGHT_)
3972                    idx = (size_t)trunc(fabs(dp_T0[arg]));
3973                    if (idx >= numvar)
3974                        fprintf(DIAG_OUT, "ADOL-C warning: index out of bounds while subscripting (ref) n=%z, idx=%z\n", numvar, idx);
3975                    arg1 = vectorloc+idx;
3976                    IF_KEEP_WRITE_TAYLOR(res,keep,k,p);
3977                    dp_T0[res] = arg1;
3978#else
3979                    fprintf(DIAG_OUT, "ADOL-C error: active subscripting does not work in safe mode, please use tight mode\n");
3980                    exit(-2);
3981#endif
3982                }
3983                break;
3984
3985            case ref_copyout:
3986                arg = get_locint_f();
3987                res = get_locint_f();
3988#if !defined(_NTIGHT_)
3989                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
3990                IF_KEEP_WRITE_TAYLOR(res,keep,k,p);
3991                dp_T0[res] = dp_T0[arg1];
3992#if defined(_INDO_)
3993                //copy_index_domain(res, arg1, ind_dom);
3994#else
3995#if !defined(_ZOS_) /* BREAK_ZOS */
3996                ASSIGN_T(Targ1,TAYLOR_BUFFER[arg1])
3997                ASSIGN_T(Tres,TAYLOR_BUFFER[res])
3998
3999                FOR_0_LE_l_LT_pk
4000                TRES_INC = TARG1_INC;
4001#endif
4002#endif
4003#else
4004                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4005                exit(-2);
4006#endif /* ALL_TOGETHER_AGAIN */
4007                break;
4008
4009            case ref_incr_a:
4010                arg = get_locint_f();
4011#if !defined(_NTIGHT_)
4012                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
4013                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p);
4014                dp_T0[arg1]++;
4015#else
4016                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4017                exit(-2);
4018#endif
4019                break;
4020
4021            case ref_decr_a:
4022                arg = get_locint_f();
4023#if !defined(_NTIGHT_)
4024                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
4025                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p);
4026                dp_T0[arg1]--;
4027#else
4028                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4029                exit(-2);
4030#endif
4031                break;
4032
4033            case ref_assign_d:
4034                arg = get_locint_f();
4035                coval = get_val_f();
4036               
4037#if !defined(_NTIGHT_)
4038                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
4039                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p)
4040                dp_T0[arg1] = coval;
4041#if defined(_INDO_)
4042                //ind_dom[arg1][0] = 0;
4043#else
4044#if !defined(_ZOS_)
4045                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
4046
4047                FOR_0_LE_l_LT_pk
4048                TARG1_INC = 0;
4049#endif
4050#endif
4051#else
4052                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4053                exit(-2);
4054#endif
4055                break;
4056
4057            case ref_assign_d_zero:
4058                arg = get_locint_f();
4059
4060#if !defined(_NTIGHT_)
4061                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
4062                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p)
4063                dp_T0[arg1] = 0.0;
4064#if defined(_INDO_)
4065                //ind_dom[arg1][0] = 0;
4066#else
4067#if !defined(_ZOS_)
4068                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
4069
4070                FOR_0_LE_l_LT_pk
4071                TARG1_INC = 0;
4072#endif
4073#endif
4074#else
4075                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4076                exit(-2);
4077#endif
4078                break;
4079
4080            case ref_assign_d_one:
4081                arg = get_locint_f();
4082
4083#if !defined(_NTIGHT_)
4084                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
4085                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p)
4086                dp_T0[arg1] = 1.0;
4087#if defined(_INDO_)
4088                //ind_dom[arg1][0] = 0;
4089#else
4090#if !defined(_ZOS_)
4091                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
4092
4093                FOR_0_LE_l_LT_pk
4094                TARG1_INC = 0;
4095#endif
4096#endif
4097#else
4098                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4099                exit(-2);
4100#endif
4101                break;
4102
4103            case ref_assign_a:           /* assign an adouble variable an    assign_a */
4104                /* adouble value. (=) */
4105                arg = get_locint_f();
4106                res = get_locint_f();
4107
4108#if !defined(_NTIGHT_)
4109                arg1 = (size_t)trunc(fabs(dp_T0[res]));
4110                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p)
4111                dp_T0[arg1] = dp_T0[arg];
4112#if defined(_INDO_)
4113                //copy_index_domain(arg1, arg, ind_dom);
4114#else
4115#if !defined(_ZOS_) /* BREAK_ZOS */
4116                ASSIGN_T(Targ,TAYLOR_BUFFER[arg])
4117                ASSIGN_T(Targ1,TAYLOR_BUFFER[arg1])
4118
4119                FOR_0_LE_l_LT_pk
4120                TARG1_INC = TARG_INC;
4121#endif
4122#endif
4123#else
4124                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4125                exit(-2);
4126#endif /* ALL_TOGETHER_AGAIN */
4127                break;
4128
4129            case ref_assign_ind:       /* assign an adouble variable an    assign_ind */
4130                /* independent double value (<<=) */
4131                arg = get_locint_f();
4132
4133
4134#if !defined(_NTIGHT_)
4135                res = (size_t)trunc(fabs(dp_T0[arg]));
4136                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4137                dp_T0[res] = basepoint[indexi];
4138#if defined(_INDO_)
4139                //ind_dom[res][0] = 1;
4140                //ind_dom[res][2] = indexi;
4141#else
4142#if !defined(_ZOS_) /* BREAK_ZOS */
4143                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
4144
4145#ifdef _INT_FOR_
4146                FOR_0_LE_l_LT_p
4147                TRES_INC = ARGUMENT(indexi,l,i);
4148#else
4149                FOR_0_LE_l_LT_p
4150                FOR_0_LE_i_LT_k
4151                TRES_INC = ARGUMENT(indexi,l,i);
4152#endif
4153#endif
4154#endif
4155#else
4156                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4157                exit(-2);
4158#endif /* ALL_TOGETHER_AGAIN */
4159                ++indexi;
4160                break;
4161
4162            case ref_eq_plus_d:            /* Add a floating point to an    eq_plus_d */
4163                /* adouble. (+=) */
4164                arg  = get_locint_f();
4165                coval = get_val_f();
4166
4167
4168#if !defined(_NTIGHT_)
4169                res = (size_t)trunc(fabs(dp_T0[arg]));
4170                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4171                dp_T0[res] += coval;
4172#else
4173                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4174                exit(-2);
4175#endif /* !_NTIGHT_ */
4176                break;
4177
4178                /*--------------------------------------------------------------------------*/
4179            case ref_eq_plus_a:             /* Add an adouble to another    eq_plus_a */
4180                /* adouble. (+=) */
4181                arg = get_locint_f();
4182                arg1 = get_locint_f();
4183
4184#if !defined(_NTIGHT_)
4185                res = (size_t)trunc(fabs(dp_T0[arg1]));
4186                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4187                dp_T0[res] += dp_T0[arg];
4188#if defined(_INDO_)
4189                //merge_2_index_domains(res, arg, ind_dom);
4190#else
4191#if !defined(_ZOS_) /* BREAK_ZOS */
4192                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
4193                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
4194
4195#ifdef _INT_FOR_
4196                FOR_0_LE_l_LT_pk
4197                TRES_INC |= TARG_INC;
4198#else
4199                FOR_0_LE_l_LT_pk
4200                TRES_INC += TARG_INC;
4201#endif
4202#endif
4203#endif
4204#else
4205                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4206                exit(-2);
4207#endif /* ALL_TOGETHER_AGAIN */
4208                break;
4209
4210            case ref_eq_min_d:       /* Subtract a floating point from an    eq_min_d */
4211                /* adouble. (-=) */
4212                arg = get_locint_f();
4213                coval = get_val_f();
4214
4215#if !defined(_NTIGHT_)
4216                res = (size_t)trunc(fabs(dp_T0[arg]));
4217                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4218                dp_T0[res] -= coval;
4219#else
4220                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4221                exit(-2);
4222#endif /* !_NTIGHT_ */
4223                break;
4224
4225                /*--------------------------------------------------------------------------*/
4226            case ref_eq_min_a:        /* Subtract an adouble from another    eq_min_a */
4227                /* adouble. (-=) */
4228                arg = get_locint_f();
4229                arg1 = get_locint_f();
4230
4231#if !defined(_NTIGHT_)
4232                res = (size_t)trunc(fabs(dp_T0[arg1]));
4233                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4234                dp_T0[res] -= dp_T0[arg];
4235#if defined(_INDO_)
4236                //merge_2_index_domains(res, arg, ind_dom);
4237#else
4238#if !defined(_ZOS_) /* BREAK_ZOS */
4239                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
4240                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
4241
4242#ifdef _INT_FOR_
4243                FOR_0_LE_l_LT_pk
4244                TRES_INC |= TARG_INC;
4245#else
4246                FOR_0_LE_l_LT_pk
4247                TRES_INC -= TARG_INC;
4248#endif
4249#endif
4250#endif
4251#else
4252                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4253                exit(-2);
4254#endif /* ALL_TOGETHER_AGAIN */
4255                break;
4256
4257            case ref_eq_mult_d:              /* Multiply an adouble by a    eq_mult_d */
4258                /* flaoting point. (*=) */
4259                arg = get_locint_f();
4260                coval = get_val_f();
4261
4262#if !defined(_NTIGHT_)
4263                res = (size_t)trunc(fabs(dp_T0[arg]));
4264                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4265                dp_T0[res] *= coval;
4266#if !defined(_INDO_)
4267#if !defined(_ZOS_) /* BREAK_ZOS */
4268#if !defined( _INT_FOR_)
4269
4270                FOR_0_LE_l_LT_pk
4271                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
4272
4273                FOR_0_LE_l_LT_pk
4274                TRES_INC *= coval;
4275#endif
4276#endif
4277#endif
4278#else
4279                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4280                exit(-2);
4281#endif /* ALL_TOGETHER_AGAIN */
4282                break;
4283
4284            case ref_eq_mult_a:       /* Multiply one adouble by another    eq_mult_a */
4285                /* (*=) */
4286                arg = get_locint_f();
4287                arg1 = get_locint_f();
4288
4289#if !defined(_NTIGHT_)
4290                res = (size_t)trunc(fabs(dp_T0[arg1]));
4291                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4292#if defined(_INDO_)
4293                //merge_2_index_domains(res, arg, ind_dom);
4294#if defined(_NONLIND_)
4295                //extend_nonlinearity_domain_binary(res, arg, ind_dom, nonl_dom);
4296#endif
4297#else
4298#if !defined(_ZOS_) /* BREAK_ZOS */
4299                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
4300                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
4301
4302                INC_pk_1(Tres)
4303                INC_pk_1(Targ)
4304
4305#ifdef _INT_FOR_
4306                FOR_p_GT_l_GE_0
4307                TRES_FODEC |= TARG_DEC;
4308#else
4309                FOR_p_GT_l_GE_0
4310                FOR_k_GT_i_GE_0
4311                { TRES_FODEC = dp_T0[res]*TARG_DEC +
4312                               TRES*dp_T0[arg];
4313                  DEC_TRES_FO
4314#ifdef _HIGHER_ORDER_
4315                  TresOP = Tres-i;
4316                  TargOP = Targ;
4317
4318                  for (j=0;j<i;j++)
4319                  *Tres += (*TresOP++) * (*TargOP--);
4320                  Tres--;
4321#endif /* _HIGHER_ORDER_ */
4322                }
4323#endif
4324#endif
4325#endif
4326                dp_T0[res] *= dp_T0[arg];
4327#else
4328                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4329                exit(-2);
4330#endif /* !_NTIGHT_ */
4331                break;
4332
4333            case ref_cond_assign:                                      /* cond_assign */
4334                arg   = get_locint_f();
4335                arg1  = get_locint_f();
4336                arg2  = get_locint_f();
4337                { 
4338                    locint ref = get_locint_f();
4339                    coval = get_val_f();
4340#if !defined(_NTIGHT_)
4341                    res   = (size_t)trunc(fabs(dp_T0[ref]));
4342
4343                    IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4344
4345                /* olvo 980924 changed order to allow reflexive ops */
4346#if defined(_INDO_)
4347                if (dp_T0[arg] > 0) {
4348                    if (coval <= 0.0)
4349                        MINDEC(ret_c,2);
4350                    dp_T0[res] = dp_T0[arg1];
4351
4352                    //combine_2_index_domains(res, arg1, arg2, ind_dom);
4353                } else {
4354                    if (coval > 0.0)
4355                        MINDEC(ret_c,2);
4356                    if (dp_T0[arg] == 0)
4357                        MINDEC(ret_c,0);
4358                    dp_T0[res] = dp_T0[arg2];
4359
4360                    //combine_2_index_domains(res, arg1, arg2, ind_dom);
4361                }
4362#else
4363#if !defined(_ZOS_) /* BREAK_ZOS */
4364                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
4365                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
4366                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
4367#endif /* ALL_TOGETHER_AGAIN */
4368
4369#ifdef _INT_FOR_
4370                coval = get_val_f();
4371
4372                if (dp_T0[arg] > 0)
4373                    FOR_0_LE_l_LT_pk
4374                    TRES_INC = TARG1_INC;
4375                else
4376                    FOR_0_LE_l_LT_pk
4377                    TRES_INC = TARG2_INC;
4378
4379                if (dp_T0[arg] > 0) {
4380                    if (coval <= 0.0)
4381                        MINDEC(ret_c,2);
4382                    dp_T0[res] = dp_T0[arg1];
4383                } else {
4384                    if (coval > 0.0)
4385                        MINDEC(ret_c,2);
4386                    if (dp_T0[arg] == 0)
4387                        MINDEC(ret_c,0);
4388                    dp_T0[res] = dp_T0[arg2];
4389                }
4390                FOR_0_LE_l_LT_pk
4391                TRES_INC = TARG1_INC | TARG2_INC;
4392#else
4393#if !defined(_ZOS_) /* BREAK_ZOS */
4394                if (dp_T0[arg] > 0)
4395                    FOR_0_LE_l_LT_pk
4396                    TRES_INC = TARG1_INC;
4397                else
4398                    FOR_0_LE_l_LT_pk
4399                    TRES_INC = TARG2_INC;
4400#endif
4401
4402                if (dp_T0[arg] > 0) {
4403                    if (coval <= 0.0)
4404                        MINDEC(ret_c,2);
4405                    dp_T0[res] = dp_T0[arg1];
4406                } else {
4407                    if (coval > 0.0)
4408                        MINDEC(ret_c,2);
4409                    if (dp_T0[arg] == 0)
4410                        MINDEC(ret_c,0);
4411                    dp_T0[res] = dp_T0[arg2];
4412                }
4413#endif
4414#endif
4415#else
4416                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4417                exit(-2);
4418#endif /* ALL_TOGETHER_AGAIN */
4419                }
4420                break;
4421
4422            case ref_cond_assign_s:                                  /* cond_assign_s */
4423                arg   = get_locint_f();
4424                arg1  = get_locint_f();
4425                arg2   = get_locint_f();
4426                coval = get_val_f();
4427
4428#if !defined(_NTIGHT_)
4429                res = (size_t)trunc(fabs(dp_T0[arg2]));
4430                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4431
4432                /* olvo 980924 changed order to allow reflexive ops */
4433#if defined(_INDO_)
4434                    //copy_index_domain(res, arg1, ind_dom);
4435#else
4436#if !defined(_ZOS_) /* BREAK_ZOS */
4437                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
4438                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
4439#endif /* ALL_TOGETHER_AGAIN */
4440
4441#ifdef _INT_FOR_
4442                coval = get_val_f();
4443
4444                if (dp_T0[arg] > 0)
4445                    FOR_0_LE_l_LT_pk
4446                    TRES_INC = TARG1_INC;
4447
4448                if (dp_T0[arg] > 0) {
4449                    if (coval <= 0.0)
4450                        MINDEC(ret_c,2);
4451                    dp_T0[res] = dp_T0[arg1];
4452                } else
4453                    if (dp_T0[arg] == 0)
4454                        MINDEC(ret_c,0);
4455#else
4456#if !defined(_ZOS_) /* BREAK_ZOS */
4457                if (dp_T0[arg] > 0)
4458                    FOR_0_LE_l_LT_pk
4459                    TRES_INC = TARG1_INC;
4460#endif
4461                if (dp_T0[arg] > 0) {
4462                    if (coval <= 0.0)
4463                        MINDEC(ret_c,2);
4464                    dp_T0[res] = dp_T0[arg1];
4465                } else
4466                    if (dp_T0[arg] == 0)
4467                        MINDEC(ret_c,0);
4468#endif
4469#endif
4470#else
4471                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4472                exit(-2);
4473#endif /* ALL_TOGETHER_AGAIN */
4474                break;
4475
4476                /****************************************************************************/
4477                /*                                                          REMAINING STUFF */
4478
4479                /*--------------------------------------------------------------------------*/
4480            case take_stock_op:                                  /* take_stock_op */
4481                size = get_locint_f();
4482                res  = get_locint_f();
4483                d    = get_val_v_f(size);
4484
4485                for (ls=0;ls<size;ls++) {
4486#if !defined(_NTIGHT_)
4487                    dp_T0[res]=*d;
4488#endif /* !_NTIGHT_ */
4489#if !defined(_INDO_)
4490#if !defined(_ZOS_) /* BREAK_ZOS */
4491                    ASSIGN_T(Tres,TAYLOR_BUFFER[res])
4492
4493                    FOR_0_LE_l_LT_pk
4494                    TRES_INC = 0;
4495
4496#endif /* ALL_TOGETHER_AGAIN */
4497                    res++;
4498#if !defined(_NTIGHT_)
4499                    d++;
4500#endif /* !_NTIGHT_ */
4501#endif
4502                }
4503                break;
4504
4505                /*--------------------------------------------------------------------------*/
4506            case death_not:                                          /* death_not */
4507                arg1=get_locint_f();
4508                arg2=get_locint_f();
4509
4510#ifdef _KEEP_
4511                if (keep) {
4512                    do {
4513                        IF_KEEP_WRITE_TAYLOR(arg2,keep,k,p)
4514                    } while(arg1 < arg2-- );
4515                }
4516#endif
4517                break;
4518
4519                /*--------------------------------------------------------------------------*/
4520#if defined(_EXTERN_) /* ZOS,  FOS, FOV up to now */
4521            case ext_diff:                       /* extern differntiated function */
4522                ADOLC_CURRENT_TAPE_INFOS.ext_diff_fct_index=get_locint_f();
4523                n=get_locint_f();
4524                m=get_locint_f();
4525                ADOLC_CURRENT_TAPE_INFOS.lowestXLoc_for = get_locint_f();
4526                ADOLC_CURRENT_TAPE_INFOS.lowestYLoc_for = get_locint_f();
4527                ADOLC_CURRENT_TAPE_INFOS.cpIndex = get_locint_f();
4528                edfct=get_ext_diff_fct(ADOLC_CURRENT_TAPE_INFOS.ext_diff_fct_index);
4529
4530                if (edfct->ADOLC_EXT_FCT_POINTER==NULL)
4531                    fail(ADOLC_EXT_DIFF_NULLPOINTER_DIFFFUNC);
4532                if (n>0) {
4533                    if (edfct->dp_x==NULL) fail(ADOLC_EXT_DIFF_NULLPOINTER_ARGUMENT);
4534#if !defined(_ZOS_)
4535                    if (ADOLC_EXT_POINTER_X==NULL) fail(ADOLC_EXT_DIFF_NULLPOINTER_ARGUMENT);
4536#endif
4537                }
4538                if (m>0) {
4539                    if (edfct->dp_y==NULL) fail(ADOLC_EXT_DIFF_NULLPOINTER_ARGUMENT);
4540#if !defined(_ZOS_)
4541                    if (ADOLC_EXT_POINTER_Y==NULL) fail(ADOLC_EXT_DIFF_NULLPOINTER_ARGUMENT);
4542#endif
4543                }
4544
4545                arg = ADOLC_CURRENT_TAPE_INFOS.lowestXLoc_for;
4546                for (loop=0; loop<n; ++loop) {
4547                    edfct->dp_x[loop]=dp_T0[arg];
4548#if !defined(_ZOS_)
4549                    ADOLC_EXT_LOOP
4550                        ADOLC_EXT_POINTER_X[loop]ADOLC_EXT_SUBSCRIPT =
4551                        TAYLOR_BUFFER[arg]ADOLC_EXT_SUBSCRIPT;
4552#endif
4553                    ++arg;
4554                }
4555                arg = ADOLC_CURRENT_TAPE_INFOS.lowestYLoc_for;
4556                for (loop=0; loop<m; ++loop) {
4557                    edfct->dp_y[loop]=dp_T0[arg];
4558#if !defined(_ZOS_)
4559                    ADOLC_EXT_LOOP
4560                        ADOLC_EXT_POINTER_Y[loop]ADOLC_EXT_SUBSCRIPT =
4561                        TAYLOR_BUFFER[arg]ADOLC_EXT_SUBSCRIPT;
4562#endif
4563                    ++arg;
4564                }
4565
4566                ext_retc = edfct->ADOLC_EXT_FCT_COMPLETE;
4567                MINDEC(ret_c, ext_retc);
4568
4569                res = ADOLC_CURRENT_TAPE_INFOS.lowestXLoc_for;
4570                for (loop=0; loop<n; ++loop) {
4571                    IF_KEEP_WRITE_TAYLOR(res, keep, k, p);
4572                    dp_T0[res]=edfct->dp_x[loop];
4573#if !defined(_ZOS_)
4574                    ADOLC_EXT_LOOP
4575                        TAYLOR_BUFFER[res]ADOLC_EXT_SUBSCRIPT =
4576                        ADOLC_EXT_POINTER_X[loop]ADOLC_EXT_SUBSCRIPT;
4577#endif
4578                    ++res;
4579                }
4580                res = ADOLC_CURRENT_TAPE_INFOS.lowestYLoc_for;
4581                for (loop=0; loop<m; ++loop) {
4582                    IF_KEEP_WRITE_TAYLOR(res, keep, k, p);
4583                    dp_T0[res]=edfct->dp_y[loop];
4584#if !defined(_ZOS_)
4585                    ADOLC_EXT_LOOP
4586                        TAYLOR_BUFFER[res]ADOLC_EXT_SUBSCRIPT =
4587                        ADOLC_EXT_POINTER_Y[loop]ADOLC_EXT_SUBSCRIPT;
4588#endif
4589                    ++res;
4590                }
4591
4592                break;
4593#endif
4594                /*--------------------------------------------------------------------------*/
4595
4596            default:                                                   /* default */
4597                /* Die here, we screwed up */
4598
4599                fprintf(DIAG_OUT,"ADOL-C fatal error in " GENERATED_FILENAME " ("
4600                        __FILE__
4601                        ") : no such operation %d\n", operation);
4602                exit(-1);
4603                break;
4604
4605        } /* endswitch */
4606
4607        /* Read the next operation */
4608        operation=get_op_f();
4609#if defined(ADOLC_DEBUG)
4610        ++countPerOperation[operation];
4611#endif /* ADOLC_DEBUG */
4612    }  /* endwhile */
4613
4614
4615#if defined(ADOLC_DEBUG)
4616    printf("\nTape contains:\n");
4617    for (v = 0; v < 256; ++v)
4618        if (countPerOperation[v] > 0)
4619            printf("operation %3d: %6d time(s) - %6d taylors written (%10.2f per operation)\n", v, countPerOperation[v], taylorPerOperation[v], (double)taylorPerOperation[v] / (double)countPerOperation[v]);
4620    printf("\n");
4621#endif /* ADOLC_DEBUG */
4622
4623#if defined(_KEEP_)
4624    if (keep) taylor_close(taylbuf);
4625#endif
4626
4627    /* clean up */
4628#if !defined (_NTIGHT_)
4629    free(dp_T0);
4630#endif /* !_NTIGHT_ */
4631#if !defined(_INDO_)
4632#if !defined(_ZOS_)
4633#   if defined(_FOS_)
4634    free(dp_T);
4635#   else
4636#if !defined (_INT_FOR_)
4637    myfree2(dpp_T);
4638    free(dp_Ttemp);
4639#endif /* !_NTIGHT_ */
4640#endif
4641#endif
4642#endif
4643#if defined(_HIGHER_ORDER_)
4644    free(dp_z);
4645#endif
4646
4647    end_sweep();
4648
4649
4650#if defined(_INDO_)
4651#if defined(_INDOPRO_)
4652    for(i=0;i<max_ind_dom;i++)
4653      {
4654        free(ind_dom[i]);
4655      }
4656    free(ind_dom);
4657#endif
4658#if defined(_NONLIND_)
4659    for( i=0; i < indcheck; i++) {
4660      traverse_crs(&nonl_dom[i],&sod[i],indcheck+1);
4661      free_tree(&nonl_dom[i],indcheck+1);
4662      crs[i] = (unsigned int*) malloc(sizeof(unsigned int) * (sod[i].entry+1));
4663      crs[i][0] = sod[i].entry;
4664      temp = sod[i].left;
4665      for( ii=1; ii <=sod[i].entry; ii++)
4666        {
4667          crs[i][ii] = temp->entry;
4668          temp1 = temp->left;
4669          free(temp);
4670          temp = temp1;
4671        }
4672    }
4673
4674    free(sod);
4675    free(nonl_dom);
4676    free(fod);
4677    free(arg_index);
4678
4679#endif
4680#endif
4681    return ret_c;
4682}
4683
4684/****************************************************************************/
4685
4686#if defined(_INDOPRO_)
4687
4688/****************************************************************************/
4689/* set operations for propagation of index domains                          */
4690
4691/*--------------------------------------------------------------------------*/
4692/* operations on index domains                                              */
4693
4694#if defined(_TIGHT_)
4695void copy_index_domain(int res, int arg, locint **ind_dom) {
4696
4697   int i;
4698
4699   if (ind_dom[arg][0] > ind_dom[res][1])
4700     {
4701       free(ind_dom[res]);
4702       ind_dom[res] = (locint *)  malloc(sizeof(locint) * 2*(ind_dom[arg][0]+1));
4703       ind_dom[res][1] = 2*ind_dom[arg][0];
4704     }
4705
4706
4707    for(i=2;i<ind_dom[arg][0]+2;i++)
4708       ind_dom[res][i] = ind_dom[arg][i];
4709    ind_dom[res][0] = ind_dom[arg][0];
4710}
4711
4712
4713void merge_2_index_domains(int res, int arg, locint **ind_dom)
4714{
4715
4716  int num,num1,num2, i,j,k,l;
4717  locint *temp_array, *arg_ind_dom, *res_ind_dom;
4718
4719  if (ind_dom[res][0] == 0)
4720    copy_index_domain(res,arg,ind_dom);
4721  else
4722    {
4723      if (res != arg)
4724     {
4725       arg_ind_dom = ind_dom[arg];
4726       res_ind_dom = ind_dom[res];
4727
4728       num  = ind_dom[res][0];
4729       num1 = arg_ind_dom[0];
4730       num2 = ind_dom[res][1];
4731
4732       if (num2 < num1+num)
4733         num2 = num1+num;
4734
4735       temp_array = (locint *)  malloc(sizeof(locint)* (num2+2));
4736       temp_array[1] = num2;
4737
4738       i = 2;
4739       j = 2;
4740       k = 2;
4741       num += 2;
4742       num1 += 2;
4743       while ((i< num) && (j < num1))
4744         {
4745           if (res_ind_dom[i] < arg_ind_dom[j])
4746          {
4747            temp_array[k] = res_ind_dom[i];
4748            i++; k++;
4749          }
4750           else
4751          {
4752            if (res_ind_dom[i] == arg_ind_dom[j])
4753              {
4754                temp_array[k] = arg_ind_dom[j];
4755                i++;j++;k++;
4756              }
4757            else
4758              {
4759                temp_array[k] = arg_ind_dom[j];
4760                j++;k++;
4761              }
4762          }
4763         }
4764       for(l = i;l<num;l++)
4765         {
4766           temp_array[k] = res_ind_dom[l];
4767           k++;
4768         }
4769       for(l = j;l<num1;l++)
4770         {
4771           temp_array[k] = arg_ind_dom[l];
4772           k++;
4773         }
4774       temp_array[0] = k-2;
4775       free(ind_dom[res]);
4776       ind_dom[res]=temp_array;
4777     }
4778    }
4779
4780
4781}
4782
4783void combine_2_index_domains(int res, int arg1, int arg2, locint **ind_dom) {
4784
4785    if (res != arg1)
4786       copy_index_domain(res, arg1, ind_dom);
4787
4788    merge_2_index_domains(res, arg2, ind_dom);
4789}
4790
4791void merge_3_index_domains(int res, int arg1, int arg2, locint **ind_dom) {
4792    merge_2_index_domains(res, arg1, ind_dom);
4793    merge_2_index_domains(res, arg2, ind_dom);
4794}
4795
4796
4797
4798#endif
4799#endif
4800
4801
4802#if defined(_NONLIND_)
4803#if defined(_TIGHT_)
4804
4805void free_tree(IndexElement* tree, int num)
4806{
4807
4808  if (tree->left != NULL)
4809    {
4810      free_tree(tree->left,num);
4811    }
4812  if (tree->right != NULL)
4813    {
4814      free_tree(tree->right,num);
4815     }
4816    {
4817      if (tree->entry == num)
4818        free(tree);
4819
4820    }
4821 
4822}
4823void traverse_crs(IndexElement* tree,  IndexElement_sod* sod, int num)
4824{
4825
4826  IndexElement_sod *temp, *temp1;
4827  int ii;
4828
4829  if (tree->left != NULL)
4830    {
4831      traverse_crs(tree->left, sod, num);
4832    }
4833  if (tree->right != NULL)
4834    {
4835      traverse_crs(tree->right, sod, num);
4836    }
4837  if (tree->entry < num)
4838    {
4839      temp = sod->left;
4840      if (temp == NULL)
4841        {
4842          temp = (struct IndexElement_sod*) malloc(sizeof(struct IndexElement_sod));
4843          temp->left = NULL;
4844          temp->entry = tree->entry;
4845          sod->entry++;
4846          sod->left=temp;
4847        }
4848      else
4849        {
4850          while ((temp->entry < tree->entry) && (temp->left != NULL))
4851            {
4852              temp1 = temp;
4853              temp = temp->left;
4854            }
4855          if (temp->left == NULL)
4856            {
4857              if(temp->entry < tree->entry)
4858                {
4859                  temp->left = (struct IndexElement_sod*) malloc(sizeof(struct IndexElement_sod));
4860                  temp = temp->left;
4861                  temp->left = NULL;
4862                  temp->entry = tree->entry;
4863                  sod->entry++;
4864                }
4865              if(temp->entry > tree->entry)
4866                {
4867                  temp->left = (struct IndexElement_sod*) malloc(sizeof(struct IndexElement_sod));
4868                  temp->left->entry = temp->entry;
4869                  temp->left->left = NULL;
4870                  temp->entry = tree->entry;
4871                  sod->entry++;
4872                }
4873            }
4874          else
4875            {
4876              if (temp->entry > tree->entry)
4877                {
4878                  temp1 = (struct IndexElement_sod*) malloc(sizeof(struct IndexElement_sod));
4879                  temp1->left = temp->left;
4880                  temp1->entry = temp->entry;
4881                  temp->entry = tree->entry;
4882                  temp->left=temp1;
4883                  sod->entry++;
4884                }
4885             
4886            }
4887        }
4888    }
4889}
4890
4891void traverse_unary(IndexElement* tree,  IndexElement* nonl_dom,  IndexElement* fodi, int num, int maxopind)
4892{
4893  IndexElement *temp;
4894
4895  if (tree->left != NULL)
4896    {
4897      traverse_unary(tree->left, nonl_dom, fodi, num, maxopind);
4898      if (tree->right != NULL)
4899        {
4900          traverse_unary(tree->right, nonl_dom, fodi, num, maxopind);
4901        }
4902     }
4903  else
4904    {
4905      if(tree->entry<maxopind)
4906        {
4907          temp = (struct IndexElement*) malloc(sizeof(struct IndexElement));
4908          temp->right = fodi;
4909          temp->left = nonl_dom[tree->entry].left;
4910          temp->entry= num;
4911          nonl_dom[tree->entry].left = temp;
4912        }
4913    }
4914}
4915
4916#endif
4917#endif
4918END_C_DECLS
Note: See TracBrowser for help on using the repository browser.