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

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

fix memory deallocation in uni5_for.c

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