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

Last change on this file since 340 was 340, checked in by utke, 8 years ago

apparently is defined also for FOV

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