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

Last change on this file since 191 was 191, checked in by awalther, 8 years ago

corrected zos_forward(1,...) to zos_forward(tag,...)

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