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

Last change on this file since 370 was 370, checked in by kulshres, 8 years ago

Merge branch '2.3.x_ISSM' into svn

This introduces the new externally differentiated functions API

From: Jean Utke <utke@…>

Please see comments in ADOL-C/include/adolc/externfcts.h for details

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

  • Property svn:keywords set to Author Date Id Revision
File size: 155.2 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 370 2012-11-22 13:18:52Z 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      switch (operation) {
1069
1070
1071                /****************************************************************************/
1072                /*                                                                  MARKERS */
1073
1074                /*--------------------------------------------------------------------------*/
1075            case end_of_op:                                          /* end_of_op */
1076                get_op_block_f();
1077                operation=get_op_f();
1078                /* Skip next operation, it's another end_of_op */
1079                break;
1080
1081                /*--------------------------------------------------------------------------*/
1082            case end_of_int:                                        /* end_of_int */
1083                get_loc_block_f();
1084                break;
1085
1086                /*--------------------------------------------------------------------------*/
1087            case end_of_val:                                        /* end_of_val */
1088               get_val_block_f();
1089                break;
1090                /*--------------------------------------------------------------------------*/
1091            case start_of_tape:                                  /* start_of_tape */
1092            case end_of_tape:                                      /* end_of_tape */
1093                break;
1094
1095
1096                /****************************************************************************/
1097                /*                                                               COMPARISON */
1098
1099                /*--------------------------------------------------------------------------*/
1100            case eq_zero:                                              /* eq_zero */
1101                arg = get_locint_f();
1102
1103#if !defined(_NTIGHT_)
1104                if (dp_T0[arg] != 0) {
1105                    if (ADOLC_GLOBAL_TAPE_VARS.branchSwitchWarning)
1106                        fprintf(DIAG_OUT,
1107                                "ADOL-C Warning: Branch switch detected in comparison "
1108                                "(operator eq_zero).\n"
1109                                "Forward sweep aborted! Retaping recommended!\n");
1110                    ret_c = -1;
1111                    operation = end_of_tape;
1112                    continue;
1113                }
1114                ret_c = 0;
1115#endif /* !_NTIGHT_ */
1116                break;
1117
1118                /*--------------------------------------------------------------------------*/
1119            case neq_zero:                                            /* neq_zero */
1120                arg = get_locint_f();
1121
1122#if !defined(_NTIGHT_)
1123                if (dp_T0[arg] == 0) {
1124                    if (ADOLC_GLOBAL_TAPE_VARS.branchSwitchWarning)
1125                        fprintf(DIAG_OUT,
1126                                "ADOL-C Warning: Branch switch detected in comparison "
1127                                "(operator neq_zero).\n"
1128                                "Forward sweep aborted! Retaping recommended!\n");
1129                    ret_c = -1;
1130                    operation = end_of_tape;
1131                    continue;
1132                }
1133#endif /* !_NTIGHT_ */
1134                break;
1135
1136                /*--------------------------------------------------------------------------*/
1137            case le_zero:                                              /* le_zero */
1138                arg = get_locint_f();
1139
1140#if !defined(_NTIGHT_)
1141                if (dp_T0[arg] > 0) {
1142                    if (ADOLC_GLOBAL_TAPE_VARS.branchSwitchWarning)
1143                        fprintf(DIAG_OUT,
1144                                "ADOL-C Warning: Branch switch detected in comparison "
1145                                "(operator le_zero).\n"
1146                                "Forward sweep aborted! Retaping recommended!\n");
1147                    ret_c = -1;
1148                    operation = end_of_tape;
1149                    continue;
1150                }
1151                if (dp_T0[arg] == 0)
1152                    ret_c = 0;
1153#endif /* !_NTIGHT_ */
1154                break;
1155
1156                /*--------------------------------------------------------------------------*/
1157            case gt_zero:                                              /* gt_zero */
1158                arg = get_locint_f();
1159
1160#if !defined(_NTIGHT_)
1161                if (dp_T0[arg] <= 0) {
1162                    if (ADOLC_GLOBAL_TAPE_VARS.branchSwitchWarning)
1163                        fprintf(DIAG_OUT,
1164                                "ADOL-C Warning: Branch switch detected in comparison "
1165                                "(operator gt_zero).\n"
1166                                "Forward sweep aborted! Retaping recommended!\n");
1167                    ret_c = -1;
1168                    operation = end_of_tape;
1169                    continue;
1170                }
1171#endif /* !_NTIGHT_ */
1172                break;
1173
1174                /*--------------------------------------------------------------------------*/
1175            case ge_zero:                                              /* ge_zero */
1176                arg = get_locint_f();
1177
1178#if !defined(_NTIGHT_)
1179                if (dp_T0[arg] < 0) {
1180                    if (ADOLC_GLOBAL_TAPE_VARS.branchSwitchWarning)
1181                        fprintf(DIAG_OUT,
1182                                "ADOL-C Warning: Branch switch detected in comparison "
1183                                "(operator ge_zero).\n"
1184                                "Forward sweep aborted! Retaping recommended!\n");
1185                    ret_c = -1;
1186                    operation = end_of_tape;
1187                    continue;
1188                }
1189                if (dp_T0[arg] == 0)
1190                    ret_c = 0;
1191#endif /* !_NTIGHT_ */
1192                break;
1193
1194                /*--------------------------------------------------------------------------*/
1195            case lt_zero:                                              /* lt_zero */
1196                arg = get_locint_f();
1197
1198#if !defined(_NTIGHT_)
1199                if (dp_T0[arg] >= 0) {
1200                    if (ADOLC_GLOBAL_TAPE_VARS.branchSwitchWarning)
1201                        fprintf(DIAG_OUT,
1202                                "ADOL-C Warning: Branch switch detected in comparison "
1203                                "(operator lt_zero).\n"
1204                                "Forward sweep aborted! Retaping recommended!\n");
1205                    ret_c = -1;
1206                    operation = end_of_tape;
1207                    continue;
1208                }
1209#endif /* !_NTIGHT_ */
1210                break;
1211
1212
1213                /****************************************************************************/
1214                /*                                                              ASSIGNMENTS */
1215
1216                /*--------------------------------------------------------------------------*/
1217            case assign_a:           /* assign an adouble variable an    assign_a */
1218                /* adouble value. (=) */
1219                arg = get_locint_f();
1220                res = get_locint_f();
1221
1222                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1223
1224
1225#if !defined(_NTIGHT_)
1226                dp_T0[res] = dp_T0[arg];
1227#endif /* !_NTIGHT_ */
1228
1229#if defined(_INDO_)
1230#if defined(_INDOPRO_)
1231                copy_index_domain(res, arg, ind_dom);
1232#endif
1233#if defined(_NONLIND_)
1234                arg_index[res] = arg_index[arg];               
1235#endif           
1236#else
1237#if !defined(_ZOS_) /* BREAK_ZOS */
1238                ASSIGN_T(Targ,TAYLOR_BUFFER[arg])
1239                ASSIGN_T(Tres,TAYLOR_BUFFER[res])
1240
1241                FOR_0_LE_l_LT_pk
1242                TRES_INC = TARG_INC;
1243#endif
1244#endif /* ALL_TOGETHER_AGAIN */
1245                break;
1246
1247                /*--------------------------------------------------------------------------*/
1248            case assign_d:            /* assign an adouble variable a    assign_d */
1249                /* double value. (=) */
1250                res   = get_locint_f();
1251                coval = get_val_f();
1252
1253                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1254
1255#if !defined(_NTIGHT_)
1256                dp_T0[res] = coval;
1257#endif /* !_NTIGHT_ */
1258
1259#if defined(_INDO_)
1260#if defined(_INDOPRO_)
1261                ind_dom[res][0]=0;
1262#endif
1263#if defined(_NONLIND_)
1264                fod[opind].entry = maxopind+2;
1265                fod[opind].left = NULL;
1266                fod[opind].right = NULL;
1267                arg_index[res] = opind++;               
1268#endif
1269#else
1270#if !defined(_ZOS_) /* BREAK_ZOS */
1271                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1272
1273                FOR_0_LE_l_LT_pk
1274                TRES_INC = 0;
1275#endif
1276#endif /* ALL_TOGETHER_AGAIN */
1277                break;
1278
1279                /*--------------------------------------------------------------------------*/
1280            case assign_d_zero:  /* assign an adouble variable a    assign_d_zero */
1281                /* double value. (0) (=) */
1282                res   = get_locint_f();
1283
1284                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1285
1286#if !defined(_NTIGHT_)
1287                dp_T0[res] = 0.0;
1288#endif /* !_NTIGHT_ */
1289
1290#if defined(_INDO_)
1291#if defined(_INDOPRO_)
1292                ind_dom[res][0]=0;
1293#endif
1294#if defined(_NONLIND_)
1295                fod[opind].entry = maxopind+2;
1296                fod[opind].left = NULL;
1297                fod[opind].right = NULL;
1298                arg_index[res] = opind++;               
1299#endif
1300#else
1301#if !defined(_ZOS_) /* BREAK_ZOS */
1302                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1303
1304                FOR_0_LE_l_LT_pk
1305                TRES_INC = 0;
1306#endif
1307#endif /* ALL_TOGETHER_AGAIN */
1308                break;
1309
1310                /*--------------------------------------------------------------------------*/
1311            case assign_d_one:    /* assign an adouble variable a    assign_d_one */
1312                /* double value. (1) (=) */
1313                res   = get_locint_f();
1314
1315                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1316
1317#if !defined(_NTIGHT_)
1318                dp_T0[res] = 1.0;
1319#endif /* !_NTIGHT_ */
1320
1321#if defined(_INDO_)
1322#if defined(_INDOPRO_)
1323                ind_dom[res][0]=0;
1324#endif
1325#if defined(_NONLIND_)
1326                fod[opind].entry = maxopind+2;
1327                fod[opind].left = NULL;
1328                fod[opind].right = NULL;
1329                arg_index[res] = opind++;               
1330#endif
1331#else
1332#if !defined(_ZOS_) /* BREAK_ZOS */
1333                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1334
1335                FOR_0_LE_l_LT_pk
1336                TRES_INC = 0;
1337
1338#endif
1339#endif /* ALL_TOGETHER_AGAIN */
1340                break;
1341
1342                /*--------------------------------------------------------------------------*/
1343            case assign_ind:       /* assign an adouble variable an    assign_ind */
1344                /* independent double value (<<=) */
1345                res = get_locint_f();
1346
1347                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1348
1349#if !defined(_NTIGHT_)
1350                dp_T0[res] = basepoint[indexi];
1351#endif /* !_NTIGHT_ */
1352
1353#if defined(_INDO_)
1354#if defined(_INDOPRO_)
1355                ind_dom[res][0] = 1;
1356                ind_dom[res][2] = indexi;
1357#endif         
1358#if defined(_NONLIND_)
1359                fod[opind].entry = indexi;
1360                fod[opind].left = NULL;
1361                fod[opind].right = NULL;
1362                arg_index[res] = opind++;               
1363#endif         
1364#else
1365#if !defined(_ZOS_) /* BREAK_ZOS */
1366                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1367
1368#ifdef _INT_FOR_
1369                FOR_0_LE_l_LT_p
1370                TRES_INC = ARGUMENT(indexi,l,i);
1371#else
1372                FOR_0_LE_l_LT_p
1373                FOR_0_LE_i_LT_k
1374                TRES_INC = ARGUMENT(indexi,l,i);
1375#endif
1376#endif
1377#endif /* ALL_TOGETHER_AGAIN */
1378                ++indexi;
1379                break;
1380
1381                /*--------------------------------------------------------------------------*/
1382            case assign_dep:           /* assign a float variable a    assign_dep */
1383                /* dependent adouble value. (>>=) */
1384                res = get_locint_f();
1385
1386#if !defined(_INDO_)
1387#if !defined(_NTIGHT_)
1388                if ( valuepoint != NULL )
1389                  valuepoint[indexd] = dp_T0[res];
1390#endif /* !_NTIGHT_ */
1391#endif
1392
1393#if defined(_INDO_)
1394#if defined(_INDOPRO_)
1395          if (ind_dom[res][0] != 0) {
1396            crs[indexd] = (unsigned int*) malloc(sizeof(unsigned int) * (ind_dom[res][0]+1));
1397            crs[indexd][0] = ind_dom[res][0];
1398            for(l=1;l<=crs[indexd][0];l++) {
1399              crs[indexd][l] = ind_dom[res][l+1];
1400            }
1401          }
1402          else {
1403            crs[indexd] = (unsigned int*) malloc(sizeof(unsigned int));
1404            crs[indexd][0] =0;
1405          }
1406#endif
1407#else
1408#if !defined(_ZOS_) /* BREAK_ZOS */
1409                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1410
1411#ifdef _INT_FOR_
1412                if (taylors != 0 )  /* ??? question: why here? */
1413                    FOR_0_LE_l_LT_p
1414                    TAYLORS(indexd,l,i) = TRES_INC;
1415#else
1416                if (taylors != 0 )  /* ??? question: why here? */
1417                    FOR_0_LE_l_LT_p
1418                    FOR_0_LE_i_LT_k
1419                    TAYLORS(indexd,l,i) = TRES_INC;
1420#endif
1421#endif
1422#endif /* ALL_TOGETHER_AGAIN */
1423                indexd++;
1424                break;
1425
1426
1427                /****************************************************************************/
1428                /*                                                   OPERATION + ASSIGNMENT */
1429
1430                /*--------------------------------------------------------------------------*/
1431            case eq_plus_d:            /* Add a floating point to an    eq_plus_d */
1432                /* adouble. (+=) */
1433                res   = get_locint_f();
1434                coval = get_val_f();
1435
1436                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1437
1438#if !defined(_NTIGHT_)
1439                dp_T0[res] += coval;
1440#endif /* !_NTIGHT_ */
1441                break;
1442
1443                /*--------------------------------------------------------------------------*/
1444            case eq_plus_a:             /* Add an adouble to another    eq_plus_a */
1445                /* adouble. (+=) */
1446                arg = get_locint_f();
1447                res = get_locint_f();
1448
1449                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1450
1451#if !defined(_NTIGHT_)
1452                dp_T0[res] += dp_T0[arg];
1453#endif /* !_NTIGHT_ */
1454
1455#if defined(_INDO_)
1456#if defined(_INDOPRO_)
1457                merge_2_index_domains(res, arg, ind_dom);
1458#endif
1459#if defined(_NONLIND_)
1460                fod[opind].entry = maxopind+2;
1461                fod[opind].left = &fod[arg_index[res]];
1462                fod[opind].right = &fod[arg_index[arg]];
1463                arg_index[res] = opind++;               
1464#endif
1465#else
1466#if !defined(_ZOS_) /* BREAK_ZOS */
1467                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1468                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
1469
1470#ifdef _INT_FOR_
1471                FOR_0_LE_l_LT_pk
1472                TRES_INC |= TARG_INC;
1473#else
1474                FOR_0_LE_l_LT_pk
1475                TRES_INC += TARG_INC;
1476#endif
1477#endif
1478#endif /* ALL_TOGETHER_AGAIN */
1479                break;
1480
1481                /*--------------------------------------------------------------------------*/
1482            case eq_min_d:       /* Subtract a floating point from an    eq_min_d */
1483                /* adouble. (-=) */
1484                res = get_locint_f();
1485                coval = get_val_f();
1486
1487                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1488
1489#if !defined(_NTIGHT_)
1490                dp_T0[res] -= coval;
1491#endif /* !_NTIGHT_ */
1492                break;
1493
1494                /*--------------------------------------------------------------------------*/
1495            case eq_min_a:        /* Subtract an adouble from another    eq_min_a */
1496                /* adouble. (-=) */
1497                arg = get_locint_f();
1498                res = get_locint_f();
1499
1500                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1501
1502#if !defined(_NTIGHT_)
1503                dp_T0[res] -= dp_T0[arg];
1504#endif /* !_NTIGHT_ */
1505
1506#if defined(_INDO_)
1507#if defined(_INDOPRO_)
1508                merge_2_index_domains(res, arg, ind_dom);
1509#endif
1510#if defined(_NONLIND_)
1511                fod[opind].entry = maxopind+2;
1512                fod[opind].left = &fod[arg_index[res]];
1513                fod[opind].right = &fod[arg_index[arg]];
1514                arg_index[res] = opind++;               
1515#endif
1516#else
1517#if !defined(_ZOS_) /* BREAK_ZOS */
1518                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1519                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
1520
1521#ifdef _INT_FOR_
1522                FOR_0_LE_l_LT_pk
1523                TRES_INC |= TARG_INC;
1524#else
1525                FOR_0_LE_l_LT_pk
1526                TRES_INC -= TARG_INC;
1527#endif
1528#endif
1529#endif /* ALL_TOGETHER_AGAIN */
1530                break;
1531
1532                /*--------------------------------------------------------------------------*/
1533            case eq_mult_d:              /* Multiply an adouble by a    eq_mult_d */
1534                /* flaoting point. (*=) */
1535                res   = get_locint_f();
1536                coval = get_val_f();
1537
1538                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1539
1540#if !defined(_NTIGHT_)
1541                dp_T0[res] *= coval;
1542#endif /* !_NTIGHT_ */
1543
1544#if !defined(_INDO_)
1545#if !defined(_ZOS_) /* BREAK_ZOS */
1546#if !defined( _INT_FOR_)
1547
1548                FOR_0_LE_l_LT_pk
1549                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1550
1551                FOR_0_LE_l_LT_pk
1552                TRES_INC *= coval;
1553#endif
1554#endif
1555#endif /* ALL_TOGETHER_AGAIN */
1556                break;
1557
1558                /*--------------------------------------------------------------------------*/
1559            case eq_mult_a:       /* Multiply one adouble by another    eq_mult_a */
1560                /* (*=) */
1561                arg = get_locint_f();
1562                res = get_locint_f();
1563
1564                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1565
1566#if defined(_INDO_)
1567#if defined(_INDOPRO_)
1568                merge_2_index_domains(res, arg, ind_dom);
1569#endif
1570#if defined(_NONLIND_)
1571                fod[opind].entry = maxopind+2;
1572                fod[opind].left = &fod[arg_index[res]];
1573                fod[opind].right = &fod[arg_index[arg]];
1574                traverse_unary(&fod[arg_index[res]], nonl_dom, &fod[arg_index[arg]], indcheck+1,maxopind+2);
1575                traverse_unary(&fod[arg_index[arg]], nonl_dom, &fod[arg_index[res]], indcheck+1,maxopind+2);
1576                arg_index[res] = opind++;               
1577#endif
1578#else
1579#if !defined(_ZOS_) /* BREAK_ZOS */
1580                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1581                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
1582
1583                INC_pk_1(Tres)
1584                INC_pk_1(Targ)
1585
1586#ifdef _INT_FOR_
1587                FOR_p_GT_l_GE_0
1588                TRES_FODEC |= TARG_DEC;
1589#else
1590                FOR_p_GT_l_GE_0
1591                FOR_k_GT_i_GE_0
1592                { TRES_FODEC = dp_T0[res]*TARG_DEC +
1593                               TRES*dp_T0[arg];
1594                  DEC_TRES_FO
1595#ifdef _HIGHER_ORDER_
1596                  TresOP = Tres-i;
1597                  TargOP = Targ;
1598
1599                  for (j=0;j<i;j++)
1600                  *Tres += (*TresOP++) * (*TargOP--);
1601                  Tres--;
1602#endif /* _HIGHER_ORDER_ */
1603                }
1604#endif
1605#endif
1606#endif /* ALL_TOGETHER_AGAIN */
1607#if !defined(_NTIGHT_)
1608               dp_T0[res] *= dp_T0[arg];
1609#endif /* !_NTIGHT_ */
1610                break;
1611
1612                /*--------------------------------------------------------------------------*/
1613            case incr_a:                        /* Increment an adouble    incr_a */
1614                res   = get_locint_f();
1615
1616                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1617
1618#if !defined(_NTIGHT_)
1619                dp_T0[res]++;
1620#endif /* !_NTIGHT_ */
1621                break;
1622
1623                /*--------------------------------------------------------------------------*/
1624            case decr_a:                        /* Increment an adouble    decr_a */
1625                res   = get_locint_f();
1626
1627                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1628
1629#if !defined(_NTIGHT_)
1630                dp_T0[res]--;
1631#endif /* !_NTIGHT_ */
1632                break;
1633
1634
1635                /****************************************************************************/
1636                /*                                                        BINARY OPERATIONS */
1637
1638                /*--------------------------------------------------------------------------*/
1639            case plus_a_a:                 /* : Add two adoubles. (+)    plus a_a */
1640                arg1 = get_locint_f();
1641                arg2 = get_locint_f();
1642                res  = get_locint_f();
1643
1644                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1645
1646#if !defined(_NTIGHT_)
1647                dp_T0[res] = dp_T0[arg1] + dp_T0[arg2];
1648#endif /* !_NTIGHT_ */
1649
1650#if defined(_INDO_)
1651#if defined(_INDOPRO_)
1652                combine_2_index_domains(res, arg1, arg2, ind_dom);
1653#endif     
1654#if defined(_NONLIND_)
1655                fod[opind].entry = maxopind+2;
1656                fod[opind].left = &fod[arg_index[arg1]];
1657                fod[opind].right = &fod[arg_index[arg2]];
1658                arg_index[res] = opind++;               
1659#endif
1660#else
1661#if !defined(_ZOS_) /* BREAK_ZOS */
1662                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
1663                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
1664                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
1665
1666#ifdef _INT_FOR_
1667                FOR_0_LE_l_LT_pk
1668                TRES_INC = TARG1_INC | TARG2_INC;
1669#else
1670                FOR_0_LE_l_LT_pk
1671                TRES_INC = TARG1_INC + TARG2_INC;
1672#endif
1673#endif
1674#endif /* ALL_TOGETHER_AGAIN */
1675                break;
1676
1677                /*--------------------------------------------------------------------------*/
1678            case plus_d_a:             /* Add an adouble and a double    plus_d_a */
1679                /* (+) */
1680                arg   = get_locint_f();
1681                res   = get_locint_f();
1682                coval = get_val_f();
1683
1684                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1685
1686#if !defined(_NTIGHT_)
1687                dp_T0[res] = dp_T0[arg] + coval;
1688#endif /* !_NTIGHT_ */
1689
1690#if defined(_INDO_)
1691#if defined(_INDOPRO_)
1692                copy_index_domain(res, arg, ind_dom);
1693#endif
1694#if defined(_NONLIND_)
1695                arg_index[res] = arg_index[arg];
1696#endif               
1697#else
1698#if !defined(_ZOS_) /* BREAK_ZOS */
1699                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1700                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
1701
1702                FOR_0_LE_l_LT_pk
1703                TRES_INC = TARG_INC;
1704#endif
1705#endif /* ALL_TOGETHER_AGAIN */
1706                break;
1707
1708                /*--------------------------------------------------------------------------*/
1709            case min_a_a:              /* Subtraction of two adoubles     min_a_a */
1710                /* (-) */
1711                arg1 = get_locint_f();
1712                arg2 = get_locint_f();
1713                res  = get_locint_f();
1714
1715                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1716
1717#if !defined(_NTIGHT_)
1718                dp_T0[res] = dp_T0[arg1] -
1719                                               dp_T0[arg2];
1720#endif /* !_NTIGHT_ */
1721
1722
1723#if defined(_INDO_)   
1724#if defined(_INDOPRO_)
1725                combine_2_index_domains(res, arg1, arg2, ind_dom);
1726#endif
1727#if defined(_NONLIND_)
1728                fod[opind].entry = maxopind+2;
1729                fod[opind].left = &fod[arg_index[arg1]];
1730                fod[opind].right = &fod[arg_index[arg2]];
1731                arg_index[res] = opind++;               
1732#endif
1733#else
1734#if !defined(_ZOS_) /* BREAK_ZOS */
1735                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
1736                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
1737                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
1738
1739#ifdef _INT_FOR_
1740                FOR_0_LE_l_LT_pk
1741                TRES_INC = TARG1_INC | TARG2_INC;
1742#else
1743                 FOR_0_LE_l_LT_pk
1744                TRES_INC = TARG1_INC - TARG2_INC;
1745#endif
1746#endif
1747#endif /* ALL_TOGETHER_AGAIN */
1748                break;
1749
1750                /*--------------------------------------------------------------------------*/
1751            case min_d_a:                /* Subtract an adouble from a    min_d_a */
1752                /* double (-) */
1753                arg =get_locint_f();
1754                res = get_locint_f();
1755                coval = get_val_f();
1756
1757                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1758
1759#if !defined(_NTIGHT_)
1760                dp_T0[res] = coval - dp_T0[arg];
1761#endif /* !_NTIGHT_ */
1762
1763#if defined(_INDO_)
1764#if defined(_INDOPRO_)
1765                copy_index_domain(res, arg, ind_dom);
1766#endif
1767#if defined(_NONLIND_)
1768                arg_index[res] = arg_index[arg];               
1769#endif
1770#else
1771#if !defined(_ZOS_) /* BREAK_ZOS */
1772                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
1773                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
1774
1775#ifdef _INT_FOR_
1776                FOR_0_LE_l_LT_pk
1777                TRES_INC = TARG_INC;
1778#else
1779                FOR_0_LE_l_LT_pk
1780                TRES_INC = -TARG_INC;
1781#endif
1782#endif
1783#endif /* ALL_TOGETHER_AGAIN */
1784                break;
1785
1786                /*--------------------------------------------------------------------------*/
1787            case mult_a_a:               /* Multiply two adoubles (*)    mult_a_a */
1788                arg1 = get_locint_f();
1789                arg2 = get_locint_f();
1790                res  = get_locint_f();
1791
1792                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1793
1794#if defined(_INDO_)
1795#if defined(_INDOPRO_)
1796                combine_2_index_domains(res, arg1, arg2, ind_dom);
1797#endif
1798#if defined(_NONLIND_)
1799                fod[opind].entry = maxopind+2;
1800                fod[opind].left = &fod[arg_index[arg1]];
1801                fod[opind].right = &fod[arg_index[arg2]];
1802                traverse_unary(&fod[arg_index[arg1]], nonl_dom, &fod[arg_index[arg2]], indcheck+1,maxopind+2);
1803                traverse_unary(&fod[arg_index[arg2]], nonl_dom, &fod[arg_index[arg1]], indcheck+1,maxopind+2);
1804                arg_index[res] = opind++;               
1805#endif
1806#else
1807#if !defined(_ZOS_) /* BREAK_ZOS */
1808                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
1809                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
1810                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
1811
1812#ifdef _INT_FOR_
1813                FOR_0_LE_l_LT_p
1814                TRES_FOINC = TARG2_INC | TARG1_INC;
1815#else
1816                /* olvo 980915 now in reverse order to allow x = x*x etc. */
1817                INC_pk_1(Tres)
1818                INC_pk_1(Targ1)
1819                INC_pk_1(Targ2)
1820
1821                FOR_p_GT_l_GE_0
1822                FOR_k_GT_i_GE_0
1823                { TRES_FODEC = dp_T0[arg1]*TARG2_DEC +
1824                               TARG1_DEC*dp_T0[arg2];
1825                  DEC_TRES_FO
1826#if defined(_HIGHER_ORDER_)
1827                  Targ1OP = Targ1-i+1;
1828                  Targ2OP = Targ2;
1829
1830                  for (j=0;j<i;j++) {
1831                  *Tres += (*Targ1OP++) * (*Targ2OP--);
1832                  }
1833                  Tres--;
1834#endif /* _HIGHER_ORDER_ */
1835            }
1836#endif
1837#endif
1838#endif /* ALL_TOGETHER_AGAIN */
1839#if !defined(_NTIGHT_)
1840                dp_T0[res] = dp_T0[arg1] *
1841                                               dp_T0[arg2];
1842#endif /* !_NTIGHT_ */
1843                break;
1844
1845                /*--------------------------------------------------------------------------*/
1846                /* olvo 991122: new op_code with recomputation */
1847            case eq_plus_prod:   /* increment a product of           eq_plus_prod */
1848                /* two adoubles (*) */
1849                arg1 = get_locint_f();
1850                arg2 = get_locint_f();
1851                res  = get_locint_f();
1852
1853#if defined(_INDO_)
1854#if defined(_INDOPRO_)
1855                merge_3_index_domains(res, arg1, arg2, ind_dom);
1856#endif
1857#if defined(_NONLIND_)
1858                // operation: v = v+u*w
1859                // first step: z = u*w, index domains
1860                fod[opind].entry = maxopind+2;
1861                fod[opind].left = &fod[arg_index[arg1]];
1862                fod[opind].right = &fod[arg_index[arg2]];
1863                // first step: z = u*w,
1864                traverse_unary(&fod[arg_index[arg1]], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
1865                traverse_unary(&fod[arg_index[arg2]], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
1866                opind++;
1867                // second step: v = v+z, index domains
1868                fod[opind].entry = maxopind+2;
1869                fod[opind].left = &fod[arg_index[res]];
1870                fod[opind].right = &fod[opind-1];
1871                // second step: v = v+z,
1872                arg_index[res] = opind++;               
1873#endif
1874#else
1875#if !defined(_ZOS_) /* BREAK_ZOS */
1876                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
1877                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
1878                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
1879
1880#ifdef _INT_FOR_
1881                FOR_0_LE_l_LT_p
1882                TRES_FOINC |= TARG2_INC | TARG1_INC;
1883#else
1884                /* olvo 980915 now in reverse order to allow x = x*x etc. */
1885                INC_pk_1(Tres)
1886                INC_pk_1(Targ1)
1887                INC_pk_1(Targ2)
1888
1889                FOR_p_GT_l_GE_0
1890                FOR_k_GT_i_GE_0
1891                { TRES_FODEC += dp_T0[arg1]*TARG2_DEC +
1892                                TARG1_DEC*dp_T0[arg2];
1893                  DEC_TRES_FO
1894#if defined(_HIGHER_ORDER_)
1895                  Targ1OP = Targ1-i+1;
1896                  Targ2OP = Targ2;
1897
1898                  for (j=0;j<i;j++)
1899                  *Tres += (*Targ1OP++) * (*Targ2OP--);
1900                  Tres--;
1901#endif /* _HIGHER_ORDER_ */
1902                }
1903#endif
1904#endif
1905#endif /* ALL_TOGETHER_AGAIN */
1906#if !defined(_NTIGHT_)
1907                dp_T0[res] += dp_T0[arg1] *
1908                                                    dp_T0[arg2];
1909#endif /* !_NTIGHT_ */
1910                break;
1911
1912                /*--------------------------------------------------------------------------*/
1913                /* olvo 991122: new op_code with recomputation */
1914            case eq_min_prod:    /* decrement a product of            eq_min_prod */
1915                /* two adoubles (*) */
1916                arg1 = get_locint_f();
1917                arg2 = get_locint_f();
1918                res  = get_locint_f();
1919
1920#if defined(_INDO_)
1921#if defined(_INDOPRO_)
1922                merge_3_index_domains(res, arg1, arg2, ind_dom);
1923#endif
1924#if defined(_NONLIND_)
1925                // operation: v = v-u*w
1926                // first step: z = u*w, index domains
1927                fod[opind].entry = maxopind+2;
1928                fod[opind].left = &fod[arg_index[arg1]];
1929                fod[opind].right = &fod[arg_index[arg2]];
1930                // first step: z = u*w,
1931                traverse_unary(&fod[arg_index[arg1]], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
1932                traverse_unary(&fod[arg_index[arg2]], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
1933                opind++;
1934                // second step: v = v-z, index domains
1935                fod[opind].entry = maxopind+2;
1936                fod[opind].left = &fod[arg_index[res]];
1937                fod[opind].right = &fod[opind-1];
1938                // second step: v = v-z,
1939                arg_index[res] = opind++;       
1940#endif
1941#else
1942#if !defined(_ZOS_) /* BREAK_ZOS */
1943                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
1944                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
1945                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
1946
1947#ifdef _INT_FOR_
1948                FOR_0_LE_l_LT_p
1949                TRES_FOINC |= TARG2_INC | TARG1_INC;
1950#else
1951                /* olvo 980915 now in reverse order to allow x = x*x etc. */
1952                INC_pk_1(Tres)
1953                INC_pk_1(Targ1)
1954                INC_pk_1(Targ2)
1955
1956                FOR_p_GT_l_GE_0
1957                FOR_k_GT_i_GE_0
1958                { TRES_FODEC -= dp_T0[arg1]*TARG2_DEC +
1959                                TARG1_DEC*dp_T0[arg2];
1960                  DEC_TRES_FO
1961#if defined(_HIGHER_ORDER_)
1962                  Targ1OP = Targ1-i+1;
1963                  Targ2OP = Targ2;
1964
1965                  for (j=0;j<i;j++)
1966                  *Tres -= (*Targ1OP++) * (*Targ2OP--);
1967                  Tres--;
1968#endif /* _HIGHER_ORDER_ */
1969                }
1970#endif
1971#endif
1972#endif /* ALL_TOGETHER_AGAIN */
1973
1974#if !defined(_NTIGHT_)
1975                dp_T0[res] -= dp_T0[arg1] *
1976                                                    dp_T0[arg2];
1977#endif /* !_NTIGHT_ */
1978                break;
1979
1980                /*--------------------------------------------------------------------------*/
1981            case mult_d_a:         /* Multiply an adouble by a double    mult_d_a */
1982                /* (*) */
1983                arg   = get_locint_f();
1984                res   = get_locint_f();
1985                coval = get_val_f();
1986
1987                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
1988
1989#if !defined(_NTIGHT_)
1990                dp_T0[res] = dp_T0[arg] * coval;
1991#endif /* !_NTIGHT_ */
1992
1993#if defined(_INDO_)
1994#if defined(_INDOPRO_)
1995                copy_index_domain(res, arg, ind_dom);
1996#endif
1997#if defined(_NONLIND_)
1998                arg_index[res] = arg_index[arg];               
1999#endif
2000#else
2001#if !defined(_ZOS_) /* BREAK_ZOS */
2002                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
2003                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
2004
2005#ifdef _INT_FOR_
2006                FOR_0_LE_l_LT_pk
2007                TRES_INC = TARG_INC;
2008#else
2009                FOR_0_LE_l_LT_pk
2010                TRES_INC = TARG_INC * coval;
2011#endif
2012#endif
2013#endif /* ALL_TOGETHER_AGAIN */
2014                break;
2015
2016                /*--------------------------------------------------------------------------*/
2017            case div_a_a:           /* Divide an adouble by an adouble    div_a_a */
2018                /* (/) */
2019                arg1 = get_locint_f();
2020                arg2 = get_locint_f();
2021                res  = get_locint_f();
2022
2023                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2024
2025#if !defined(_NTIGHT_)
2026#if !defined(_ZOS_) /* BREAK_ZOS */
2027                divs = 1.0 / dp_T0[arg2];
2028#endif /* ALL_TOGETHER_AGAIN */
2029
2030                dp_T0[res] = dp_T0[arg1] /
2031                                               dp_T0[arg2];
2032#endif /* !_NTIGHT_ */
2033
2034#if defined(_INDO_)
2035#if defined(_INDOPRO_)
2036                combine_2_index_domains(res, arg1, arg2, ind_dom);
2037#endif
2038#if defined(_NONLIND_)
2039                fod[opind].entry = maxopind+2;
2040                fod[opind].left = &fod[arg_index[arg1]];
2041                fod[opind].right = &fod[arg_index[arg2]];
2042                traverse_unary(&fod[arg_index[arg1]], nonl_dom, &fod[arg_index[arg2]], indcheck+1,maxopind+2);
2043                traverse_unary(&fod[arg_index[arg2]], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2044                arg_index[res] = opind++;               
2045#endif
2046#else
2047#if !defined(_ZOS_) /* BREAK_ZOS */
2048                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2049                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2050                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2051
2052#ifdef _INT_FOR_
2053                FOR_0_LE_l_LT_p
2054                TRES_FOINC = TARG1_INC | TARG2_FOINC;
2055#else
2056                FOR_0_LE_l_LT_p
2057                FOR_0_LE_i_LT_k
2058                { /* olvo 980922 changed order to allow x = y/x */
2059#if defined(_HIGHER_ORDER_)
2060                    zOP      = dp_z+i;
2061                    (*zOP--) = -(*Targ2) * divs;
2062#endif /* _HIGHER_ORDER_ */
2063
2064                    TRES_FOINC = TARG1_INC * divs + dp_T0[res] *
2065                                 (-TARG2_INC * divs);
2066
2067#if defined(_HIGHER_ORDER_)
2068                    TresOP = Tres-i;
2069
2070                    for (j=0;j<i;j++)
2071                    *Tres += (*TresOP++) * (*zOP--);
2072                    Tres++;
2073#endif /* _HIGHER_ORDER_ */
2074                }
2075#endif
2076#endif
2077#endif /* ALL_TOGETHER_AGAIN */
2078                break;
2079
2080            /*--------------------------------------------------------------------------*/
2081        case div_d_a:             /* Division double - adouble (/)    div_d_a */
2082            arg   = get_locint_f();
2083                res   = get_locint_f();
2084                coval = get_val_f();
2085
2086                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2087
2088                /* olvo 980922 necessary for reverse */
2089                if (arg == res) {
2090                    IF_KEEP_WRITE_TAYLOR(arg,keep,k,p)
2091                }
2092
2093#if !defined(_NTIGHT_)
2094#if !defined(_ZOS_) /* BREAK_ZOS */
2095                divs = 1.0 / dp_T0[arg];
2096#endif /* ALL_TOGETHER_AGAIN */
2097
2098                dp_T0[res] = coval / dp_T0[arg];
2099#endif /* !_NTIGHT_ */
2100
2101#if defined(_INDO_)
2102#if defined(_INDOPRO_)
2103                copy_index_domain(res, arg, ind_dom);
2104#endif
2105#if defined(_NONLIND_)
2106                fod[opind].entry = maxopind+2;
2107                fod[opind].left = &fod[arg_index[arg]];
2108                fod[opind].right = NULL;
2109                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2110                arg_index[res] = opind++;               
2111#endif
2112#else
2113#if !defined(_ZOS_) /* BREAK_ZOS */
2114                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
2115                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
2116
2117#ifdef _INT_FOR_
2118                FOR_0_LE_l_LT_p
2119                TRES_FOINC = TARG_FOINC;
2120#else
2121                FOR_0_LE_l_LT_p
2122                FOR_0_LE_i_LT_k
2123                { /* olvo 980922 changed order to allow x = d/x */
2124#if defined(_HIGHER_ORDER_)
2125                    zOP      = dp_z+i;
2126                    (*zOP--) = -(*Targ) * divs;
2127#endif /* _HIGHER_ORDER_ */
2128
2129                    TRES_FOINC = dp_T0[res] * (-TARG_INC * divs);
2130
2131#if defined(_HIGHER_ORDER_)
2132                    TresOP = Tres-i;
2133
2134                    for (j=0;j<i;j++)
2135                    *Tres += (*TresOP++) * (*zOP--);
2136                    Tres++;
2137#endif /* _HIGHER_ORDER_ */
2138                }
2139#endif
2140#endif
2141#endif /* ALL_TOGETHER_AGAIN */
2142                break;
2143
2144
2145            /****************************************************************************/
2146            /*                                                         SIGN  OPERATIONS */
2147
2148            /*--------------------------------------------------------------------------*/
2149        case pos_sign_a:                                        /* pos_sign_a */
2150            arg   = get_locint_f();
2151                res   = get_locint_f();
2152
2153                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2154
2155#if !defined(_NTIGHT_)
2156                dp_T0[res] = dp_T0[arg];
2157#endif /* !_NTIGHT_ */
2158
2159#if defined(_INDO_)
2160#if defined(_INDOPRO_)
2161                copy_index_domain(res, arg, ind_dom);
2162#endif
2163#if defined(_NONLIND_)
2164                arg_index[res] = arg_index[arg];                           
2165#endif 
2166#else
2167#if !defined(_ZOS_) /* BREAK_ZOS */
2168                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
2169                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
2170
2171                FOR_0_LE_l_LT_pk
2172                TRES_INC = TARG_INC;
2173#endif
2174#endif /* ALL_TOGETHER_AGAIN */
2175                break;
2176
2177                /*--------------------------------------------------------------------------*/
2178            case neg_sign_a:                                        /* neg_sign_a */
2179                arg   = get_locint_f();
2180                res   = get_locint_f();
2181
2182                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2183
2184#if !defined(_NTIGHT_)
2185                dp_T0[res] = -dp_T0[arg];
2186#endif /* !_NTIGHT_ */
2187
2188#if defined(_INDO_)
2189#if defined(_INDOPRO_)
2190                copy_index_domain(res, arg, ind_dom);
2191#endif
2192#if defined(_NONLIND_)
2193                arg_index[res] = arg_index[arg];                           
2194#endif
2195#else
2196#if !defined(_ZOS_) /* BREAK_ZOS */
2197                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
2198                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
2199
2200#ifdef _INT_FOR_
2201                FOR_0_LE_l_LT_pk
2202                TRES_INC = TARG_INC;
2203#else
2204                FOR_0_LE_l_LT_pk
2205                TRES_INC = -TARG_INC;
2206#endif
2207#endif
2208#endif /* ALL_TOGETHER_AGAIN */
2209                break;
2210
2211
2212                /****************************************************************************/
2213                /*                                                         UNARY OPERATIONS */
2214
2215                /*--------------------------------------------------------------------------*/
2216            case exp_op:                          /* exponent operation    exp_op */
2217                arg = get_locint_f();
2218                res = get_locint_f();
2219
2220                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2221
2222#if !defined(_NTIGHT_)
2223                dp_T0[res] = exp(dp_T0[arg]);
2224#endif /* !_NTIGHT_ */
2225
2226                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2227
2228#if defined(_INDO_)
2229#if defined(_INDOPRO_)
2230                copy_index_domain(res, arg, ind_dom);
2231#endif
2232#if defined(_NONLIND_)
2233                fod[opind].entry = maxopind+2;
2234                fod[opind].left = &fod[arg_index[arg]];
2235                fod[opind].right = NULL;
2236                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2237                arg_index[res] = opind++;               
2238#endif
2239#else
2240#if !defined(_ZOS_) /* BREAK_ZOS */
2241                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
2242                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
2243
2244#ifdef _INT_FOR_
2245                FOR_0_LE_l_LT_p
2246                TRES_FOINC = TARG_FOINC;
2247#else
2248                FOR_0_LE_l_LT_p
2249                FOR_0_LE_i_LT_k
2250                { /* olvo 980915 changed order to allow x = exp(x) */
2251#if defined(_HIGHER_ORDER_)
2252                    zOP      = dp_z+i;
2253                    (*zOP--) = (i+1) * (*Targ);
2254#endif /* _HIGHER_ORDER_ */
2255
2256                    TRES_FOINC = dp_T0[res] * TARG_INC;
2257
2258#if defined(_HIGHER_ORDER_)
2259                    TresOP = Tres-i;
2260
2261                    *Tres *= (i+1);
2262                    for (j=0;j<i;j++)
2263                    *Tres += (*TresOP++) * (*zOP--);
2264                    *Tres++ /= (i+1); /* important only for i>0 */
2265#endif /* _HIGHER_ORDER_ */
2266                }
2267
2268#endif
2269#endif
2270#endif /* ALL_TOGETHER_AGAIN */
2271                break;
2272
2273            /*--------------------------------------------------------------------------*/
2274        case sin_op:                              /* sine operation    sin_op */
2275                arg1 = get_locint_f();
2276                arg2 = get_locint_f();
2277                res  = get_locint_f();
2278
2279                IF_KEEP_WRITE_TAYLOR(arg2,keep,k,p) /* olvo 980710 covalue */
2280                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2281
2282#if !defined(_NTIGHT_)
2283                /* Note: always arg2 != arg1 */
2284                dp_T0[arg2] = cos(dp_T0[arg1]);
2285                dp_T0[res]  = sin(dp_T0[arg1]);
2286#endif /* !_NTIGHT_ */
2287
2288                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2289
2290#if defined(_INDO_)
2291#if defined(_INDOPRO_)
2292                copy_index_domain(res, arg1, ind_dom);
2293#endif
2294#if defined(_NONLIND_)
2295                fod[opind].entry = maxopind+2;
2296                fod[opind].left = &fod[arg_index[arg1]];
2297                fod[opind].right = NULL;
2298                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2299                arg_index[res] = opind++;               
2300#endif
2301#else
2302#if !defined(_ZOS_) /* BREAK_ZOS */
2303                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2304                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2305                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2306
2307#ifdef _INT_FOR_
2308                FOR_0_LE_l_LT_p
2309                { /* olvo 980923 changed order to allow x = sin(x) */
2310                    TARG2_FOINC =  TARG1;
2311                    TRES_FOINC  =  TARG1_FOINC;
2312            }
2313#else
2314                FOR_0_LE_l_LT_p
2315                FOR_0_LE_i_LT_k
2316                { /* olvo 980921 changed order to allow x = sin(x) */
2317#if defined(_HIGHER_ORDER_)
2318                    zOP      = dp_z+i;
2319                    (*zOP--) = (i+1) * (*Targ1);
2320#endif /* _HIGHER_ORDER_ */
2321
2322                    /* Note: always arg2 != arg1 */
2323                    TARG2_FOINC = -dp_T0[res]  * TARG1;
2324                    TRES_FOINC  =  dp_T0[arg2] * TARG1_INC;
2325
2326#if defined(_HIGHER_ORDER_)
2327                    TresOP  = Tres-i;
2328                    Targ2OP = Targ2-i;
2329
2330                    *Tres  *= (i+1);
2331                    *Targ2 *= (i+1);
2332                    for (j=0;j<i;j++) {
2333                    *Tres  += (*Targ2OP++) * (*zOP);
2334                        *Targ2 -= (*TresOP++)  * (*zOP--);
2335                    }
2336                    *Targ2++ /= (i+1);
2337                    *Tres++  /= (i+1);
2338#endif /* _HIGHER_ORDER_ */
2339            }
2340#endif
2341#endif
2342#endif /* ALL_TOGETHER_AGAIN */
2343                break;
2344
2345                /*--------------------------------------------------------------------------*/
2346            case cos_op:                            /* cosine operation    cos_op */
2347                arg1 = get_locint_f();
2348                arg2 = get_locint_f();
2349                res  = get_locint_f();
2350
2351                IF_KEEP_WRITE_TAYLOR(arg2,keep,k,p) /* olvo 980710 covalue */
2352                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2353
2354#if !defined(_NTIGHT_)
2355                /* Note: always arg2 != arg1 */
2356                dp_T0[arg2] = sin(dp_T0[arg1]);
2357                dp_T0[res]  = cos(dp_T0[arg1]);
2358#endif /* !_NTIGHT_ */
2359
2360                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2361
2362#if defined(_INDO_)
2363#if defined(_INDOPRO_)
2364                copy_index_domain(res, arg1, ind_dom);
2365#endif
2366#if defined(_NONLIND_)
2367                fod[opind].entry = maxopind+2;
2368                fod[opind].left = &fod[arg_index[arg1]];
2369                fod[opind].right = NULL;
2370                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2371                arg_index[res] = opind++;               
2372#endif
2373
2374#else
2375#if !defined(_ZOS_) /* BREAK_ZOS */
2376                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2377                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2378                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2379
2380#ifdef _INT_FOR_
2381                FOR_0_LE_l_LT_p
2382                { /* olvo 980923 changed order to allow x = cos(x) */
2383                    TARG2_FOINC = TARG1;
2384                    TRES_FOINC  = TARG1_FOINC;
2385            }
2386#else
2387                FOR_0_LE_l_LT_p
2388                FOR_0_LE_i_LT_k
2389                { /* olvo 980921 changed order to allow x = cos(x) */
2390#if defined(_HIGHER_ORDER_)
2391                    zOP      = dp_z+i;
2392                    (*zOP--) = (i+1) * (*Targ1);
2393#endif /* _HIGHER_ORDER_ */
2394
2395                    /* Note: always arg2 != arg1 */
2396                    TARG2_FOINC =  dp_T0[res]  * TARG1;
2397                    TRES_FOINC  = -dp_T0[arg2] * TARG1_INC;
2398
2399#if defined(_HIGHER_ORDER_)
2400                    TresOP  = Tres-i;
2401                    Targ2OP = Targ2-i;
2402
2403                    *Tres  *= (i+1);
2404                    *Targ2 *= (i+1);
2405                    for (j=0;j<i;j++) {
2406                    *Tres  -= (*Targ2OP++) * (*zOP);
2407                        *Targ2 += (*TresOP++)  * (*zOP--);
2408                    }
2409                    *Targ2++ /= (i+1);
2410                    *Tres++  /= (i+1);
2411#endif /* _HIGHER_ORDER_ */
2412            }
2413#endif
2414#endif
2415#endif /* ALL_TOGETHER_AGAIN */
2416                break;
2417
2418                /*--------------------------------------------------------------------------*/
2419            case atan_op:                                              /* atan_op */
2420                arg1 = get_locint_f();
2421                arg2 = get_locint_f();
2422                res  = get_locint_f();
2423
2424                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2425
2426#if !defined(_NTIGHT_)
2427                dp_T0[res]=atan(dp_T0[arg1]);
2428#endif /* !_NTIGHT_ */
2429
2430                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2431
2432#if defined(_INDO_)
2433#if defined(_INDOPRO_)
2434                copy_index_domain(res, arg1, ind_dom);
2435#endif
2436#if defined(_NONLIND_)
2437                fod[opind].entry = maxopind+2;
2438                fod[opind].left = &fod[arg_index[arg1]];
2439                fod[opind].right = NULL;
2440                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2441                arg_index[res] = opind++;               
2442#endif
2443#else
2444#if !defined(_ZOS_) /* BREAK_ZOS */
2445                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2446                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2447#ifdef _INT_FOR_
2448                FOR_0_LE_l_LT_p
2449                TRES_FOINC = TARG1_FOINC;
2450#else
2451                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2452
2453                FOR_0_LE_l_LT_p
2454                { FOR_0_LE_i_LT_k
2455                  { /* olvo 980921 changed order to allow x = atan(x) */
2456#if defined(_HIGHER_ORDER_)
2457                      zOP      = dp_z+i;
2458                      (*zOP--) = (i+1) * (*Targ1);
2459#endif /* _HIGHER_ORDER_ */
2460
2461                      TRES_FOINC = dp_T0[arg2] * TARG1_INC;
2462
2463#if defined(_HIGHER_ORDER_)
2464                      Targ2OP = Targ2;
2465
2466                      *Tres *= (i+1);
2467                      for (j=0;j<i;j++)
2468                      *Tres  += (*Targ2OP++) * (*zOP--);
2469                      *Tres++ /= (i+1);
2470#endif /* _HIGHER_ORDER_ */
2471                  }
2472                  HOV_INC(Targ2, k)
2473                }
2474#endif
2475#endif
2476#endif /* ALL_TOGETHER_AGAIN */
2477                break;
2478
2479            /*--------------------------------------------------------------------------*/
2480        case asin_op:                                              /* asin_op */
2481            arg1 = get_locint_f();
2482                arg2 = get_locint_f();
2483                res  = get_locint_f();
2484
2485                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2486
2487#if !defined(_NTIGHT_)
2488                dp_T0[res] = asin(dp_T0[arg1]);
2489#endif /* !_NTIGHT_ */
2490
2491                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2492
2493#if defined(_INDO_)
2494#if defined(_INDOPRO_)
2495                copy_index_domain(res, arg1, ind_dom);
2496#endif
2497#if defined(_NONLIND_)
2498                fod[opind].entry = maxopind+2;
2499                fod[opind].left = &fod[arg_index[arg1]];
2500                fod[opind].right = NULL;
2501                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2502                arg_index[res] = opind++;               
2503#endif
2504#else
2505#if !defined(_ZOS_) /* BREAK_ZOS */
2506                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2507                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2508#ifdef _INT_FOR_
2509                FOR_0_LE_l_LT_p
2510                TRES_FOINC = TARG1_FOINC;
2511#else
2512                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2513
2514                if (dp_T0[arg1] == 1.0)
2515                    FOR_0_LE_l_LT_p
2516                    { FOR_0_LE_i_LT_k
2517                      if (TARG1 > 0.0) {
2518                      r0 = make_nan();
2519                          VEC_INC(Targ1, k-i)
2520                          BREAK_FOR_I
2521                      } else
2522                          if (TARG1 < 0.0) {
2523                          r0 = make_inf();
2524                              VEC_INC(Targ1, k-i)
2525                              BREAK_FOR_I
2526                          } else {
2527                              r0 = 0.0;
2528                              Targ1++;
2529                          }
2530                  TRES = r0;
2531                  VEC_INC(Tres, k)
2532            } else
2533                    if (dp_T0[arg1] == -1.0)
2534                        FOR_0_LE_l_LT_p
2535                        { FOR_0_LE_i_LT_k
2536                          if (TARG1 > 0.0) {
2537                          r0 = make_inf();
2538                              VEC_INC(Targ1, k-i)
2539                              BREAK_FOR_I
2540                          } else
2541                              if (TARG1 < 0.0) {
2542                              r0 = make_nan();
2543                                  VEC_INC(Targ1, k-i)
2544                                  BREAK_FOR_I
2545                              } else {
2546                                  r0 = 0.0;
2547                                  Targ1++;
2548                              }
2549                  TRES = r0;
2550                  VEC_INC(Tres, k)
2551                } else
2552                        FOR_0_LE_l_LT_p {
2553                            FOR_0_LE_i_LT_k
2554                            { /* olvo 980921 changed order to allow x = asin(x) */
2555#if defined(_HIGHER_ORDER_)
2556                                zOP      = dp_z+i;
2557                                (*zOP--) = (i+1) * (*Targ1);
2558#endif /* _HIGHER_ORDER_ */
2559
2560                                TRES_FOINC = dp_T0[arg2] * TARG1_INC;
2561
2562#if defined(_HIGHER_ORDER_)
2563                                Targ2OP = Targ2;
2564
2565                                *Tres *= (i+1);
2566                                for (j=0;j<i;j++)
2567                                *Tres += (*Targ2OP++) * (*zOP--);
2568                                *Tres++ /= (i+1);
2569#endif /* _HIGHER_ORDER_ */
2570                            }
2571                            HOV_INC(Targ2, k)
2572                        }
2573#endif
2574#endif
2575#endif /* ALL_TOGETHER_AGAIN */
2576                        break;
2577
2578            /*--------------------------------------------------------------------------*/
2579        case acos_op:                                              /* acos_op */
2580                arg1 = get_locint_f();
2581                arg2 = get_locint_f();
2582                res  = get_locint_f();
2583
2584                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2585
2586#if !defined(_NTIGHT_)
2587                dp_T0[res] = acos(dp_T0[arg1]);
2588#endif /* !_NTIGHT_ */
2589
2590                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2591
2592#if defined(_INDO_)
2593#if defined(_INDOPRO_)
2594                copy_index_domain(res, arg1, ind_dom);
2595#endif
2596#if defined(_NONLIND_)
2597                fod[opind].entry = maxopind+2;
2598                fod[opind].left = &fod[arg_index[arg1]];
2599                fod[opind].right = NULL;
2600                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2601                arg_index[res] = opind++;               
2602#endif
2603#else
2604#if !defined(_ZOS_) /* BREAK_ZOS */
2605                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2606                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2607#ifdef _INT_FOR_
2608                FOR_0_LE_l_LT_p
2609                TRES_FOINC = TARG1_FOINC;
2610#else
2611                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2612
2613                if (dp_T0[arg1] == 1.0)
2614                    FOR_0_LE_l_LT_p
2615                    { FOR_0_LE_i_LT_k
2616                      if (TARG1 > 0.0) {
2617                      r0 = make_nan();
2618                          VEC_INC(Targ1, k-i)
2619                          BREAK_FOR_I
2620                      } else
2621                          if (TARG1 < 0.0) {
2622                          r0 = -make_inf();
2623                              VEC_INC(Targ1, k-i)
2624                              BREAK_FOR_I
2625                          } else {
2626                              r0 = 0.0;
2627                              Targ1++;
2628                          }
2629                  TRES = r0;
2630                  VEC_INC(Tres, k)
2631            } else
2632                    if (dp_T0[arg1] == -1.0)
2633                        FOR_0_LE_l_LT_p
2634                        { FOR_0_LE_i_LT_k
2635                          if (TARG1 > 0.0) {
2636                          r0 = -make_inf();
2637                              VEC_INC(Targ1, k-i)
2638                              BREAK_FOR_I
2639                          } else
2640                              if (TARG1 < 0.0) {
2641                              r0 = make_nan();
2642                                  VEC_INC(Targ1, k-i)
2643                                  BREAK_FOR_I
2644                              } else {
2645                                  r0 = 0.0;
2646                                  Targ1++;
2647                              }
2648                  TRES = r0;
2649                  VEC_INC(Tres, k)
2650                } else
2651                        FOR_0_LE_l_LT_p {
2652                            FOR_0_LE_i_LT_k
2653                            { /* olvo 980921 changed order to allow x = acos(x) */
2654#if defined(_HIGHER_ORDER_)
2655                                zOP      = dp_z+i;
2656                                (*zOP--) = (i+1) * (*Targ1);
2657#endif /* _HIGHER_ORDER_ */
2658
2659                                TRES_FOINC = dp_T0[arg2] * TARG1_INC;
2660
2661#if defined(_HIGHER_ORDER_)
2662                                Targ2OP = Targ2;
2663
2664                                *Tres *= (i+1);
2665                                for (j=0;j<i;j++)
2666                                *Tres += (*Targ2OP++) * (*zOP--);
2667                                *Tres++ /= (i+1);
2668#endif /* _HIGHER_ORDER_ */
2669                            }
2670                            HOV_INC(Targ2, k)
2671                        }
2672#endif
2673#endif
2674#endif /* ALL_TOGETHER_AGAIN */
2675                        break;
2676
2677#ifdef ATRIG_ERF
2678
2679            /*--------------------------------------------------------------------------*/
2680        case asinh_op:                                            /* asinh_op */
2681                arg1 = get_locint_f();
2682                arg2 = get_locint_f();
2683                res  = get_locint_f();
2684
2685                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2686
2687#if !defined(_NTIGHT_)
2688                dp_T0[res] = asinh(dp_T0[arg1]);
2689#endif /* !_NTIGHT_ */
2690
2691                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2692
2693#if defined(_INDO_)
2694#if defined(_INDOPRO_)
2695                copy_index_domain(res, arg1, ind_dom);
2696#endif
2697#if defined(_NONLIND_)
2698                fod[opind].entry = maxopind+2;
2699                fod[opind].left = &fod[arg_index[arg1]];
2700                fod[opind].right = NULL;
2701                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2702                arg_index[res] = opind++;               
2703#endif
2704#else
2705#if !defined(_ZOS_) /* BREAK_ZOS */
2706                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2707                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2708#ifdef _INT_FOR_
2709                FOR_0_LE_l_LT_p
2710                TRES_FOINC = TARG1_FOINC;
2711#else
2712                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2713
2714                FOR_0_LE_l_LT_p
2715                { FOR_0_LE_i_LT_k
2716                  { /* olvo 980921 changed order to allow x = asinh(x) */
2717#if defined(_HIGHER_ORDER_)
2718                      zOP      = dp_z+i;
2719                      (*zOP--) = (i+1) * (*Targ1);
2720#endif /* _HIGHER_ORDER_ */
2721
2722                      TRES_FOINC = dp_T0[arg2] * TARG1_INC;
2723
2724#if defined(_HIGHER_ORDER_)
2725                      Targ2OP = Targ2;
2726
2727                      *Tres *= (i+1);
2728                      for (j=0;j<i;j++)
2729                      *Tres += (*Targ2OP++) * (*zOP--);
2730                      *Tres++ /= (i+1);
2731#endif /* _HIGHER_ORDER_ */
2732                  }
2733                  HOV_INC(Targ2, k)
2734                }
2735#endif
2736#endif
2737#endif /* ALL_TOGETHER_AGAIN */
2738                break;
2739
2740            /*--------------------------------------------------------------------------*/
2741        case acosh_op:                                           /* acosh_op */
2742                arg1 = get_locint_f();
2743                arg2 = get_locint_f();
2744                res  = get_locint_f();
2745
2746                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2747
2748#if !defined(_NTIGHT_)
2749                dp_T0[res] = acosh(dp_T0[arg1]);
2750#endif /* !_NTIGHT_ */
2751
2752                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2753
2754#if defined(_INDO_)
2755#if defined(_INDOPRO_)
2756                copy_index_domain(res, arg1, ind_dom);
2757#endif
2758#if defined(_NONLIND_)
2759                fod[opind].entry = maxopind+2;
2760                fod[opind].left = &fod[arg_index[arg1]];
2761                fod[opind].right = NULL;
2762                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2763                arg_index[res] = opind++;               
2764#endif
2765#else
2766#if !defined(_ZOS_) /* BREAK_ZOS */
2767                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2768                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2769#ifdef _INT_FOR_
2770                FOR_0_LE_l_LT_p
2771                TRES_FOINC = TARG1_FOINC;
2772#else
2773                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2774
2775                if (dp_T0[arg1] == 1.0)
2776                    FOR_0_LE_l_LT_p
2777                    { FOR_0_LE_i_LT_k
2778                      if (TARG1 > 0.0) {
2779                      r0 = make_inf();
2780                          VEC_INC(Targ1, k-i)
2781                          BREAK_FOR_I
2782                      } else
2783                          if (TARG1 < 0.0) {
2784                          r0 = make_nan();
2785                              VEC_INC(Targ1, k-i)
2786                              BREAK_FOR_I
2787                          } else {
2788                              r0 = 0.0;
2789                              Targ1++;
2790                          }
2791                  TRES_INC = r0;
2792#if defined(_HIGHER_ORDER_)
2793                  for (i=1;i<k;i++)
2794                  *Tres++ = make_nan();
2795#endif /* _HIGHER_ORDER_ */
2796                } else
2797                    FOR_0_LE_l_LT_p {
2798                        FOR_0_LE_i_LT_k
2799                        { /* olvo 980921 changed order to allow x = acosh(x) */
2800#if defined(_HIGHER_ORDER_)
2801                            zOP      = dp_z+i;
2802                            (*zOP--) = (i+1) * (*Targ1);
2803#endif /* _HIGHER_ORDER_ */
2804
2805                            TRES_FOINC = dp_T0[arg2] * TARG1_INC;
2806
2807#if defined(_HIGHER_ORDER_)
2808                            Targ2OP = Targ2;
2809
2810                            *Tres *= (i+1);
2811                            for (j=0;j<i;j++)
2812                                *Tres += (*Targ2OP++) * (*zOP--);
2813                                *Tres++ /= (i+1);
2814#endif /* _HIGHER_ORDER_ */
2815                            }
2816                            HOV_INC(Targ2, k)
2817                        }
2818#endif
2819#endif
2820#endif /* ALL_TOGETHER_AGAIN */
2821                        break;
2822
2823            /*--------------------------------------------------------------------------*/
2824        case atanh_op:                                            /* atanh_op */
2825                arg1 = get_locint_f();
2826                arg2 = get_locint_f();
2827                res  = get_locint_f();
2828
2829                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2830
2831#if !defined(_NTIGHT_)
2832                dp_T0[res] = atanh(dp_T0[arg1]);
2833#endif /* !_NTIGHT_ */
2834
2835                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2836
2837#if defined(_INDO_)
2838#if defined(_INDOPRO_)
2839                copy_index_domain(res, arg1, ind_dom);
2840#endif
2841#if defined(_NONLIND_)
2842                fod[opind].entry = maxopind+2;
2843                fod[opind].left = &fod[arg_index[arg1]];
2844                fod[opind].right = NULL;
2845                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2846                arg_index[res] = opind++;               
2847#endif
2848#else
2849#if !defined(_ZOS_) /* BREAK_ZOS */
2850                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
2851                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
2852#ifdef _INT_FOR_
2853                FOR_0_LE_l_LT_p
2854                TRES_FOINC = TARG1_FOINC;
2855#else
2856                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
2857
2858                if (dp_T0[arg1] == 1.0)
2859                    FOR_0_LE_l_LT_p
2860                    { FOR_0_LE_i_LT_k
2861                      if (TARG1 > 0.0) {
2862                      r0 = make_nan();
2863                          VEC_INC(Targ1, k-i)
2864                          BREAK_FOR_I
2865                      } else
2866                          if (TARG1 < 0.0) {
2867                          r0 = make_inf();
2868                              VEC_INC(Targ1, k-i)
2869                              BREAK_FOR_I
2870                          } else {
2871                              r0 = 0.0;
2872                              Targ1++;
2873                          }
2874                  TRES_INC = r0;
2875#if defined(_HIGHER_ORDER_)
2876                  for (i=1;i<k;i++)
2877                  *Tres++ = make_nan();
2878#endif /* _HIGHER_ORDER_ */
2879                } else
2880                    if (dp_T0[arg1] == -1.0)
2881                            FOR_0_LE_l_LT_p
2882                            { FOR_0_LE_i_LT_k
2883                              if (TARG1 > 0.0) {
2884                              r0 = make_inf();
2885                                  VEC_INC(Targ1, k-i)
2886                                  BREAK_FOR_I
2887                              } else
2888                                  if (TARG1 < 0.0) {
2889                                  r0 = make_nan();
2890                                      VEC_INC(Targ1, k-i)
2891                                      BREAK_FOR_I
2892                                  } else {
2893                                      r0 = 0.0;
2894                                      Targ1++;
2895                                  }
2896                  TRES_INC = r0;
2897#if defined(_HIGHER_ORDER_)
2898                  for (i=1;i<k;i++)
2899                  *Tres++ = make_nan();
2900#endif /* _HIGHER_ORDER_ */
2901                        } else
2902                            FOR_0_LE_l_LT_p {
2903                                FOR_0_LE_i_LT_k
2904                                { /* olvo 980921 changed order to allow x = atanh(x) */
2905#if defined(_HIGHER_ORDER_)
2906                                    zOP      = dp_z+i;
2907                                    (*zOP--) = (i+1) * (*Targ1);
2908#endif /* _HIGHER_ORDER_ */
2909
2910                                    TRES_FOINC = dp_T0[arg2] * TARG1_INC;
2911
2912#if defined(_HIGHER_ORDER_)
2913                                    Targ2OP = Targ2;
2914
2915                                    *Tres *= (i+1);
2916                                    for (j=0;j<i;j++)
2917                                        *Tres += (*Targ2OP++) * (*zOP--);
2918                                        *Tres++ /= (i+1);
2919#endif /* _HIGHER_ORDER_ */
2920                                    }
2921                                    HOV_INC(Targ2, k)
2922                                }
2923#endif
2924#endif
2925#endif /* ALL_TOGETHER_AGAIN */
2926                                break;
2927
2928            /*--------------------------------------------------------------------------*/
2929        case erf_op:                                                /* erf_op */
2930                arg1 = get_locint_f();
2931                arg2 = get_locint_f();
2932                res  = get_locint_f();
2933
2934                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2935
2936#if !defined(_NTIGHT_)
2937                dp_T0[res] = erf(dp_T0[arg1]);
2938#endif /* !_NTIGHT_ */
2939
2940                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
2941#if defined(_INDO_)
2942#if defined(_INDOPRO_)
2943                copy_index_domain(res, arg1, ind_dom);
2944#endif
2945#if defined(_NONLIND_)
2946                fod[opind].entry = maxopind+2;
2947                fod[opind].left = &fod[arg_index[arg1]];
2948                fod[opind].right = NULL;
2949                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
2950                arg_index[res] = opind++;               
2951#endif       
2952#else
2953#if !defined(_ZOS_) /* BREAK_ZOS */
2954                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
2955                ASSIGN_T(Targ1,TAYLOR_BUFFER[arg1])
2956#ifdef _INT_FOR_
2957                FOR_0_LE_l_LT_p
2958                TRES_FOINC = TARG1_FOINC;
2959#else
2960                ASSIGN_T(Targ2,TAYLOR_BUFFER[arg2])
2961
2962                FOR_0_LE_l_LT_p
2963                { FOR_0_LE_i_LT_k
2964                  { /* olvo 980921 changed order to allow x = erf(x) */
2965#if defined(_HIGHER_ORDER_)
2966                      zOP      = dp_z+i;
2967                      (*zOP--) = (i+1) * (*Targ1);
2968#endif /* _HIGHER_ORDER_ */
2969
2970                      TRES_FOINC = dp_T0[arg2] * TARG1_INC;
2971
2972#if defined(_HIGHER_ORDER_)
2973                      Targ2OP = Targ2;
2974
2975                      *Tres *= (i+1);
2976                      for (j=0;j<i;j++)
2977                      *Tres += (*Targ2OP++) * (*zOP--);
2978                      *Tres++ /= (i+1);
2979#endif /* _HIGHER_ORDER_ */
2980                  }
2981                  HOV_INC(Targ2, k)
2982                }
2983#endif
2984#endif
2985#endif /* ALL_TOGETHER_AGAIN */
2986                break;
2987
2988#endif
2989
2990            /*--------------------------------------------------------------------------*/
2991        case log_op:                                                /* log_op */
2992                arg = get_locint_f();
2993                res = get_locint_f();
2994
2995                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
2996
2997#if defined(_INDO_)
2998#if defined(_INDOPRO_)
2999                copy_index_domain(res, arg, ind_dom);
3000#endif
3001#if defined(_NONLIND_)
3002                fod[opind].entry = maxopind+2;
3003                fod[opind].left = &fod[arg_index[arg]];
3004                fod[opind].right = NULL;
3005                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
3006                arg_index[res] = opind++;               
3007#endif
3008#else
3009#if !defined(_ZOS_) /* BREAK_ZOS */
3010                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
3011                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
3012
3013#ifdef _INT_FOR_
3014                FOR_0_LE_l_LT_p
3015                TRES_FOINC = TARG_INC;
3016#else
3017                divs = 1.0 / dp_T0[arg];
3018                FOR_0_LE_l_LT_p
3019                { if (dp_T0[arg] == 0.0) {
3020                  TargOP = Targ;
3021                  FOR_0_LE_i_LT_k
3022                  { if (*TargOP++ < 0.0) {
3023                        divs = make_nan();
3024                            BREAK_FOR_I
3025                        }
3026                      }
3027                  }
3028
3029                  /* olvo 980921 changed order to allow x = log(x) */
3030                  FOR_0_LE_i_LT_k
3031                  { TRES_FOINC = TARG_INC * divs;
3032#if defined(_HIGHER_ORDER_)
3033                    TresOP = Tres - i;
3034                    zOP    = dp_z+i;
3035
3036                    (*zOP--) = *Tres;
3037                    (*Tres) *= i+1;
3038                    for (j=0;j<i;j++)
3039                    (*Tres) -= (*zOP--) * (*TresOP++) * (j+1);
3040                    *Tres++ /= i+1;
3041#endif /* _HIGHER_ORDER_ */
3042                  }
3043                }
3044#endif
3045#endif
3046#endif /* ALL_TOGETHER_AGAIN */
3047#if !defined(_NTIGHT_)
3048                dp_T0[res] = log(dp_T0[arg]);
3049#endif /* !_NTIGHT_ */
3050
3051                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
3052                break;
3053
3054                /*--------------------------------------------------------------------------*/
3055            case pow_op:                                                /* pow_op */
3056                arg   = get_locint_f();
3057                res   = get_locint_f();
3058
3059                coval = get_val_f();
3060
3061                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3062
3063                /* olvo 980921 necessary for reverse */
3064                if (arg == res) {
3065                    IF_KEEP_WRITE_TAYLOR(arg,keep,k,p)
3066                }
3067
3068#if !defined(_NTIGHT_)
3069
3070#ifndef _ZOS_ /* BREAK_ZOS */
3071#if !defined(_INT_FOR_)
3072                T0arg   = dp_T0[arg];
3073#endif
3074#endif /* ALL_TOGETHER_AGAIN */
3075
3076                dp_T0[res] =
3077                    pow(dp_T0[arg], coval);
3078#endif /* !_NTIGHT_ */
3079
3080                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
3081#if defined(_INDO_)
3082#if defined(_INDOPRO_)
3083                copy_index_domain(res, arg, ind_dom);
3084#endif
3085#if defined(_NONLIND_)
3086                fod[opind].entry = maxopind+2;
3087                fod[opind].left = &fod[arg_index[arg]];
3088                fod[opind].right = NULL;
3089                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
3090                arg_index[res] = opind++;               
3091#endif
3092#else
3093#ifndef _ZOS_ /* BREAK_ZOS */
3094                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
3095                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
3096
3097#ifdef _INT_FOR_
3098                FOR_0_LE_l_LT_p
3099                TRES_FOINC = TARG_INC;
3100#else
3101                if (T0arg == 0.0) {
3102                    if (coval <= 0.0)
3103                        FOR_0_LE_l_LT_pk
3104                        TRES_INC = make_nan();
3105                    else {
3106                        /* coval not a whole number */
3107                        if (coval - floor(coval) != 0) {
3108                            FOR_0_LE_l_LT_p
3109                            {
3110                                i = 0;
3111                                FOR_0_LE_i_LT_k
3112                                {
3113                                    if (coval - i > 1)
3114                                    TRES_INC = 0;
3115                                    if ((coval - i < 1) && (coval - i > 0))
3116                                        TRES_INC = make_inf();
3117                                        if (coval - i < 0)
3118                                            TRES_INC = make_nan();
3119                                        }
3120                                    }
3121                                } else {
3122                        if (coval == 1) {
3123                                FOR_0_LE_l_LT_pk
3124                                TRES_INC = TARG_INC;
3125                            } else
3126                                /* coval is an int > 1 */
3127                                /* the following is not efficient but at least it works */
3128                                /* it reformulates x^n into x* ... *x n times */
3129                            {
3130                                INC_pk_1(Targ)
3131                                INC_pk_1(Tres)
3132
3133                                FOR_p_GT_l_GE_0
3134                                {
3135                                    FOR_k_GT_i_GE_0
3136                                    {
3137                                        *Tres = 0;
3138                                        DEC_TRES_FO
3139#if defined(_HIGHER_ORDER_)
3140                                        if (i == k-1) {
3141                                        zOP = dp_z+k-1;
3142                                        for(j=k-1;j>=0;j--) {
3143                                                (*zOP--) = (*Targ--);
3144                                            }
3145                                        }
3146                                        for (j=0;j<i;j++) {
3147                                        *Tres += dp_z[j] *
3148                                                     dp_z[i-j-1];
3149                                        }
3150                                        Tres--;
3151#endif /* _HIGHER_ORDER_ */
3152                                    }
3153                                }
3154                                for(ii=3;ii<=coval;ii++) {
3155                                    ASSIGN_T(Tres, TAYLOR_BUFFER[res])
3156                                    ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
3157                                    INC_pk_1(Targ)
3158                                    INC_pk_1(Tres)
3159
3160                                    FOR_p_GT_l_GE_0
3161                                    {
3162                                        FOR_k_GT_i_GE_0
3163                                        {
3164                                            *Tres = 0;
3165                                            DEC_TRES_FO
3166#if defined(_HIGHER_ORDER_)
3167                                            TresOP = Tres-i;
3168                                            for (j=0;j<i;j++)
3169                                            *Tres += TresOP[j] * dp_z[i-j-1];
3170                                            Tres--;
3171#endif /* _HIGHER_ORDER_ */
3172                                        }
3173                                    }
3174                                }
3175                        }
3176                    }
3177                }
3178            } else {
3179                r0 = 1.0 / T0arg;
3180                FOR_0_LE_l_LT_p
3181                FOR_0_LE_i_LT_k
3182                { /* olvo 980921 changed order to allow x = pow(x,n) */
3183#ifdef _HIGHER_ORDER_
3184                    zOP      = dp_z+i;
3185                    (*zOP--) = (*Targ) * r0;
3186#endif /* _HIGHER_ORDER_ */
3187
3188                    TRES_FOINC = dp_T0[res] *
3189                                 TARG_INC * coval * r0;
3190
3191#ifdef _HIGHER_ORDER_
3192                    TresOP = Tres-i;
3193
3194                    (*Tres) *= i+1;
3195                    y = coval*i -1;
3196                    for (j=0;j<i;j++) {
3197                        *Tres += (*TresOP++) * (*zOP--) * y;
3198                            y -= coval + 1;
3199                        }
3200                        *Tres++ /= (i+1);
3201#endif /* _HIGHER_ORDER_ */
3202                    }
3203                }
3204#endif
3205#endif
3206#endif /* ALL_TOGETHER_AGAIN */
3207                break;
3208
3209                /*--------------------------------------------------------------------------*/
3210            case sqrt_op:                                              /* sqrt_op */
3211                arg = get_locint_f();
3212                res = get_locint_f();
3213
3214                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3215
3216#if !defined(_NTIGHT_)
3217                dp_T0[res] = sqrt(dp_T0[arg]);
3218#endif /* !_NTIGHT_ */
3219
3220                ADOLC_OPENMP_RESTORE_THREAD_NUMBER;
3221
3222#if defined(_INDO_)
3223#if defined(_INDOPRO_)
3224                copy_index_domain(res, arg, ind_dom);
3225#endif
3226#if defined(_NONLIND_)
3227                fod[opind].entry = maxopind+2;
3228                fod[opind].left = &fod[arg_index[arg]];
3229                fod[opind].right = NULL;
3230                traverse_unary(&fod[opind], nonl_dom, &fod[opind], indcheck+1,maxopind+2);
3231                arg_index[res] = opind++;               
3232#endif
3233#else
3234#if !defined(_ZOS_) /* BREAK_ZOS */
3235                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
3236                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
3237
3238#ifdef _INT_FOR_
3239                FOR_0_LE_l_LT_p
3240                TRES_FOINC = TARG_INC;
3241#else
3242                FOR_0_LE_l_LT_p
3243                { TargOP = Targ;
3244                  if (dp_T0[arg] == 0.0)
3245                  /* Note: <=> dp_T0[res] == 0.0 */
3246              { r0 = 0.0;
3247                  FOR_0_LE_i_LT_k
3248                  { if (TARG>0.0) {
3249                        r0 = make_inf();
3250                            VEC_INC(Targ, k-i)
3251                            BREAK_FOR_I
3252                        } else
3253                            if (TARG<0.0) {
3254                            r0 = make_nan();
3255                                VEC_INC(Targ, k-i)
3256                                BREAK_FOR_I
3257                            } else
3258                                Targ++;
3259                              }
3260                          }
3261                  else {
3262                      r0 = 0.5/dp_T0[res];
3263                  }
3264                  Targ = TargOP;
3265
3266                  even = 1;
3267                  FOR_0_LE_i_LT_k
3268                  { TRES_FOINC = r0 * TARG_INC;
3269#if defined(_HIGHER_ORDER_)
3270                    TresOP  = Tres-i;
3271                    TresOP2 = Tres-1;
3272
3273                    x = 0;
3274                    for (j=1;2*j-1<i;j++)
3275                    x += (*TresOP++) * (*TresOP2--);
3276                    x *= 2;
3277                    if (!even)
3278                        x += (*TresOP) * (*TresOP2); /* !!! */
3279                        even = !even;
3280                        *Tres++ -= r0*x;
3281#endif /* _HIGHER_ORDER_ */
3282                      }
3283                    }
3284#endif
3285#endif
3286#endif /* ALL_TOGETHER_AGAIN */
3287                    break;
3288
3289            /*--------------------------------------------------------------------------*/
3290        case gen_quad:                                            /* gen_quad */
3291            arg1 = get_locint_f();
3292                arg2 = get_locint_f();
3293                res  = get_locint_f();
3294
3295#if !defined(_NTIGHT_)
3296                if (get_val_f()!=dp_T0[arg1]) {
3297                    fprintf(DIAG_OUT,
3298                            "ADOL-C Warning: forward sweep aborted; tape invalid!\n");
3299                    IF_KEEP_TAYLOR_CLOSE
3300                    end_sweep();
3301                    return -2;
3302                }
3303#endif /* !_NTIGHT_ */
3304
3305                coval = get_val_f();
3306
3307                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3308
3309#if !defined(_NTIGHT_)
3310                dp_T0[res] = coval;
3311#endif /* !_NTIGHT_ */
3312
3313#if defined(_INDO_)
3314               fprintf(DIAG_OUT,
3315                    "ADOL-C Warning: forward sweep aborted; sparse mode not available for gen_quad!\n");
3316               end_sweep();
3317               return -2;
3318#else
3319#if !defined(_ZOS_) /* BREAK_ZOS */
3320                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
3321                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
3322#ifdef _INT_FOR_
3323                FOR_0_LE_l_LT_p
3324                TRES_FOINC = TARG1_FOINC;
3325#else
3326                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
3327
3328                FOR_0_LE_l_LT_p
3329                { FOR_0_LE_i_LT_k
3330                  { /* olvo 980922 changed order to allow x = gen_quad(x) */
3331#if defined(_HIGHER_ORDER_)
3332                      zOP      = dp_z+i;
3333                      (*zOP--) = (i+1) * (*Targ1);
3334#endif /* _HIGHER_ORDER_ */
3335
3336                      TRES_FOINC = dp_T0[arg2] * TARG1_INC;
3337
3338#if defined(_HIGHER_ORDER_)
3339                      Targ2OP = Targ2;
3340
3341                      *Tres *= (i+1);
3342                      for (j=0;j<i;j++)
3343                      *Tres += (*Targ2OP++) * (*zOP--);
3344                      *Tres++ /= (i+1);
3345#endif /* _HIGHER_ORDER_ */
3346                  }
3347                  HOV_INC(Targ2, k)
3348                }
3349#endif
3350#endif
3351#endif /* ALL_TOGETHER_AGAIN */
3352                break;
3353
3354            /*--------------------------------------------------------------------------*/
3355        case min_op:                                                /* min_op */
3356            arg1  = get_locint_f();
3357                arg2  = get_locint_f();
3358                res   = get_locint_f();
3359                coval = get_val_f();
3360
3361                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3362
3363#if !defined(_NTIGHT_)
3364                /* olvo 980923 changed order to allow x = min(x,y) etc. */
3365
3366                /* olvo/mitev 980721 return value (taken from below) */
3367                if (dp_T0[arg1] > dp_T0[arg2]) {
3368                    if (coval)
3369                        MINDEC(ret_c,2);
3370                } else
3371                    if (dp_T0[arg1] < dp_T0[arg2]) {
3372                        if (!coval)
3373                            MINDEC(ret_c,2);
3374                    } else
3375                        if (arg1 != arg2)
3376                            MINDEC(ret_c,1);
3377#endif /* !_NTIGHT_ */
3378
3379#if defined (_INDO_)
3380#if defined (_INDOPRO_)
3381#if defined (_TIGHT_)
3382                    if (dp_T0[arg1] < dp_T0[arg2])
3383                        copy_index_domain(res, arg1, ind_dom);
3384                    else {
3385                        if (dp_T0[arg1] > dp_T0[arg2])
3386                            copy_index_domain(res, arg2, ind_dom);
3387                        else
3388                            combine_2_index_domains(res, arg1, arg2, ind_dom);
3389                    }
3390#else
3391                    combine_2_index_domains(res, arg1, arg2, ind_dom);
3392#endif
3393#endif
3394#if defined(_NONLIND_)
3395#ifdef _TIGHT_
3396                    if (dp_T0[arg1] < dp_T0[arg2])
3397                      {
3398                        fod[opind].entry = maxopind+2;
3399                        fod[opind].left = &fod[arg_index[arg1]];
3400                        fod[opind].right = NULL;
3401                        arg_index[res] = opind++;               
3402                      }           
3403                    else {
3404                        if (dp_T0[arg1] > dp_T0[arg2])
3405                          {
3406                            fod[opind].entry = maxopind+2;
3407                            fod[opind].left = &fod[arg_index[arg2]];
3408                            fod[opind].right = NULL;
3409                            arg_index[res] = opind++;           
3410
3411                          }               
3412                        else
3413                          {
3414                            fod[opind].entry = maxopind+2;
3415                            fod[opind].left = &fod[arg_index[arg1]];
3416                            fod[opind].right = &fod[arg_index[arg2]];
3417                            arg_index[res] = opind++;           
3418                          }
3419                    }
3420#else
3421                    fod[opind].entry = maxopind+2;
3422                    fod[opind].left = &fod[arg_index[arg1]];
3423                    fod[opind].right = &fod[arg_index[arg2]];
3424                    arg_index[res] = opind++;           
3425                    arg_index[res] = opind++;           
3426#endif
3427#endif
3428#else
3429#if !defined(_ZOS_) /* BREAK_ZOS */
3430                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
3431                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
3432                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
3433
3434#ifdef _INT_FOR_
3435#ifdef _TIGHT_
3436                Tqo = NULL;
3437                if (dp_T0[arg1] > dp_T0[arg2])
3438                    Tqo = Targ2;
3439                else
3440                    if (dp_T0[arg1] < dp_T0[arg2])
3441                        Tqo = Targ1;
3442
3443                FOR_0_LE_l_LT_p
3444                { Targ = Tqo;
3445                  if (Targ == NULL) /* e.g. T0[arg1] == T0[arg2] */
3446                    { Targ1OP = Targ1;
3447                      Targ2OP = Targ2;
3448                      if (TARG1 > TARG2)
3449                          Targ = Targ2OP;
3450                      else
3451                          if (TARG1 < TARG2)
3452                              Targ = Targ1OP;
3453                      Targ1++;
3454                      Targ2++;
3455                      if (Targ == NULL) /* e.g. both are equal */
3456                          Targ = Targ1OP;
3457                  }
3458
3459                  TRES_INC = TARG_INC;
3460
3461                  if (Tqo)
3462                    Tqo++;
3463                }
3464
3465                dp_T0[res] = MIN_ADOLC(dp_T0[arg1], dp_T0[arg2]);
3466#endif /* _TIGHT_ */
3467#ifdef _NTIGHT_
3468                TRES_INC = TARG1_INC | TARG2_INC;
3469#endif /* _NTIGHT_ */
3470#else
3471                Tqo = NULL;
3472                if (dp_T0[arg1] > dp_T0[arg2])
3473                    Tqo = Targ2;
3474                else
3475                    if (dp_T0[arg1] < dp_T0[arg2])
3476                        Tqo = Targ1;
3477
3478                FOR_0_LE_l_LT_p
3479                { Targ = Tqo;
3480                  if (Targ == NULL) /* e.g. dp_T0[arg1] ==
3481                                                                                 dp_T0[arg2] */
3482              { Targ1OP = Targ1;
3483                  Targ2OP = Targ2;
3484                  FOR_0_LE_i_LT_k
3485                  { if (TARG1 > TARG2) {
3486                        Targ = Targ2OP;
3487                        VEC_INC(Targ1, k-i)
3488                            VEC_INC(Targ2, k-i)
3489                            BREAK_FOR_I
3490                        } else
3491                            if (TARG1 < TARG2) {
3492                            Targ = Targ1OP;
3493                            VEC_INC(Targ1, k-i)
3494                                VEC_INC(Targ2, k-i)
3495                                BREAK_FOR_I
3496                            }
3497                        Targ1++;
3498                        Targ2++;
3499                      }
3500                      if (Targ == NULL) /* e.g. both are equal */
3501                          Targ = Targ1OP;
3502                  }
3503
3504                  FOR_0_LE_i_LT_k
3505                  TRES_INC = TARG_INC;
3506
3507                  if (Tqo) {
3508                  VEC_INC(Tqo, k)
3509                  }
3510            }
3511#endif
3512#endif
3513#endif /* ALL_TOGETHER_AGAIN */
3514#if !defined(_NTIGHT_)
3515                dp_T0[res] =
3516                    MIN_ADOLC( dp_T0[arg1],
3517                               dp_T0[arg2] );
3518#endif /* !_NTIGHT_ */
3519                break;
3520
3521                /*--------------------------------------------------------------------------*/
3522            case abs_val:                                              /* abs_val */
3523                arg   = get_locint_f();
3524                res   = get_locint_f();
3525                coval = get_val_f();
3526
3527                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3528
3529#if !defined(_NTIGHT_)
3530                /* olvo 980923 changed order to allow x = min(x,y) etc. */
3531
3532                /* olvo/mitev 980721 ec n3l (taken from below) */
3533                if (dp_T0[arg] < 0.0) {
3534                    if (coval)
3535                        MINDEC(ret_c,2);
3536                } else
3537                    if (dp_T0[arg] > 0.0) {
3538                        if (!coval)
3539                            MINDEC(ret_c,2);
3540                    }
3541#endif /* !_NTIGHT_ */
3542
3543#if defined(_INDO_)
3544#if defined(_INDOPRO_)
3545                copy_index_domain(res, arg, ind_dom);
3546#endif
3547#if defined(_NONLIND_)
3548                arg_index[res] = arg_index[arg];               
3549#endif
3550#else
3551#if !defined(_ZOS_) /* BREAK_ZOS */
3552                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
3553                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
3554
3555#ifdef _INT_FOR_
3556#ifdef _TIGHT_
3557                y = 0.0;
3558                if (dp_T0[arg] != 0.0) {
3559                    if (dp_T0[arg] < 0.0)
3560                        y = -1.0;
3561                    else
3562                        y = 1.0;
3563                }
3564                FOR_0_LE_l_LT_p
3565                { if ((y == 0.0) && (TARG != 0.0))
3566                  MINDEC(ret_c,1);
3567
3568                  TRES_INC = TARG_INC;
3569                }
3570
3571                dp_T0[res] = fabs(dp_T0[arg]);
3572#endif /* _TIGHT_ */
3573#ifdef _NTIGHT_
3574                FOR_0_LE_l_LT_p
3575                TRES_INC = TARG_INC;
3576#endif /* _NTIGHT_ */
3577#else
3578                y = 0.0;
3579                if (dp_T0[arg] != 0.0) {
3580                    if (dp_T0[arg] < 0.0)
3581                        y = -1.0;
3582                    else
3583                        y = 1.0;
3584                }
3585
3586                FOR_0_LE_l_LT_p
3587                { x = y;
3588                  FOR_0_LE_i_LT_k
3589                  { if ((x == 0.0) && (TARG != 0.0)) {
3590                    MINDEC(ret_c,1);
3591                        if (TARG < 0.0)
3592                            x = -1.0;
3593                        else
3594                            x = 1.0;
3595                    }
3596                    TRES_INC = x * TARG_INC;
3597              }
3598            }
3599#endif
3600#endif
3601#endif /* ALL_TOGETHER_AGAIN */
3602#if !defined(_NTIGHT_)
3603                dp_T0[res] = fabs(dp_T0[arg]);
3604#endif /* !_NTIGHT_ */
3605                break;
3606
3607                /*--------------------------------------------------------------------------*/
3608            case ceil_op:                                              /* ceil_op */
3609                arg   = get_locint_f();
3610                res   = get_locint_f();
3611                coval = get_val_f();
3612
3613                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3614
3615#if !defined(_NTIGHT_)
3616                dp_T0[res]=ceil(dp_T0[arg]);
3617                /* olvo/mitev 980721 ec n2l (taken from below) */
3618                if (coval != dp_T0[res])
3619                    MINDEC(ret_c,2);
3620#endif /* !_NTIGHT_ */
3621
3622#if defined(_INDO_)
3623#if defined(_INDOPRO_)
3624#ifdef _TIGHT_
3625                ind_dom[res][0] = 0;
3626#else
3627                copy_index_domain(res, arg, ind_dom);
3628#endif /* _TIGHT_ */
3629#endif
3630#if defined(_NONLIND_)
3631                arg_index[res] = arg_index[arg];               
3632#endif
3633#else
3634#if !defined(_ZOS_) /* BREAK_ZOS */
3635                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
3636
3637                FOR_0_LE_l_LT_pk
3638                TRES_INC = 0.0;
3639#endif
3640#endif /* ALL_TOGETHER_AGAIN */
3641                break;
3642
3643                /*--------------------------------------------------------------------------*/
3644            case floor_op:                 /* Compute ceil of adouble    floor_op */
3645                arg   = get_locint_f();
3646                res   = get_locint_f();
3647                coval = get_val_f();
3648
3649                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3650
3651#if !defined(_NTIGHT_)
3652                dp_T0[res] = floor(dp_T0[arg]);
3653                /* olvo/mitev 980721 ec n2l (taken from below) */
3654                if (coval != dp_T0[res])
3655                    MINDEC(ret_c,2);
3656#endif /* !_NTIGHT_ */
3657
3658#if defined(_INDO_)
3659#if defined(_INDOPRO_)
3660#ifdef _TIGHT_
3661                ind_dom[res][0] = 0;
3662#else
3663                copy_index_domain(res, arg, ind_dom);
3664#endif /* _TIGHT_ */
3665#endif
3666#if defined(_NONLIND_)
3667                arg_index[res] = arg_index[arg];               
3668#endif
3669#else
3670#if !defined(_ZOS_) /* BREAK_ZOS */
3671                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
3672
3673                FOR_0_LE_l_LT_pk
3674                TRES_INC = 0.0;
3675#endif
3676#endif /* ALL_TOGETHER_AGAIN */
3677                break;
3678
3679
3680                /****************************************************************************/
3681                /*                                                             CONDITIONALS */
3682
3683                /*--------------------------------------------------------------------------*/
3684            case cond_assign:                                      /* cond_assign */
3685                arg   = get_locint_f();
3686                arg1  = get_locint_f();
3687                arg2  = get_locint_f();
3688                res   = get_locint_f();
3689                coval = get_val_f();
3690
3691                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3692
3693                /* olvo 980924 changed order to allow reflexive ops */
3694#if defined (_INDO_)
3695#if defined (_INDOPRO_)
3696#if defined (_TIGHT_)
3697                if (dp_T0[arg] > 0) {
3698                    if (coval <= 0.0)
3699                        MINDEC(ret_c,2);
3700                    dp_T0[res] = dp_T0[arg1];
3701
3702                    copy_index_domain(res, arg1, ind_dom);
3703
3704                } else {
3705                    if (coval > 0.0)
3706                        MINDEC(ret_c,2);
3707                    if (dp_T0[arg] == 0)
3708                        MINDEC(ret_c,0);
3709                    dp_T0[res] = dp_T0[arg2];
3710                    copy_index_domain(res, arg2, ind_dom);
3711                }
3712#else
3713                    combine_2_index_domains(res, arg1, arg2, ind_dom);
3714#endif
3715#endif
3716#if defined (_NONLIND_)
3717#ifdef _TIGHT_
3718                if (dp_T0[arg] > 0) {
3719                    if (coval <= 0.0)
3720                        MINDEC(ret_c,2);
3721                    dp_T0[res] = dp_T0[arg1];
3722
3723                    arg_index[res] = arg_index[arg1];           
3724                } else {
3725                    if (coval > 0.0)
3726                        MINDEC(ret_c,2);
3727                    if (dp_T0[arg] == 0)
3728                        MINDEC(ret_c,0);
3729                    dp_T0[res] = dp_T0[arg2];
3730
3731                    arg_index[res] = arg_index[arg2];           
3732                }
3733
3734#else
3735               arg_index[res] = opind++;               
3736#endif
3737#endif
3738#else
3739#if !defined(_ZOS_) /* BREAK_ZOS */
3740                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
3741                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
3742                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
3743#endif /* ALL_TOGETHER_AGAIN */
3744
3745#ifdef _INT_FOR_
3746#ifdef _TIGHT_
3747                coval = get_val_f();
3748
3749                if (dp_T0[arg] > 0)
3750                    FOR_0_LE_l_LT_pk
3751                    TRES_INC = TARG1_INC;
3752                else
3753                    FOR_0_LE_l_LT_pk
3754                    TRES_INC = TARG2_INC;
3755
3756                if (dp_T0[arg] > 0) {
3757                    if (coval <= 0.0)
3758                        MINDEC(ret_c,2);
3759                    dp_T0[res] = dp_T0[arg1];
3760                } else {
3761                    if (coval > 0.0)
3762                        MINDEC(ret_c,2);
3763                    if (dp_T0[arg] == 0)
3764                        MINDEC(ret_c,0);
3765                    dp_T0[res] = dp_T0[arg2];
3766                }
3767#endif /* _TIGHT_ */
3768#ifdef _NTIGHT_
3769                FOR_0_LE_l_LT_pk
3770                TRES_INC = TARG1_INC | TARG2_INC;
3771#endif /* _NTIGHT_ */
3772#else
3773#if !defined(_ZOS_) /* BREAK_ZOS */
3774                if (dp_T0[arg] > 0)
3775                    FOR_0_LE_l_LT_pk
3776                    TRES_INC = TARG1_INC;
3777                else
3778                    FOR_0_LE_l_LT_pk
3779                    TRES_INC = TARG2_INC;
3780#endif
3781
3782                if (dp_T0[arg] > 0) {
3783                    if (coval <= 0.0)
3784                        MINDEC(ret_c,2);
3785                    dp_T0[res] = dp_T0[arg1];
3786                } else {
3787                    if (coval > 0.0)
3788                        MINDEC(ret_c,2);
3789                    if (dp_T0[arg] == 0)
3790                        MINDEC(ret_c,0);
3791                    dp_T0[res] = dp_T0[arg2];
3792                }
3793#endif
3794#endif /* ALL_TOGETHER_AGAIN */
3795                break;
3796
3797                /*--------------------------------------------------------------------------*/
3798            case cond_assign_s:                                  /* cond_assign_s */
3799                arg   = get_locint_f();
3800                arg1  = get_locint_f();
3801                res   = get_locint_f();
3802                coval = get_val_f();
3803
3804                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
3805
3806                /* olvo 980924 changed order to allow reflexive ops */
3807#if defined(_INDO_)
3808#if defined(_INDOPRO_)
3809                copy_index_domain(res, arg1, ind_dom);
3810#endif
3811#if defined(_NONLIND_)
3812                arg_index[res] = arg_index[arg1];               
3813#endif
3814#else
3815#if !defined(_ZOS_) /* BREAK_ZOS */
3816                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
3817                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
3818#endif /* ALL_TOGETHER_AGAIN */
3819
3820#ifdef _INT_FOR_
3821#ifdef _TIGHT_
3822                coval = get_val_f();
3823
3824                if (dp_T0[arg] > 0)
3825#endif /* _TIGHT_ */
3826                    FOR_0_LE_l_LT_pk
3827                    TRES_INC = TARG1_INC;
3828
3829#ifdef _TIGHT_
3830                if (dp_T0[arg] > 0) {
3831                    if (coval <= 0.0)
3832                        MINDEC(ret_c,2);
3833                    dp_T0[res] = dp_T0[arg1];
3834                } else
3835                    if (dp_T0[arg] == 0)
3836                        MINDEC(ret_c,0);
3837#endif /* _TIGHT_ */
3838#else
3839#if !defined(_ZOS_) /* BREAK_ZOS */
3840                if (dp_T0[arg] > 0)
3841                    FOR_0_LE_l_LT_pk
3842                    TRES_INC = TARG1_INC;
3843#endif
3844                if (dp_T0[arg] > 0) {
3845                    if (coval <= 0.0)
3846                        MINDEC(ret_c,2);
3847                    dp_T0[res] = dp_T0[arg1];
3848                } else
3849                    if (dp_T0[arg] == 0)
3850                        MINDEC(ret_c,0);
3851#endif
3852#endif /* ALL_TOGETHER_AGAIN */
3853                break;
3854
3855
3856                /*--------------------------------------------------------------------------*/
3857                /* NEW CONDITIONALS */
3858                /*--------------------------------------------------------------------------*/
3859            case neq_a_a:
3860            case eq_a_a:
3861            case le_a_a:
3862            case ge_a_a:
3863            case lt_a_a:
3864            case gt_a_a:
3865                coval = get_val_f();
3866                arg = get_locint_f();
3867                arg1 = get_locint_f();
3868                res = get_locint_f();
3869#if !defined(_NTIGHT_)
3870                {
3871                    revreal retval = -1;
3872                    const char* opname = "";
3873                    switch (operation) {
3874                    case neq_a_a:
3875                        retval = (revreal)(dp_T0[arg] != dp_T0[arg1]);
3876                        opname = "neq_a_a";
3877                        break;
3878                    case eq_a_a:
3879                        retval = (revreal)(dp_T0[arg] == dp_T0[arg1]);
3880                        opname = "eq_a_a";
3881                        break;
3882                    case ge_a_a:
3883                        retval = (revreal)(dp_T0[arg] >= dp_T0[arg1]);
3884                        opname = "ge_a_a";
3885                        break;
3886                    case le_a_a:
3887                        retval = (revreal)(dp_T0[arg] <= dp_T0[arg1]);
3888                        opname = "le_a_a";
3889                        break;
3890                    case gt_a_a:
3891                        retval = (revreal)(dp_T0[arg] > dp_T0[arg1]);
3892                        opname = "gt_a_a";
3893                        break;
3894                    case lt_a_a:
3895                        retval = (revreal)(dp_T0[arg] < dp_T0[arg1]);
3896                        opname = "lt_a_a";
3897                        break;
3898                    }
3899                    if (retval != coval && ADOLC_GLOBAL_TAPE_VARS.branchSwitchWarning)
3900                        fprintf(DIAG_OUT,
3901                                "ADOL-C Warning: Branch switch detected in comparison "
3902                                "(operator %s).\n"
3903                                "Results may be unpredictable! Retaping recommended!\n",opname);
3904                    dp_T0[res] = retval;
3905                }
3906#endif
3907#if defined(_INDO_)
3908#if defined(_INDOPRO_)
3909                ind_dom[res][0]=0;
3910#endif
3911#if defined(_NONLIND_)
3912                fod[opind].entry = maxopind+2;
3913                fod[opind].left = NULL;
3914                fod[opind].right = NULL;
3915                arg_index[res] = opind++;
3916#endif
3917#else
3918#if !defined(_ZOS_) /* BREAK_ZOS */
3919                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
3920
3921                FOR_0_LE_l_LT_pk
3922                TRES_INC = 0;
3923#endif
3924#endif /* ALL_TOGETHER_AGAIN */
3925
3926                break;
3927
3928                /*--------------------------------------------------------------------------*/
3929            case subscript:
3930                coval = get_val_f();
3931                arg = get_locint_f();
3932                res = get_locint_f();
3933                {
3934                    size_t cnt, idx, numvar = (size_t)trunc(fabs(coval));
3935                    locint vectorloc;
3936                    vectorloc = get_locint_f();
3937#if !defined(_NTIGHT_)
3938                    idx = (size_t)trunc(fabs(dp_T0[arg]));
3939                    if (idx >= numvar)
3940                        fprintf(DIAG_OUT, "ADOL-C warning: index out of bounds while subscripting n=%z, idx=%z\n", numvar, idx);
3941                    arg1 = vectorloc+idx;
3942                    IF_KEEP_WRITE_TAYLOR(res,keep,k,p);
3943                    dp_T0[res] = dp_T0[arg1];
3944#if defined(_INDO_)
3945#if defined(_INDOPRO_)
3946                    copy_index_domain(res, arg1, ind_dom);
3947#endif
3948#if defined(_NONLIND_)
3949                    arg_index[res] = arg_index[arg1];
3950#endif
3951#else
3952#if !defined(_ZOS_) /* BREAK_ZOS */
3953                    ASSIGN_T(Targ1,TAYLOR_BUFFER[arg1])
3954                    ASSIGN_T(Tres,TAYLOR_BUFFER[res])
3955
3956                    FOR_0_LE_l_LT_pk
3957                    TRES_INC = TARG1_INC;
3958#endif
3959#endif
3960#else
3961                    fprintf(DIAG_OUT, "ADOL-C error: active subscripting does not work in safe mode, please use tight mode\n");
3962                    exit(-2);
3963#endif /* ALL_TOGETHER_AGAIN */
3964                }
3965                break;
3966
3967            case subscript_ref:
3968                coval = get_val_f();
3969                arg = get_locint_f();
3970                res = get_locint_f();
3971                {
3972                    size_t cnt, idx, numvar = (size_t)trunc(fabs(coval));
3973                    locint vectorloc;
3974                    vectorloc = get_locint_f();
3975#if !defined(_NTIGHT_)
3976                    idx = (size_t)trunc(fabs(dp_T0[arg]));
3977                    if (idx >= numvar)
3978                        fprintf(DIAG_OUT, "ADOL-C warning: index out of bounds while subscripting (ref) n=%z, idx=%z\n", numvar, idx);
3979                    arg1 = vectorloc+idx;
3980                    IF_KEEP_WRITE_TAYLOR(res,keep,k,p);
3981                    dp_T0[res] = arg1;
3982#else
3983                    fprintf(DIAG_OUT, "ADOL-C error: active subscripting does not work in safe mode, please use tight mode\n");
3984                    exit(-2);
3985#endif
3986                }
3987                break;
3988
3989            case ref_copyout:
3990                arg = get_locint_f();
3991                res = get_locint_f();
3992#if !defined(_NTIGHT_)
3993                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
3994                IF_KEEP_WRITE_TAYLOR(res,keep,k,p);
3995                dp_T0[res] = dp_T0[arg1];
3996#if defined(_INDO_)
3997#if defined(_INDOPRO_)
3998                copy_index_domain(res, arg1, ind_dom);
3999#endif
4000#if defined(_NONLIND_)
4001                arg_index[res] = arg_index[arg1];
4002#endif
4003#else
4004#if !defined(_ZOS_) /* BREAK_ZOS */
4005                ASSIGN_T(Targ1,TAYLOR_BUFFER[arg1])
4006                ASSIGN_T(Tres,TAYLOR_BUFFER[res])
4007
4008                FOR_0_LE_l_LT_pk
4009                TRES_INC = TARG1_INC;
4010#endif
4011#endif
4012#else
4013                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4014                exit(-2);
4015#endif /* ALL_TOGETHER_AGAIN */
4016                break;
4017
4018            case ref_incr_a:
4019                arg = get_locint_f();
4020#if !defined(_NTIGHT_)
4021                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
4022                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p);
4023                dp_T0[arg1]++;
4024#else
4025                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4026                exit(-2);
4027#endif
4028                break;
4029
4030            case ref_decr_a:
4031                arg = get_locint_f();
4032#if !defined(_NTIGHT_)
4033                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
4034                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p);
4035                dp_T0[arg1]--;
4036#else
4037                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4038                exit(-2);
4039#endif
4040                break;
4041
4042            case ref_assign_d:
4043                arg = get_locint_f();
4044                coval = get_val_f();
4045               
4046#if !defined(_NTIGHT_)
4047                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
4048                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p)
4049                dp_T0[arg1] = coval;
4050#if defined(_INDO_)
4051#if defined(_INDOPRO_)
4052                ind_dom[arg1][0] = 0;
4053#endif
4054#if defined(_NONLIND_)
4055                fod[opind].entry = maxopind+2;
4056                fod[opind].left = NULL;
4057                fod[opind].right = NULL;
4058                arg_index[arg1] = opind++;
4059#endif
4060#else
4061#if !defined(_ZOS_)
4062                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
4063
4064                FOR_0_LE_l_LT_pk
4065                TARG1_INC = 0;
4066#endif
4067#endif
4068#else
4069                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4070                exit(-2);
4071#endif
4072                break;
4073
4074            case ref_assign_d_zero:
4075                arg = get_locint_f();
4076
4077#if !defined(_NTIGHT_)
4078                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
4079                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p)
4080                dp_T0[arg1] = 0.0;
4081#if defined(_INDO_)
4082#if defined(_INDOPRO_)
4083                ind_dom[arg1][0] = 0;
4084#endif
4085#if defined(_NONLIND_)
4086                fod[opind].entry = maxopind+2;
4087                fod[opind].left = NULL;
4088                fod[opind].right = NULL;
4089                arg_index[arg1] = opind++;
4090#endif
4091#else
4092#if !defined(_ZOS_)
4093                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
4094
4095                FOR_0_LE_l_LT_pk
4096                TARG1_INC = 0;
4097#endif
4098#endif
4099#else
4100                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4101                exit(-2);
4102#endif
4103                break;
4104
4105            case ref_assign_d_one:
4106                arg = get_locint_f();
4107
4108#if !defined(_NTIGHT_)
4109                arg1 = (size_t)trunc(fabs(dp_T0[arg]));
4110                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p)
4111                dp_T0[arg1] = 1.0;
4112#if defined(_INDO_)
4113#if defined(_INDOPRO_)
4114                ind_dom[arg1][0] = 0;
4115#endif
4116#if defined(_NONLIND_)
4117                fod[opind].entry = maxopind+2;
4118                fod[opind].left = NULL;
4119                fod[opind].right = NULL;
4120                arg_index[arg1] = opind++;
4121#endif
4122#else
4123#if !defined(_ZOS_)
4124                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
4125
4126                FOR_0_LE_l_LT_pk
4127                TARG1_INC = 0;
4128#endif
4129#endif
4130#else
4131                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4132                exit(-2);
4133#endif
4134                break;
4135
4136            case ref_assign_a:           /* assign an adouble variable an    assign_a */
4137                /* adouble value. (=) */
4138                arg = get_locint_f();
4139                res = get_locint_f();
4140
4141#if !defined(_NTIGHT_)
4142                arg1 = (size_t)trunc(fabs(dp_T0[res]));
4143                IF_KEEP_WRITE_TAYLOR(arg1,keep,k,p)
4144                dp_T0[arg1] = dp_T0[arg];
4145#if defined(_INDO_)
4146#if defined(_INDOPRO_)
4147                copy_index_domain(arg1, arg, ind_dom);
4148#endif
4149#if defined(_NONLIND_)
4150                arg_index[arg1] = arg_index[arg];
4151#endif
4152#else
4153#if !defined(_ZOS_) /* BREAK_ZOS */
4154                ASSIGN_T(Targ,TAYLOR_BUFFER[arg])
4155                ASSIGN_T(Targ1,TAYLOR_BUFFER[arg1])
4156
4157                FOR_0_LE_l_LT_pk
4158                TARG1_INC = TARG_INC;
4159#endif
4160#endif
4161#else
4162                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4163                exit(-2);
4164#endif /* ALL_TOGETHER_AGAIN */
4165                break;
4166
4167            case ref_assign_ind:       /* assign an adouble variable an    assign_ind */
4168                /* independent double value (<<=) */
4169                arg = get_locint_f();
4170
4171
4172#if !defined(_NTIGHT_)
4173                res = (size_t)trunc(fabs(dp_T0[arg]));
4174                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4175                dp_T0[res] = basepoint[indexi];
4176#if defined(_INDO_)
4177#if defined(_INDOPRO_)
4178                ind_dom[res][0] = 1;
4179                ind_dom[res][2] = indexi;
4180#endif
4181#if defined(_NONLIND_)
4182                fod[opind].entry = indexi;
4183                fod[opind].left = NULL;
4184                fod[opind].right = NULL;
4185                arg_index[res] = opind++;
4186#endif
4187#else
4188#if !defined(_ZOS_) /* BREAK_ZOS */
4189                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
4190
4191#ifdef _INT_FOR_
4192                FOR_0_LE_l_LT_p
4193                TRES_INC = ARGUMENT(indexi,l,i);
4194#else
4195                FOR_0_LE_l_LT_p
4196                FOR_0_LE_i_LT_k
4197                TRES_INC = ARGUMENT(indexi,l,i);
4198#endif
4199#endif
4200#endif
4201#else
4202                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4203                exit(-2);
4204#endif /* ALL_TOGETHER_AGAIN */
4205                ++indexi;
4206                break;
4207
4208            case ref_eq_plus_d:            /* Add a floating point to an    eq_plus_d */
4209                /* adouble. (+=) */
4210                arg  = get_locint_f();
4211                coval = get_val_f();
4212
4213
4214#if !defined(_NTIGHT_)
4215                res = (size_t)trunc(fabs(dp_T0[arg]));
4216                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4217                dp_T0[res] += coval;
4218#else
4219                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4220                exit(-2);
4221#endif /* !_NTIGHT_ */
4222                break;
4223
4224                /*--------------------------------------------------------------------------*/
4225            case ref_eq_plus_a:             /* Add an adouble to another    eq_plus_a */
4226                /* adouble. (+=) */
4227                arg = get_locint_f();
4228                arg1 = get_locint_f();
4229
4230#if !defined(_NTIGHT_)
4231                res = (size_t)trunc(fabs(dp_T0[arg1]));
4232                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4233                dp_T0[res] += dp_T0[arg];
4234#if defined(_INDO_)
4235#if defined(_INDOPRO_)
4236                merge_2_index_domains(res, arg, ind_dom);
4237#endif
4238#if defined(_NONLIND_)
4239                fod[opind].entry = maxopind+2;
4240                fod[opind].left = &fod[arg_index[res]];
4241                fod[opind].right = &fod[arg_index[arg]];
4242                arg_index[res] = opind++;
4243#endif
4244#else
4245#if !defined(_ZOS_) /* BREAK_ZOS */
4246                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
4247                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
4248
4249#ifdef _INT_FOR_
4250                FOR_0_LE_l_LT_pk
4251                TRES_INC |= TARG_INC;
4252#else
4253                FOR_0_LE_l_LT_pk
4254                TRES_INC += TARG_INC;
4255#endif
4256#endif
4257#endif
4258#else
4259                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4260                exit(-2);
4261#endif /* ALL_TOGETHER_AGAIN */
4262                break;
4263
4264            case ref_eq_min_d:       /* Subtract a floating point from an    eq_min_d */
4265                /* adouble. (-=) */
4266                arg = get_locint_f();
4267                coval = get_val_f();
4268
4269#if !defined(_NTIGHT_)
4270                res = (size_t)trunc(fabs(dp_T0[arg]));
4271                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4272                dp_T0[res] -= coval;
4273#else
4274                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4275                exit(-2);
4276#endif /* !_NTIGHT_ */
4277                break;
4278
4279                /*--------------------------------------------------------------------------*/
4280            case ref_eq_min_a:        /* Subtract an adouble from another    eq_min_a */
4281                /* adouble. (-=) */
4282                arg = get_locint_f();
4283                arg1 = get_locint_f();
4284
4285#if !defined(_NTIGHT_)
4286                res = (size_t)trunc(fabs(dp_T0[arg1]));
4287                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4288                dp_T0[res] -= dp_T0[arg];
4289#if defined(_INDO_)
4290#if defined(_INDOPRO_)
4291                merge_2_index_domains(res, arg, ind_dom);
4292#endif
4293#if defined(_NONLIND_)
4294                fod[opind].entry = maxopind+2;
4295                fod[opind].left = &fod[arg_index[res]];
4296                fod[opind].right = &fod[arg_index[arg]];
4297                arg_index[res] = opind++;
4298#endif
4299#else
4300#if !defined(_ZOS_) /* BREAK_ZOS */
4301                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
4302                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
4303
4304#ifdef _INT_FOR_
4305                FOR_0_LE_l_LT_pk
4306                TRES_INC |= TARG_INC;
4307#else
4308                FOR_0_LE_l_LT_pk
4309                TRES_INC -= TARG_INC;
4310#endif
4311#endif
4312#endif
4313#else
4314                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4315                exit(-2);
4316#endif /* ALL_TOGETHER_AGAIN */
4317                break;
4318
4319            case ref_eq_mult_d:              /* Multiply an adouble by a    eq_mult_d */
4320                /* flaoting point. (*=) */
4321                arg = get_locint_f();
4322                coval = get_val_f();
4323
4324#if !defined(_NTIGHT_)
4325                res = (size_t)trunc(fabs(dp_T0[arg]));
4326                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4327                dp_T0[res] *= coval;
4328#if !defined(_INDO_)
4329#if !defined(_ZOS_) /* BREAK_ZOS */
4330#if !defined( _INT_FOR_)
4331
4332                FOR_0_LE_l_LT_pk
4333                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
4334
4335                FOR_0_LE_l_LT_pk
4336                TRES_INC *= coval;
4337#endif
4338#endif
4339#endif
4340#else
4341                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4342                exit(-2);
4343#endif /* ALL_TOGETHER_AGAIN */
4344                break;
4345
4346            case ref_eq_mult_a:       /* Multiply one adouble by another    eq_mult_a */
4347                /* (*=) */
4348                arg = get_locint_f();
4349                arg1 = get_locint_f();
4350
4351#if !defined(_NTIGHT_)
4352                res = (size_t)trunc(fabs(dp_T0[arg1]));
4353                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4354#if defined(_INDO_)
4355#if defined(_INDOPRO_)
4356                merge_2_index_domains(res, arg, ind_dom);
4357#endif
4358#if defined(_NONLIND_)
4359                fod[opind].entry = maxopind+2;
4360                fod[opind].left = &fod[arg_index[res]];
4361                fod[opind].right = &fod[arg_index[arg]];
4362                traverse_unary(&fod[arg_index[res]], nonl_dom, &fod[arg_index[arg]], indcheck+1,maxopind+2);
4363                traverse_unary(&fod[arg_index[arg]], nonl_dom, &fod[arg_index[res]], indcheck+1,maxopind+2);
4364                arg_index[res] = opind++;
4365#endif
4366#else
4367#if !defined(_ZOS_) /* BREAK_ZOS */
4368                ASSIGN_T(Tres, TAYLOR_BUFFER[res])
4369                ASSIGN_T(Targ, TAYLOR_BUFFER[arg])
4370
4371                INC_pk_1(Tres)
4372                INC_pk_1(Targ)
4373
4374#ifdef _INT_FOR_
4375                FOR_p_GT_l_GE_0
4376                TRES_FODEC |= TARG_DEC;
4377#else
4378                FOR_p_GT_l_GE_0
4379                FOR_k_GT_i_GE_0
4380                { TRES_FODEC = dp_T0[res]*TARG_DEC +
4381                               TRES*dp_T0[arg];
4382                  DEC_TRES_FO
4383#ifdef _HIGHER_ORDER_
4384                  TresOP = Tres-i;
4385                  TargOP = Targ;
4386
4387                  for (j=0;j<i;j++)
4388                  *Tres += (*TresOP++) * (*TargOP--);
4389                  Tres--;
4390#endif /* _HIGHER_ORDER_ */
4391                }
4392#endif
4393#endif
4394#endif
4395                dp_T0[res] *= dp_T0[arg];
4396#else
4397                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4398                exit(-2);
4399#endif /* !_NTIGHT_ */
4400                break;
4401
4402            case ref_cond_assign:                                      /* cond_assign */
4403                arg   = get_locint_f();
4404                arg1  = get_locint_f();
4405                arg2  = get_locint_f();
4406                { 
4407                    locint ref = get_locint_f();
4408                    coval = get_val_f();
4409#if !defined(_NTIGHT_)
4410                    res   = (size_t)trunc(fabs(dp_T0[ref]));
4411
4412                    IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4413
4414                /* olvo 980924 changed order to allow reflexive ops */
4415#if defined(_INDO_)
4416                if (dp_T0[arg] > 0) {
4417                    if (coval <= 0.0)
4418                        MINDEC(ret_c,2);
4419                    dp_T0[res] = dp_T0[arg1];
4420
4421#if defined(_INDOPRO_)
4422                    copy_index_domain(res, arg1, ind_dom);
4423#endif
4424#if defined(_NONLIND_)
4425                    arg_index[res] = arg_index[arg1];
4426#endif
4427                } else {
4428                    if (coval > 0.0)
4429                        MINDEC(ret_c,2);
4430                    if (dp_T0[arg] == 0)
4431                        MINDEC(ret_c,0);
4432                    dp_T0[res] = dp_T0[arg2];
4433
4434#if defined(_INDOPRO_)
4435                    copy_index_domain(res, arg2, ind_dom);
4436#endif
4437#if defined(_NONLIND_)
4438                    arg_index[res] = arg_index[arg2];
4439#endif
4440                }
4441#else
4442#if !defined(_ZOS_) /* BREAK_ZOS */
4443                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
4444                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
4445                ASSIGN_T(Targ2, TAYLOR_BUFFER[arg2])
4446#endif /* ALL_TOGETHER_AGAIN */
4447
4448#ifdef _INT_FOR_
4449                coval = get_val_f();
4450
4451                if (dp_T0[arg] > 0)
4452                    FOR_0_LE_l_LT_pk
4453                    TRES_INC = TARG1_INC;
4454                else
4455                    FOR_0_LE_l_LT_pk
4456                    TRES_INC = TARG2_INC;
4457
4458                if (dp_T0[arg] > 0) {
4459                    if (coval <= 0.0)
4460                        MINDEC(ret_c,2);
4461                    dp_T0[res] = dp_T0[arg1];
4462                } else {
4463                    if (coval > 0.0)
4464                        MINDEC(ret_c,2);
4465                    if (dp_T0[arg] == 0)
4466                        MINDEC(ret_c,0);
4467                    dp_T0[res] = dp_T0[arg2];
4468                }
4469                FOR_0_LE_l_LT_pk
4470                TRES_INC = TARG1_INC | TARG2_INC;
4471#else
4472#if !defined(_ZOS_) /* BREAK_ZOS */
4473                if (dp_T0[arg] > 0)
4474                    FOR_0_LE_l_LT_pk
4475                    TRES_INC = TARG1_INC;
4476                else
4477                    FOR_0_LE_l_LT_pk
4478                    TRES_INC = TARG2_INC;
4479#endif
4480
4481                if (dp_T0[arg] > 0) {
4482                    if (coval <= 0.0)
4483                        MINDEC(ret_c,2);
4484                    dp_T0[res] = dp_T0[arg1];
4485                } else {
4486                    if (coval > 0.0)
4487                        MINDEC(ret_c,2);
4488                    if (dp_T0[arg] == 0)
4489                        MINDEC(ret_c,0);
4490                    dp_T0[res] = dp_T0[arg2];
4491                }
4492#endif
4493#endif
4494#else
4495                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4496                exit(-2);
4497#endif /* ALL_TOGETHER_AGAIN */
4498                }
4499                break;
4500
4501            case ref_cond_assign_s:                                  /* cond_assign_s */
4502                arg   = get_locint_f();
4503                arg1  = get_locint_f();
4504                arg2   = get_locint_f();
4505                coval = get_val_f();
4506
4507#if !defined(_NTIGHT_)
4508                res = (size_t)trunc(fabs(dp_T0[arg2]));
4509                IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
4510
4511                /* olvo 980924 changed order to allow reflexive ops */
4512#if defined(_INDO_)
4513                if (dp_T0[arg] > 0) {
4514#if defined(_INDOPRO_)
4515                    copy_index_domain(res, arg1, ind_dom);
4516#endif
4517#if defined(_NONLIND_)
4518                    arg_index[res] = arg_index[arg1];
4519#endif
4520                }
4521#else
4522#if !defined(_ZOS_) /* BREAK_ZOS */
4523                ASSIGN_T(Tres,  TAYLOR_BUFFER[res])
4524                ASSIGN_T(Targ1, TAYLOR_BUFFER[arg1])
4525#endif /* ALL_TOGETHER_AGAIN */
4526
4527#ifdef _INT_FOR_
4528                coval = get_val_f();
4529
4530                if (dp_T0[arg] > 0)
4531                    FOR_0_LE_l_LT_pk
4532                    TRES_INC = TARG1_INC;
4533
4534                if (dp_T0[arg] > 0) {
4535                    if (coval <= 0.0)
4536                        MINDEC(ret_c,2);
4537                    dp_T0[res] = dp_T0[arg1];
4538                } else
4539                    if (dp_T0[arg] == 0)
4540                        MINDEC(ret_c,0);
4541#else
4542#if !defined(_ZOS_) /* BREAK_ZOS */
4543                if (dp_T0[arg] > 0)
4544                    FOR_0_LE_l_LT_pk
4545                    TRES_INC = TARG1_INC;
4546#endif
4547                if (dp_T0[arg] > 0) {
4548                    if (coval <= 0.0)
4549                        MINDEC(ret_c,2);
4550                    dp_T0[res] = dp_T0[arg1];
4551                } else
4552                    if (dp_T0[arg] == 0)
4553                        MINDEC(ret_c,0);
4554#endif
4555#endif
4556#else
4557                fprintf(DIAG_OUT, "ADOL-C error: active vector element referencing does not work in safe mode, please use tight mode\n");
4558                exit(-2);
4559#endif /* ALL_TOGETHER_AGAIN */
4560                break;
4561
4562                /****************************************************************************/
4563                /*                                                          REMAINING STUFF */
4564
4565                /*--------------------------------------------------------------------------*/
4566            case take_stock_op:                                  /* take_stock_op */
4567                size = get_locint_f();
4568                res  = get_locint_f();
4569                d    = get_val_v_f(size);
4570
4571                for (ls=0;ls<size;ls++) {
4572#if !defined(_NTIGHT_)
4573                    dp_T0[res]=*d;
4574#endif /* !_NTIGHT_ */
4575#if !defined(_INDO_)
4576#if !defined(_ZOS_) /* BREAK_ZOS */
4577                    ASSIGN_T(Tres,TAYLOR_BUFFER[res])
4578
4579                    FOR_0_LE_l_LT_pk
4580                    TRES_INC = 0;
4581
4582#endif /* ALL_TOGETHER_AGAIN */
4583                    res++;
4584#if !defined(_NTIGHT_)
4585                    d++;
4586#endif /* !_NTIGHT_ */
4587#endif
4588                }
4589                break;
4590
4591                /*--------------------------------------------------------------------------*/
4592            case death_not:                                          /* death_not */
4593                arg1=get_locint_f();
4594                arg2=get_locint_f();
4595
4596#ifdef _KEEP_
4597                if (keep) {
4598                    do {
4599                        IF_KEEP_WRITE_TAYLOR(arg2,keep,k,p)
4600                    } while(arg1 < arg2-- );
4601                }
4602#endif
4603                break;
4604
4605                /*--------------------------------------------------------------------------*/
4606#if defined(_EXTERN_) /* ZOS,  FOS, FOV up to now */
4607            case ext_diff:                       /* extern differntiated function */
4608                ADOLC_CURRENT_TAPE_INFOS.ext_diff_fct_index=get_locint_f();
4609                n=get_locint_f();
4610                m=get_locint_f();
4611                ADOLC_CURRENT_TAPE_INFOS.lowestXLoc_for = get_locint_f();
4612                ADOLC_CURRENT_TAPE_INFOS.lowestYLoc_for = get_locint_f();
4613                ADOLC_CURRENT_TAPE_INFOS.cpIndex = get_locint_f();
4614                edfct=get_ext_diff_fct(ADOLC_CURRENT_TAPE_INFOS.ext_diff_fct_index);
4615
4616                if (edfct->ADOLC_EXT_FCT_POINTER==NULL)
4617                    fail(ADOLC_EXT_DIFF_NULLPOINTER_DIFFFUNC);
4618                if (n>0) {
4619                    if (edfct->dp_x==NULL) fail(ADOLC_EXT_DIFF_NULLPOINTER_ARGUMENT);
4620#if !defined(_ZOS_)
4621                    if (ADOLC_EXT_POINTER_X==NULL) fail(ADOLC_EXT_DIFF_NULLPOINTER_ARGUMENT);
4622#endif
4623                }
4624                if (m>0) {
4625                    if (edfct->dp_y==NULL) fail(ADOLC_EXT_DIFF_NULLPOINTER_ARGUMENT);
4626#if !defined(_ZOS_)
4627                    if (ADOLC_EXT_POINTER_Y==NULL) fail(ADOLC_EXT_DIFF_NULLPOINTER_ARGUMENT);
4628#endif
4629                }
4630
4631                arg = ADOLC_CURRENT_TAPE_INFOS.lowestXLoc_for;
4632                for (loop=0; loop<n; ++loop) {
4633                    if (edfct->dp_x_changes) {
4634                      IF_KEEP_WRITE_TAYLOR(arg, keep, k, p);
4635                    }
4636                    edfct->dp_x[loop]=dp_T0[arg];
4637#if !defined(_ZOS_)
4638                    ADOLC_EXT_LOOP
4639                        ADOLC_EXT_POINTER_X[loop]ADOLC_EXT_SUBSCRIPT =
4640                        TAYLOR_BUFFER[arg]ADOLC_EXT_SUBSCRIPT;
4641#endif
4642                    ++arg;
4643                }
4644                arg = ADOLC_CURRENT_TAPE_INFOS.lowestYLoc_for;
4645                for (loop=0; loop<m; ++loop) {
4646                    if (edfct->dp_y_priorRequired) {
4647                      IF_KEEP_WRITE_TAYLOR(arg, keep, k, p);
4648                    }
4649                    edfct->dp_y[loop]=dp_T0[arg];
4650#if !defined(_ZOS_)
4651                    ADOLC_EXT_LOOP
4652                        ADOLC_EXT_POINTER_Y[loop]ADOLC_EXT_SUBSCRIPT =
4653                        TAYLOR_BUFFER[arg]ADOLC_EXT_SUBSCRIPT;
4654#endif
4655                    ++arg;
4656                }
4657
4658                ext_retc = edfct->ADOLC_EXT_FCT_COMPLETE;
4659                MINDEC(ret_c, ext_retc);
4660
4661                res = ADOLC_CURRENT_TAPE_INFOS.lowestXLoc_for;
4662                for (loop=0; loop<n; ++loop) {
4663                    dp_T0[res]=edfct->dp_x[loop];
4664#if !defined(_ZOS_)
4665                    ADOLC_EXT_LOOP
4666                        TAYLOR_BUFFER[res]ADOLC_EXT_SUBSCRIPT =
4667                        ADOLC_EXT_POINTER_X[loop]ADOLC_EXT_SUBSCRIPT;
4668#endif
4669                    ++res;
4670                }
4671                res = ADOLC_CURRENT_TAPE_INFOS.lowestYLoc_for;
4672                for (loop=0; loop<m; ++loop) {
4673                    dp_T0[res]=edfct->dp_y[loop];
4674#if !defined(_ZOS_)
4675                    ADOLC_EXT_LOOP
4676                        TAYLOR_BUFFER[res]ADOLC_EXT_SUBSCRIPT =
4677                        ADOLC_EXT_POINTER_Y[loop]ADOLC_EXT_SUBSCRIPT;
4678#endif
4679                    ++res;
4680                }
4681
4682                break;
4683#endif
4684
4685                /*--------------------------------------------------------------------------*/
4686
4687            default:                                                   /* default */
4688                /* Die here, we screwed up */
4689
4690                fprintf(DIAG_OUT,"ADOL-C fatal error in " GENERATED_FILENAME " ("
4691                        __FILE__
4692                        ") : no such operation %d\n", operation);
4693                exit(-1);
4694                break;
4695
4696        } /* endswitch */
4697
4698        /* Read the next operation */
4699        operation=get_op_f();
4700#if defined(ADOLC_DEBUG)
4701        ++countPerOperation[operation];
4702#endif /* ADOLC_DEBUG */
4703    }  /* endwhile */
4704
4705
4706#if defined(ADOLC_DEBUG)
4707    printf("\nTape contains:\n");
4708    for (v = 0; v < 256; ++v)
4709        if (countPerOperation[v] > 0)
4710            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]);
4711    printf("\n");
4712#endif /* ADOLC_DEBUG */
4713
4714#if defined(_KEEP_)
4715    if (keep) taylor_close(taylbuf);
4716#endif
4717
4718    /* clean up */
4719#if !defined (_NTIGHT_)
4720    free(dp_T0);
4721#endif /* !_NTIGHT_ */
4722#if !defined(_INDO_)
4723#if !defined(_ZOS_)
4724#   if defined(_FOS_)
4725    free(dp_T);
4726#   else
4727#if !defined (_INT_FOR_)
4728    myfree2(dpp_T);
4729    free(dp_Ttemp);
4730#endif /* !_NTIGHT_ */
4731#endif
4732#endif
4733#endif
4734#if defined(_HIGHER_ORDER_)
4735    free(dp_z);
4736#endif
4737
4738    end_sweep();
4739
4740
4741#if defined(_INDO_)
4742#if defined(_INDOPRO_)
4743    for(i=0;i<max_ind_dom;i++)
4744      {
4745        free(ind_dom[i]);
4746      }
4747    free(ind_dom);
4748#endif
4749#if defined(_NONLIND_)
4750    for( i=0; i < indcheck; i++) {
4751      traverse_crs(&nonl_dom[i],&sod[i],indcheck+1);
4752      free_tree(&nonl_dom[i],indcheck+1);
4753      crs[i] = (unsigned int*) malloc(sizeof(unsigned int) * (sod[i].entry+1));
4754      crs[i][0] = sod[i].entry;
4755      temp = sod[i].left;
4756      for( ii=1; ii <=sod[i].entry; ii++)
4757        {
4758          crs[i][ii] = temp->entry;
4759          temp1 = temp->left;
4760          free(temp);
4761          temp = temp1;
4762        }
4763    }
4764
4765    free(sod);
4766    free(nonl_dom);
4767    free(fod);
4768    free(arg_index);
4769
4770#endif
4771#endif
4772    return ret_c;
4773}
4774
4775/****************************************************************************/
4776
4777#if defined(_INDOPRO_)
4778
4779/****************************************************************************/
4780/* set operations for propagation of index domains                          */
4781
4782/*--------------------------------------------------------------------------*/
4783/* operations on index domains                                              */
4784
4785#if defined(_TIGHT_)
4786void copy_index_domain(int res, int arg, locint **ind_dom) {
4787
4788   int i;
4789
4790   if (ind_dom[arg][0] > ind_dom[res][1])
4791     {
4792       free(ind_dom[res]);
4793       ind_dom[res] = (locint *)  malloc(sizeof(locint) * 2*(ind_dom[arg][0]+1));
4794       ind_dom[res][1] = 2*ind_dom[arg][0];
4795     }
4796
4797
4798    for(i=2;i<ind_dom[arg][0]+2;i++)
4799       ind_dom[res][i] = ind_dom[arg][i];
4800    ind_dom[res][0] = ind_dom[arg][0];
4801}
4802
4803
4804void merge_2_index_domains(int res, int arg, locint **ind_dom)
4805{
4806
4807  int num,num1,num2, i,j,k,l;
4808  locint *temp_array, *arg_ind_dom, *res_ind_dom;
4809
4810  if (ind_dom[res][0] == 0)
4811    copy_index_domain(res,arg,ind_dom);
4812  else
4813    {
4814      if (res != arg)
4815     {
4816       arg_ind_dom = ind_dom[arg];
4817       res_ind_dom = ind_dom[res];
4818
4819       num  = ind_dom[res][0];
4820       num1 = arg_ind_dom[0];
4821       num2 = ind_dom[res][1];
4822
4823       if (num2 < num1+num)
4824         num2 = num1+num;
4825
4826       temp_array = (locint *)  malloc(sizeof(locint)* (num2+2));
4827       temp_array[1] = num2;
4828
4829       i = 2;
4830       j = 2;
4831       k = 2;
4832       num += 2;
4833       num1 += 2;
4834       while ((i< num) && (j < num1))
4835         {
4836           if (res_ind_dom[i] < arg_ind_dom[j])
4837          {
4838            temp_array[k] = res_ind_dom[i];
4839            i++; k++;
4840          }
4841           else
4842          {
4843            if (res_ind_dom[i] == arg_ind_dom[j])
4844              {
4845                temp_array[k] = arg_ind_dom[j];
4846                i++;j++;k++;
4847              }
4848            else
4849              {
4850                temp_array[k] = arg_ind_dom[j];
4851                j++;k++;
4852              }
4853          }
4854         }
4855       for(l = i;l<num;l++)
4856         {
4857           temp_array[k] = res_ind_dom[l];
4858           k++;
4859         }
4860       for(l = j;l<num1;l++)
4861         {
4862           temp_array[k] = arg_ind_dom[l];
4863           k++;
4864         }
4865       temp_array[0] = k-2;
4866       free(ind_dom[res]);
4867       ind_dom[res]=temp_array;
4868     }
4869    }
4870
4871
4872}
4873
4874void combine_2_index_domains(int res, int arg1, int arg2, locint **ind_dom) {
4875
4876    if (res != arg1)
4877       copy_index_domain(res, arg1, ind_dom);
4878
4879    merge_2_index_domains(res, arg2, ind_dom);
4880}
4881
4882void merge_3_index_domains(int res, int arg1, int arg2, locint **ind_dom) {
4883    merge_2_index_domains(res, arg1, ind_dom);
4884    merge_2_index_domains(res, arg2, ind_dom);
4885}
4886
4887
4888
4889#endif
4890#endif
4891
4892
4893#if defined(_NONLIND_)
4894#if defined(_TIGHT_)
4895
4896void free_tree(IndexElement* tree, int num)
4897{
4898
4899  if (tree->left != NULL)
4900    {
4901      free_tree(tree->left,num);
4902    }
4903  if (tree->right != NULL)
4904    {
4905      free_tree(tree->right,num);
4906     }
4907    {
4908      if (tree->entry == num)
4909        free(tree);
4910
4911    }
4912 
4913}
4914void traverse_crs(IndexElement* tree,  IndexElement_sod* sod, int num)
4915{
4916
4917  IndexElement_sod *temp, *temp1;
4918  int ii;
4919
4920  if (tree->left != NULL)
4921    {
4922      traverse_crs(tree->left, sod, num);
4923    }
4924  if (tree->right != NULL)
4925    {
4926      traverse_crs(tree->right, sod, num);
4927    }
4928  if (tree->entry < num)
4929    {
4930      temp = sod->left;
4931      if (temp == NULL)
4932        {
4933          temp = (struct IndexElement_sod*) malloc(sizeof(struct IndexElement_sod));
4934          temp->left = NULL;
4935          temp->entry = tree->entry;
4936          sod->entry++;
4937          sod->left=temp;
4938        }
4939      else
4940        {
4941          while ((temp->entry < tree->entry) && (temp->left != NULL))
4942            {
4943              temp1 = temp;
4944              temp = temp->left;
4945            }
4946          if (temp->left == NULL)
4947            {
4948              if(temp->entry < tree->entry)
4949                {
4950                  temp->left = (struct IndexElement_sod*) malloc(sizeof(struct IndexElement_sod));
4951                  temp = temp->left;
4952                  temp->left = NULL;
4953                  temp->entry = tree->entry;
4954                  sod->entry++;
4955                }
4956              if(temp->entry > tree->entry)
4957                {
4958                  temp->left = (struct IndexElement_sod*) malloc(sizeof(struct IndexElement_sod));
4959                  temp->left->entry = temp->entry;
4960                  temp->left->left = NULL;
4961                  temp->entry = tree->entry;
4962                  sod->entry++;
4963                }
4964            }
4965          else
4966            {
4967              if (temp->entry > tree->entry)
4968                {
4969                  temp1 = (struct IndexElement_sod*) malloc(sizeof(struct IndexElement_sod));
4970                  temp1->left = temp->left;
4971                  temp1->entry = temp->entry;
4972                  temp->entry = tree->entry;
4973                  temp->left=temp1;
4974                  sod->entry++;
4975                }
4976             
4977            }
4978        }
4979    }
4980}
4981
4982void traverse_unary(IndexElement* tree,  IndexElement* nonl_dom,  IndexElement* fodi, int num, int maxopind)
4983{
4984  IndexElement *temp;
4985
4986  if (tree->left != NULL)
4987    {
4988      traverse_unary(tree->left, nonl_dom, fodi, num, maxopind);
4989      if (tree->right != NULL)
4990        {
4991          traverse_unary(tree->right, nonl_dom, fodi, num, maxopind);
4992        }
4993     }
4994  else
4995    {
4996      if(tree->entry<maxopind)
4997        {
4998          temp = (struct IndexElement*) malloc(sizeof(struct IndexElement));
4999          temp->right = fodi;
5000          temp->left = nonl_dom[tree->entry].left;
5001          temp->entry= num;
5002          nonl_dom[tree->entry].left = temp;
5003        }
5004    }
5005}
5006
5007#endif
5008#endif
5009END_C_DECLS
Note: See TracBrowser for help on using the repository browser.