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

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

reactivate index domain propagation for active subscripts and references

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

  • Property svn:keywords set to Author Date Id Revision
File size: 155.1 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 362 2012-11-14 15:14:21Z 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#if defined(_INDOPRO_)
3947                    copy_index_domain(res, arg1, ind_dom);
3948#endif
3949#if defined(_NONLIND_)
3950                    arg_index[res] = arg_index[arg1];
3951#endif
3952#else
3953#if !defined(_ZOS_) /* BREAK_ZOS */
3954                    ASSIGN_T(Targ1,TAYLOR_BUFFER[arg1])
3955                    ASSIGN_T(Tres,TAYLOR_BUFFER[res])
3956
3957                    FOR_0_LE_l_LT_pk
3958                    TRES_INC = TARG1_INC;
3959#endif
3960#endif
3961#else
3962                    fprintf(DIAG_OUT, "ADOL-C error: active subscripting does not work in safe mode, please use tight mode\n");
3963                    exit(-2);
3964#endif /* ALL_TOGETHER_AGAIN */
3965                }
3966                break;
3967
3968            case subscript_ref:
3969                coval = get_val_f();
3970                arg = get_locint_f();
3971                res = get_locint_f();
3972                {
3973                    size_t cnt, idx, numvar = (size_t)trunc(fabs(coval));
3974                    locint vectorloc;
3975                    vectorloc = get_locint_f();
3976#if !defined(_NTIGHT_)
3977                    idx = (size_t)trunc(fabs(dp_T0[arg]));
3978                    if (idx >= numvar)
3979                        fprintf(DIAG_OUT, "ADOL-C warning: index out of bounds while subscripting (ref) n=%z, idx=%z\n", numvar, idx);
3980                    arg1 = vectorloc+idx;
3981                    IF_KEEP_WRITE_TAYLOR(res,keep,k,p);
3982                    dp_T0[res] = arg1;
3983#else
3984                    fprintf(DIAG_OUT, "ADOL-C error: active subscripting does not work in safe mode, please use tight mode\n");
3985                    exit(-2);
3986#endif
3987                }
3988                break;
3989
3990            case ref_copyout:
3991                arg = get_locint_f();
3992                res = get_locint_f();
3993#if !defined(_NTIGHT_)
3994                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
3995                IF_KEEP_WRITE_TAYLOR(res,keep,k,p);
3996                dp_T0[res] = dp_T0[arg1];
3997#if defined(_INDO_)
3998#if defined(_INDOPRO_)
3999                copy_index_domain(res, arg1, ind_dom);
4000#endif
4001#if defined(_NONLIND_)
4002                arg_index[res] = arg_index[arg1];
4003#endif
4004#else
4005#if !defined(_ZOS_) /* BREAK_ZOS */
4006                ASSIGN_T(Targ1,TAYLOR_BUFFER[arg1])
4007                ASSIGN_T(Tres,TAYLOR_BUFFER[res])
4008
4009                FOR_0_LE_l_LT_pk
4010                TRES_INC = TARG1_INC;
4011#endif
4012#endif
4013#else
4014                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4015                exit(-2);
4016#endif /* ALL_TOGETHER_AGAIN */
4017                break;
4018
4019            case ref_incr_a:
4020                arg = get_locint_f();
4021#if !defined(_NTIGHT_)
4022                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
4023                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p);
4024                dp_T0[arg1]++;
4025#else
4026                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4027                exit(-2);
4028#endif
4029                break;
4030
4031            case ref_decr_a:
4032                arg = get_locint_f();
4033#if !defined(_NTIGHT_)
4034                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
4035                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p);
4036                dp_T0[arg1]--;
4037#else
4038                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4039                exit(-2);
4040#endif
4041                break;
4042
4043            case ref_assign_d:
4044                arg = get_locint_f();
4045                coval = get_val_f();
4046               
4047#if !defined(_NTIGHT_)
4048                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
4049                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p)
4050                dp_T0[arg1] = coval;
4051#if defined(_INDO_)
4052#if defined(_INDOPRO_)
4053                ind_dom[arg1][0] = 0;
4054#endif
4055#if defined(_NONLIND_)
4056                fod[opind].entry = maxopind+2;
4057                fod[opind].left = NULL;
4058                fod[opind].right = NULL;
4059                arg_index[arg1] = opind++;
4060#endif
4061#else
4062#if !defined(_ZOS_)
4063                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
4064
4065                FOR_0_LE_l_LT_pk
4066                TARG1_INC = 0;
4067#endif
4068#endif
4069#else
4070                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4071                exit(-2);
4072#endif
4073                break;
4074
4075            case ref_assign_d_zero:
4076                arg = get_locint_f();
4077
4078#if !defined(_NTIGHT_)
4079                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
4080                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p)
4081                dp_T0[arg1] = 0.0;
4082#if defined(_INDO_)
4083#if defined(_INDOPRO_)
4084                ind_dom[arg1][0] = 0;
4085#endif
4086#if defined(_NONLIND_)
4087                fod[opind].entry = maxopind+2;
4088                fod[opind].left = NULL;
4089                fod[opind].right = NULL;
4090                arg_index[arg1] = opind++;
4091#endif
4092#else
4093#if !defined(_ZOS_)
4094                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
4095
4096                FOR_0_LE_l_LT_pk
4097                TARG1_INC = 0;
4098#endif
4099#endif
4100#else
4101                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4102                exit(-2);
4103#endif
4104                break;
4105
4106            case ref_assign_d_one:
4107                arg = get_locint_f();
4108
4109#if !defined(_NTIGHT_)
4110                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
4111                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p)
4112                dp_T0[arg1] = 1.0;
4113#if defined(_INDO_)
4114#if defined(_INDOPRO_)
4115                ind_dom[arg1][0] = 0;
4116#endif
4117#if defined(_NONLIND_)
4118                fod[opind].entry = maxopind+2;
4119                fod[opind].left = NULL;
4120                fod[opind].right = NULL;
4121                arg_index[arg1] = opind++;
4122#endif
4123#else
4124#if !defined(_ZOS_)
4125                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
4126
4127                FOR_0_LE_l_LT_pk
4128                TARG1_INC = 0;
4129#endif
4130#endif
4131#else
4132                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4133                exit(-2);
4134#endif
4135                break;
4136
4137            case ref_assign_a:           /* assign an adouble variable an    assign_a */
4138                /* adouble value. (=) */
4139                arg = get_locint_f();
4140                res = get_locint_f();
4141
4142#if !defined(_NTIGHT_)
4143                arg1 = (size_t)trunc(fabs(dp_T0[res]));
4144                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p)
4145                dp_T0[arg1] = dp_T0[arg];
4146#if defined(_INDO_)
4147#if defined(_INDOPRO_)
4148                copy_index_domain(arg1, arg, ind_dom);
4149#endif
4150#if defined(_NONLIND_)
4151                arg_index[arg1] = arg_index[arg];
4152#endif
4153#else
4154#if !defined(_ZOS_) /* BREAK_ZOS */
4155                ASSIGN_T(Targ,TAYLOR_BUFFER[arg])
4156                ASSIGN_T(Targ1,TAYLOR_BUFFER[arg1])
4157
4158                FOR_0_LE_l_LT_pk
4159                TARG1_INC = TARG_INC;
4160#endif
4161#endif
4162#else
4163                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4164                exit(-2);
4165#endif /* ALL_TOGETHER_AGAIN */
4166                break;
4167
4168            case ref_assign_ind:       /* assign an adouble variable an    assign_ind */
4169                /* independent double value (<<=) */
4170                arg = get_locint_f();
4171
4172
4173#if !defined(_NTIGHT_)
4174                res = (size_t)trunc(fabs(dp_T0[arg]));
4175                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4176                dp_T0[res] = basepoint[indexi];
4177#if defined(_INDO_)
4178#if defined(_INDOPRO_)
4179                ind_dom[res][0] = 1;
4180                ind_dom[res][2] = indexi;
4181#endif
4182#if defined(_NONLIND_)
4183                fod[opind].entry = indexi;
4184                fod[opind].left = NULL;
4185                fod[opind].right = NULL;
4186                arg_index[res] = opind++;
4187#endif
4188#else
4189#if !defined(_ZOS_) /* BREAK_ZOS */
4190                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
4191
4192#ifdef _INT_FOR_
4193                FOR_0_LE_l_LT_p
4194                TRES_INC = ARGUMENT(indexi,l,i);
4195#else
4196                FOR_0_LE_l_LT_p
4197                FOR_0_LE_i_LT_k
4198                TRES_INC = ARGUMENT(indexi,l,i);
4199#endif
4200#endif
4201#endif
4202#else
4203                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4204                exit(-2);
4205#endif /* ALL_TOGETHER_AGAIN */
4206                ++indexi;
4207                break;
4208
4209            case ref_eq_plus_d:            /* Add a floating point to an    eq_plus_d */
4210                /* adouble. (+=) */
4211                arg  = get_locint_f();
4212                coval = get_val_f();
4213
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_plus_a:             /* Add an adouble to another    eq_plus_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#if defined(_INDOPRO_)
4237                merge_2_index_domains(res, arg, ind_dom);
4238#endif
4239#if defined(_NONLIND_)
4240                fod[opind].entry = maxopind+2;
4241                fod[opind].left = &fod[arg_index[res]];
4242                fod[opind].right = &fod[arg_index[arg]];
4243                arg_index[res] = opind++;
4244#endif
4245#else
4246#if !defined(_ZOS_) /* BREAK_ZOS */
4247                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
4248                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
4249
4250#ifdef _INT_FOR_
4251                FOR_0_LE_l_LT_pk
4252                TRES_INC |= TARG_INC;
4253#else
4254                FOR_0_LE_l_LT_pk
4255                TRES_INC += TARG_INC;
4256#endif
4257#endif
4258#endif
4259#else
4260                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4261                exit(-2);
4262#endif /* ALL_TOGETHER_AGAIN */
4263                break;
4264
4265            case ref_eq_min_d:       /* Subtract a floating point from an    eq_min_d */
4266                /* adouble. (-=) */
4267                arg = get_locint_f();
4268                coval = get_val_f();
4269
4270#if !defined(_NTIGHT_)
4271                res = (size_t)trunc(fabs(dp_T0[arg]));
4272                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4273                dp_T0[res] -= coval;
4274#else
4275                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4276                exit(-2);
4277#endif /* !_NTIGHT_ */
4278                break;
4279
4280                /*--------------------------------------------------------------------------*/
4281            case ref_eq_min_a:        /* Subtract an adouble from another    eq_min_a */
4282                /* adouble. (-=) */
4283                arg = get_locint_f();
4284                arg1 = get_locint_f();
4285
4286#if !defined(_NTIGHT_)
4287                res = (size_t)trunc(fabs(dp_T0[arg1]));
4288                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4289                dp_T0[res] -= dp_T0[arg];
4290#if defined(_INDO_)
4291#if defined(_INDOPRO_)
4292                merge_2_index_domains(res, arg, ind_dom);
4293#endif
4294#if defined(_NONLIND_)
4295                fod[opind].entry = maxopind+2;
4296                fod[opind].left = &fod[arg_index[res]];
4297                fod[opind].right = &fod[arg_index[arg]];
4298                arg_index[res] = opind++;
4299#endif
4300#else
4301#if !defined(_ZOS_) /* BREAK_ZOS */
4302                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
4303                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
4304
4305#ifdef _INT_FOR_
4306                FOR_0_LE_l_LT_pk
4307                TRES_INC |= TARG_INC;
4308#else
4309                FOR_0_LE_l_LT_pk
4310                TRES_INC -= TARG_INC;
4311#endif
4312#endif
4313#endif
4314#else
4315                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4316                exit(-2);
4317#endif /* ALL_TOGETHER_AGAIN */
4318                break;
4319
4320            case ref_eq_mult_d:              /* Multiply an adouble by a    eq_mult_d */
4321                /* flaoting point. (*=) */
4322                arg = get_locint_f();
4323                coval = get_val_f();
4324
4325#if !defined(_NTIGHT_)
4326                res = (size_t)trunc(fabs(dp_T0[arg]));
4327                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4328                dp_T0[res] *= coval;
4329#if !defined(_INDO_)
4330#if !defined(_ZOS_) /* BREAK_ZOS */
4331#if !defined( _INT_FOR_)
4332
4333                FOR_0_LE_l_LT_pk
4334                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
4335
4336                FOR_0_LE_l_LT_pk
4337                TRES_INC *= coval;
4338#endif
4339#endif
4340#endif
4341#else
4342                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4343                exit(-2);
4344#endif /* ALL_TOGETHER_AGAIN */
4345                break;
4346
4347            case ref_eq_mult_a:       /* Multiply one adouble by another    eq_mult_a */
4348                /* (*=) */
4349                arg = get_locint_f();
4350                arg1 = get_locint_f();
4351
4352#if !defined(_NTIGHT_)
4353                res = (size_t)trunc(fabs(dp_T0[arg1]));
4354                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4355#if defined(_INDO_)
4356#if defined(_INDOPRO_)
4357                merge_2_index_domains(res, arg, ind_dom);
4358#endif
4359#if defined(_NONLIND_)
4360                fod[opind].entry = maxopind+2;
4361                fod[opind].left = &fod[arg_index[res]];
4362                fod[opind].right = &fod[arg_index[arg]];
4363                traverse_unary(&fod[arg_index[res]], nonl_dom, &fod[arg_index[arg]], indcheck+1,maxopind+2);
4364                traverse_unary(&fod[arg_index[arg]], nonl_dom, &fod[arg_index[res]], indcheck+1,maxopind+2);
4365                arg_index[res] = opind++;
4366#endif
4367#else
4368#if !defined(_ZOS_) /* BREAK_ZOS */
4369                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
4370                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
4371
4372                INC_pk_1(Tres)
4373                INC_pk_1(Targ)
4374
4375#ifdef _INT_FOR_
4376                FOR_p_GT_l_GE_0
4377                TRES_FODEC |= TARG_DEC;
4378#else
4379                FOR_p_GT_l_GE_0
4380                FOR_k_GT_i_GE_0
4381                { TRES_FODEC = dp_T0[res]*TARG_DEC +
4382                               TRES*dp_T0[arg];
4383                  DEC_TRES_FO
4384#ifdef _HIGHER_ORDER_
4385                  TresOP = Tres-i;
4386                  TargOP = Targ;
4387
4388                  for (j=0;j<i;j++)
4389                  *Tres += (*TresOP++) * (*TargOP--);
4390                  Tres--;
4391#endif /* _HIGHER_ORDER_ */
4392                }
4393#endif
4394#endif
4395#endif
4396                dp_T0[res] *= dp_T0[arg];
4397#else
4398                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4399                exit(-2);
4400#endif /* !_NTIGHT_ */
4401                break;
4402
4403            case ref_cond_assign:                                      /* cond_assign */
4404                arg   = get_locint_f();
4405                arg1  = get_locint_f();
4406                arg2  = get_locint_f();
4407                { 
4408                    locint ref = get_locint_f();
4409                    coval = get_val_f();
4410#if !defined(_NTIGHT_)
4411                    res   = (size_t)trunc(fabs(dp_T0[ref]));
4412
4413                    IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4414
4415                /* olvo 980924 changed order to allow reflexive ops */
4416#if defined(_INDO_)
4417                if (dp_T0[arg] > 0) {
4418                    if (coval <= 0.0)
4419                        MINDEC(ret_c,2);
4420                    dp_T0[res] = dp_T0[arg1];
4421
4422#if defined(_INDOPRO_)
4423                    copy_index_domains(res, arg1, ind_dom);
4424#endif
4425#if defined(_NONLIND_)
4426                    arg_index[res] = arg_index[arg1];
4427#endif
4428                } else {
4429                    if (coval > 0.0)
4430                        MINDEC(ret_c,2);
4431                    if (dp_T0[arg] == 0)
4432                        MINDEC(ret_c,0);
4433                    dp_T0[res] = dp_T0[arg2];
4434
4435#if defined(_INDOPRO_)
4436                    copy_index_domains(res, arg2, ind_dom);
4437#endif
4438#if defined(_NONLIND_)
4439                    arg_index[res] = arg_index[arg2];
4440#endif
4441                }
4442#else
4443#if !defined(_ZOS_) /* BREAK_ZOS */
4444                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
4445                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
4446                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
4447#endif /* ALL_TOGETHER_AGAIN */
4448
4449#ifdef _INT_FOR_
4450                coval = get_val_f();
4451
4452                if (dp_T0[arg] > 0)
4453                    FOR_0_LE_l_LT_pk
4454                    TRES_INC = TARG1_INC;
4455                else
4456                    FOR_0_LE_l_LT_pk
4457                    TRES_INC = TARG2_INC;
4458
4459                if (dp_T0[arg] > 0) {
4460                    if (coval <= 0.0)
4461                        MINDEC(ret_c,2);
4462                    dp_T0[res] = dp_T0[arg1];
4463                } else {
4464                    if (coval > 0.0)
4465                        MINDEC(ret_c,2);
4466                    if (dp_T0[arg] == 0)
4467                        MINDEC(ret_c,0);
4468                    dp_T0[res] = dp_T0[arg2];
4469                }
4470                FOR_0_LE_l_LT_pk
4471                TRES_INC = TARG1_INC | TARG2_INC;
4472#else
4473#if !defined(_ZOS_) /* BREAK_ZOS */
4474                if (dp_T0[arg] > 0)
4475                    FOR_0_LE_l_LT_pk
4476                    TRES_INC = TARG1_INC;
4477                else
4478                    FOR_0_LE_l_LT_pk
4479                    TRES_INC = TARG2_INC;
4480#endif
4481
4482                if (dp_T0[arg] > 0) {
4483                    if (coval <= 0.0)
4484                        MINDEC(ret_c,2);
4485                    dp_T0[res] = dp_T0[arg1];
4486                } else {
4487                    if (coval > 0.0)
4488                        MINDEC(ret_c,2);
4489                    if (dp_T0[arg] == 0)
4490                        MINDEC(ret_c,0);
4491                    dp_T0[res] = dp_T0[arg2];
4492                }
4493#endif
4494#endif
4495#else
4496                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4497                exit(-2);
4498#endif /* ALL_TOGETHER_AGAIN */
4499                }
4500                break;
4501
4502            case ref_cond_assign_s:                                  /* cond_assign_s */
4503                arg   = get_locint_f();
4504                arg1  = get_locint_f();
4505                arg2   = get_locint_f();
4506                coval = get_val_f();
4507
4508#if !defined(_NTIGHT_)
4509                res = (size_t)trunc(fabs(dp_T0[arg2]));
4510                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4511
4512                /* olvo 980924 changed order to allow reflexive ops */
4513#if defined(_INDO_)
4514                if (dp_T0[arg] > 0) {
4515#if defined(_INDOPRO_)
4516                    copy_index_domain(res, arg1, ind_dom);
4517#endif
4518#if defined(_NONLIND_)
4519                    arg_index[res] = arg_index[arg1];
4520#endif
4521                }
4522#else
4523#if !defined(_ZOS_) /* BREAK_ZOS */
4524                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
4525                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
4526#endif /* ALL_TOGETHER_AGAIN */
4527
4528#ifdef _INT_FOR_
4529                coval = get_val_f();
4530
4531                if (dp_T0[arg] > 0)
4532                    FOR_0_LE_l_LT_pk
4533                    TRES_INC = TARG1_INC;
4534
4535                if (dp_T0[arg] > 0) {
4536                    if (coval <= 0.0)
4537                        MINDEC(ret_c,2);
4538                    dp_T0[res] = dp_T0[arg1];
4539                } else
4540                    if (dp_T0[arg] == 0)
4541                        MINDEC(ret_c,0);
4542#else
4543#if !defined(_ZOS_) /* BREAK_ZOS */
4544                if (dp_T0[arg] > 0)
4545                    FOR_0_LE_l_LT_pk
4546                    TRES_INC = TARG1_INC;
4547#endif
4548                if (dp_T0[arg] > 0) {
4549                    if (coval <= 0.0)
4550                        MINDEC(ret_c,2);
4551                    dp_T0[res] = dp_T0[arg1];
4552                } else
4553                    if (dp_T0[arg] == 0)
4554                        MINDEC(ret_c,0);
4555#endif
4556#endif
4557#else
4558                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4559                exit(-2);
4560#endif /* ALL_TOGETHER_AGAIN */
4561                break;
4562
4563                /****************************************************************************/
4564                /*                                                          REMAINING STUFF */
4565
4566                /*--------------------------------------------------------------------------*/
4567            case take_stock_op:                                  /* take_stock_op */
4568                size = get_locint_f();
4569                res  = get_locint_f();
4570                d    = get_val_v_f(size);
4571
4572                for (ls=0;ls<size;ls++) {
4573#if !defined(_NTIGHT_)
4574                    dp_T0[res]=*d;
4575#endif /* !_NTIGHT_ */
4576#if !defined(_INDO_)
4577#if !defined(_ZOS_) /* BREAK_ZOS */
4578                    ASSIGN_T(Tres,TAYLOR_BUFFER[res])
4579
4580                    FOR_0_LE_l_LT_pk
4581                    TRES_INC = 0;
4582
4583#endif /* ALL_TOGETHER_AGAIN */
4584                    res++;
4585#if !defined(_NTIGHT_)
4586                    d++;
4587#endif /* !_NTIGHT_ */
4588#endif
4589                }
4590                break;
4591
4592                /*--------------------------------------------------------------------------*/
4593            case death_not:                                          /* death_not */
4594                arg1=get_locint_f();
4595                arg2=get_locint_f();
4596
4597#ifdef _KEEP_
4598                if (keep) {
4599                    do {
4600                        IF_KEEP_WRITE_TAYLOR(arg2,keep,k,p)
4601                    } while(arg1 < arg2-- );
4602                }
4603#endif
4604                break;
4605
4606                /*--------------------------------------------------------------------------*/
4607#if defined(_EXTERN_) /* ZOS,  FOS, FOV up to now */
4608            case ext_diff:                       /* extern differntiated function */
4609                ADOLC_CURRENT_TAPE_INFOS.ext_diff_fct_index=get_locint_f();
4610                n=get_locint_f();
4611                m=get_locint_f();
4612                ADOLC_CURRENT_TAPE_INFOS.lowestXLoc_for = get_locint_f();
4613                ADOLC_CURRENT_TAPE_INFOS.lowestYLoc_for = get_locint_f();
4614                ADOLC_CURRENT_TAPE_INFOS.cpIndex = get_locint_f();
4615                edfct=get_ext_diff_fct(ADOLC_CURRENT_TAPE_INFOS.ext_diff_fct_index);
4616
4617                if (edfct->ADOLC_EXT_FCT_POINTER==NULL)
4618                    fail(ADOLC_EXT_DIFF_NULLPOINTER_DIFFFUNC);
4619                if (n>0) {
4620                    if (edfct->dp_x==NULL) fail(ADOLC_EXT_DIFF_NULLPOINTER_ARGUMENT);
4621#if !defined(_ZOS_)
4622                    if (ADOLC_EXT_POINTER_X==NULL) fail(ADOLC_EXT_DIFF_NULLPOINTER_ARGUMENT);
4623#endif
4624                }
4625                if (m>0) {
4626                    if (edfct->dp_y==NULL) fail(ADOLC_EXT_DIFF_NULLPOINTER_ARGUMENT);
4627#if !defined(_ZOS_)
4628                    if (ADOLC_EXT_POINTER_Y==NULL) fail(ADOLC_EXT_DIFF_NULLPOINTER_ARGUMENT);
4629#endif
4630                }
4631
4632                arg = ADOLC_CURRENT_TAPE_INFOS.lowestXLoc_for;
4633                for (loop=0; loop<n; ++loop) {
4634                    edfct->dp_x[loop]=dp_T0[arg];
4635#if !defined(_ZOS_)
4636                    ADOLC_EXT_LOOP
4637                        ADOLC_EXT_POINTER_X[loop]ADOLC_EXT_SUBSCRIPT =
4638                        TAYLOR_BUFFER[arg]ADOLC_EXT_SUBSCRIPT;
4639#endif
4640                    ++arg;
4641                }
4642                arg = ADOLC_CURRENT_TAPE_INFOS.lowestYLoc_for;
4643                for (loop=0; loop<m; ++loop) {
4644                    edfct->dp_y[loop]=dp_T0[arg];
4645#if !defined(_ZOS_)
4646                    ADOLC_EXT_LOOP
4647                        ADOLC_EXT_POINTER_Y[loop]ADOLC_EXT_SUBSCRIPT =
4648                        TAYLOR_BUFFER[arg]ADOLC_EXT_SUBSCRIPT;
4649#endif
4650                    ++arg;
4651                }
4652
4653                ext_retc = edfct->ADOLC_EXT_FCT_COMPLETE;
4654                MINDEC(ret_c, ext_retc);
4655
4656                res = ADOLC_CURRENT_TAPE_INFOS.lowestXLoc_for;
4657                for (loop=0; loop<n; ++loop) {
4658                    IF_KEEP_WRITE_TAYLOR(res, keep, k, p);
4659                    dp_T0[res]=edfct->dp_x[loop];
4660#if !defined(_ZOS_)
4661                    ADOLC_EXT_LOOP
4662                        TAYLOR_BUFFER[res]ADOLC_EXT_SUBSCRIPT =
4663                        ADOLC_EXT_POINTER_X[loop]ADOLC_EXT_SUBSCRIPT;
4664#endif
4665                    ++res;
4666                }
4667                res = ADOLC_CURRENT_TAPE_INFOS.lowestYLoc_for;
4668                for (loop=0; loop<m; ++loop) {
4669                    IF_KEEP_WRITE_TAYLOR(res, keep, k, p);
4670                    dp_T0[res]=edfct->dp_y[loop];
4671#if !defined(_ZOS_)
4672                    ADOLC_EXT_LOOP
4673                        TAYLOR_BUFFER[res]ADOLC_EXT_SUBSCRIPT =
4674                        ADOLC_EXT_POINTER_Y[loop]ADOLC_EXT_SUBSCRIPT;
4675#endif
4676                    ++res;
4677                }
4678
4679                break;
4680#endif
4681                /*--------------------------------------------------------------------------*/
4682
4683            default:                                                   /* default */
4684                /* Die here, we screwed up */
4685
4686                fprintf(DIAG_OUT,"ADOL-C fatal error in " GENERATED_FILENAME " ("
4687                        __FILE__
4688                        ") : no such operation %d\n", operation);
4689                exit(-1);
4690                break;
4691
4692        } /* endswitch */
4693
4694        /* Read the next operation */
4695        operation=get_op_f();
4696#if defined(ADOLC_DEBUG)
4697        ++countPerOperation[operation];
4698#endif /* ADOLC_DEBUG */
4699    }  /* endwhile */
4700
4701
4702#if defined(ADOLC_DEBUG)
4703    printf("\nTape contains:\n");
4704    for (v = 0; v < 256; ++v)
4705        if (countPerOperation[v] > 0)
4706            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]);
4707    printf("\n");
4708#endif /* ADOLC_DEBUG */
4709
4710#if defined(_KEEP_)
4711    if (keep) taylor_close(taylbuf);
4712#endif
4713
4714    /* clean up */
4715#if !defined (_NTIGHT_)
4716    free(dp_T0);
4717#endif /* !_NTIGHT_ */
4718#if !defined(_INDO_)
4719#if !defined(_ZOS_)
4720#   if defined(_FOS_)
4721    free(dp_T);
4722#   else
4723#if !defined (_INT_FOR_)
4724    myfree2(dpp_T);
4725    free(dp_Ttemp);
4726#endif /* !_NTIGHT_ */
4727#endif
4728#endif
4729#endif
4730#if defined(_HIGHER_ORDER_)
4731    free(dp_z);
4732#endif
4733
4734    end_sweep();
4735
4736
4737#if defined(_INDO_)
4738#if defined(_INDOPRO_)
4739    for(i=0;i<max_ind_dom;i++)
4740      {
4741        free(ind_dom[i]);
4742      }
4743    free(ind_dom);
4744#endif
4745#if defined(_NONLIND_)
4746    for( i=0; i < indcheck; i++) {
4747      traverse_crs(&nonl_dom[i],&sod[i],indcheck+1);
4748      free_tree(&nonl_dom[i],indcheck+1);
4749      crs[i] = (unsigned int*) malloc(sizeof(unsigned int) * (sod[i].entry+1));
4750      crs[i][0] = sod[i].entry;
4751      temp = sod[i].left;
4752      for( ii=1; ii <=sod[i].entry; ii++)
4753        {
4754          crs[i][ii] = temp->entry;
4755          temp1 = temp->left;
4756          free(temp);
4757          temp = temp1;
4758        }
4759    }
4760
4761    free(sod);
4762    free(nonl_dom);
4763    free(fod);
4764    free(arg_index);
4765
4766#endif
4767#endif
4768    return ret_c;
4769}
4770
4771/****************************************************************************/
4772
4773#if defined(_INDOPRO_)
4774
4775/****************************************************************************/
4776/* set operations for propagation of index domains                          */
4777
4778/*--------------------------------------------------------------------------*/
4779/* operations on index domains                                              */
4780
4781#if defined(_TIGHT_)
4782void copy_index_domain(int res, int arg, locint **ind_dom) {
4783
4784   int i;
4785
4786   if (ind_dom[arg][0] > ind_dom[res][1])
4787     {
4788       free(ind_dom[res]);
4789       ind_dom[res] = (locint *)  malloc(sizeof(locint) * 2*(ind_dom[arg][0]+1));
4790       ind_dom[res][1] = 2*ind_dom[arg][0];
4791     }
4792
4793
4794    for(i=2;i<ind_dom[arg][0]+2;i++)
4795       ind_dom[res][i] = ind_dom[arg][i];
4796    ind_dom[res][0] = ind_dom[arg][0];
4797}
4798
4799
4800void merge_2_index_domains(int res, int arg, locint **ind_dom)
4801{
4802
4803  int num,num1,num2, i,j,k,l;
4804  locint *temp_array, *arg_ind_dom, *res_ind_dom;
4805
4806  if (ind_dom[res][0] == 0)
4807    copy_index_domain(res,arg,ind_dom);
4808  else
4809    {
4810      if (res != arg)
4811     {
4812       arg_ind_dom = ind_dom[arg];
4813       res_ind_dom = ind_dom[res];
4814
4815       num  = ind_dom[res][0];
4816       num1 = arg_ind_dom[0];
4817       num2 = ind_dom[res][1];
4818
4819       if (num2 < num1+num)
4820         num2 = num1+num;
4821
4822       temp_array = (locint *)  malloc(sizeof(locint)* (num2+2));
4823       temp_array[1] = num2;
4824
4825       i = 2;
4826       j = 2;
4827       k = 2;
4828       num += 2;
4829       num1 += 2;
4830       while ((i< num) && (j < num1))
4831         {
4832           if (res_ind_dom[i] < arg_ind_dom[j])
4833          {
4834            temp_array[k] = res_ind_dom[i];
4835            i++; k++;
4836          }
4837           else
4838          {
4839            if (res_ind_dom[i] == arg_ind_dom[j])
4840              {
4841                temp_array[k] = arg_ind_dom[j];
4842                i++;j++;k++;
4843              }
4844            else
4845              {
4846                temp_array[k] = arg_ind_dom[j];
4847                j++;k++;
4848              }
4849          }
4850         }
4851       for(l = i;l<num;l++)
4852         {
4853           temp_array[k] = res_ind_dom[l];
4854           k++;
4855         }
4856       for(l = j;l<num1;l++)
4857         {
4858           temp_array[k] = arg_ind_dom[l];
4859           k++;
4860         }
4861       temp_array[0] = k-2;
4862       free(ind_dom[res]);
4863       ind_dom[res]=temp_array;
4864     }
4865    }
4866
4867
4868}
4869
4870void combine_2_index_domains(int res, int arg1, int arg2, locint **ind_dom) {
4871
4872    if (res != arg1)
4873       copy_index_domain(res, arg1, ind_dom);
4874
4875    merge_2_index_domains(res, arg2, ind_dom);
4876}
4877
4878void merge_3_index_domains(int res, int arg1, int arg2, locint **ind_dom) {
4879    merge_2_index_domains(res, arg1, ind_dom);
4880    merge_2_index_domains(res, arg2, ind_dom);
4881}
4882
4883
4884
4885#endif
4886#endif
4887
4888
4889#if defined(_NONLIND_)
4890#if defined(_TIGHT_)
4891
4892void free_tree(IndexElement* tree, int num)
4893{
4894
4895  if (tree->left != NULL)
4896    {
4897      free_tree(tree->left,num);
4898    }
4899  if (tree->right != NULL)
4900    {
4901      free_tree(tree->right,num);
4902     }
4903    {
4904      if (tree->entry == num)
4905        free(tree);
4906
4907    }
4908 
4909}
4910void traverse_crs(IndexElement* tree,  IndexElement_sod* sod, int num)
4911{
4912
4913  IndexElement_sod *temp, *temp1;
4914  int ii;
4915
4916  if (tree->left != NULL)
4917    {
4918      traverse_crs(tree->left, sod, num);
4919    }
4920  if (tree->right != NULL)
4921    {
4922      traverse_crs(tree->right, sod, num);
4923    }
4924  if (tree->entry < num)
4925    {
4926      temp = sod->left;
4927      if (temp == NULL)
4928        {
4929          temp = (struct IndexElement_sod*) malloc(sizeof(struct IndexElement_sod));
4930          temp->left = NULL;
4931          temp->entry = tree->entry;
4932          sod->entry++;
4933          sod->left=temp;
4934        }
4935      else
4936        {
4937          while ((temp->entry < tree->entry) && (temp->left != NULL))
4938            {
4939              temp1 = temp;
4940              temp = temp->left;
4941            }
4942          if (temp->left == NULL)
4943            {
4944              if(temp->entry < tree->entry)
4945                {
4946                  temp->left = (struct IndexElement_sod*) malloc(sizeof(struct IndexElement_sod));
4947                  temp = temp->left;
4948                  temp->left = NULL;
4949                  temp->entry = tree->entry;
4950                  sod->entry++;
4951                }
4952              if(temp->entry > tree->entry)
4953                {
4954                  temp->left = (struct IndexElement_sod*) malloc(sizeof(struct IndexElement_sod));
4955                  temp->left->entry = temp->entry;
4956                  temp->left->left = NULL;
4957                  temp->entry = tree->entry;
4958                  sod->entry++;
4959                }
4960            }
4961          else
4962            {
4963              if (temp->entry > tree->entry)
4964                {
4965                  temp1 = (struct IndexElement_sod*) malloc(sizeof(struct IndexElement_sod));
4966                  temp1->left = temp->left;
4967                  temp1->entry = temp->entry;
4968                  temp->entry = tree->entry;
4969                  temp->left=temp1;
4970                  sod->entry++;
4971                }
4972             
4973            }
4974        }
4975    }
4976}
4977
4978void traverse_unary(IndexElement* tree,  IndexElement* nonl_dom,  IndexElement* fodi, int num, int maxopind)
4979{
4980  IndexElement *temp;
4981
4982  if (tree->left != NULL)
4983    {
4984      traverse_unary(tree->left, nonl_dom, fodi, num, maxopind);
4985      if (tree->right != NULL)
4986        {
4987          traverse_unary(tree->right, nonl_dom, fodi, num, maxopind);
4988        }
4989     }
4990  else
4991    {
4992      if(tree->entry<maxopind)
4993        {
4994          temp = (struct IndexElement*) malloc(sizeof(struct IndexElement));
4995          temp->right = fodi;
4996          temp->left = nonl_dom[tree->entry].left;
4997          temp->entry= num;
4998          nonl_dom[tree->entry].left = temp;
4999        }
5000    }
5001}
5002
5003#endif
5004#endif
5005END_C_DECLS
Note: See TracBrowser for help on using the repository browser.