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

Last change on this file since 56 was 56, checked in by awalther, 10 years ago

const declaration for preparation of ipopt interface

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