source: stable/2.2/ADOL-C/src/uni5_for.c @ 230

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

Merge stable version with revision 229

  • Property svn:keywords set to Author Date Id Revision
File size: 124.7 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 230 2011-06-08 12:12:19Z kulshres $
8
9 Contents: Contains the routines :
10           zos_forward (zero-order-scalar forward mode):      define _ZOS_
11           fos_forward (first-order-scalar forward mode):     define _FOS_
12           hos_forward (higher-order-scalar forward mode):    define _HOS_
13           fov_forward (first-order-vector forward mode):     define _FOV_
14           hov_forward (higher-order-vector forward mode):    define _HOV_
15           hov_wk_forward (higher-order-vector forward mode): define _HOV_WK_
16           int_forward_safe:                                  define _INT_FOR_ and _NTIGHT__
17
18           Uses the preprocessor to compile the 7 different object files
19           with/without "keep" parameter:                     define _KEEP_
20
21 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz,
22               Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel,
23               Benjamin Letschert
24
25 This file is part of ADOL-C. This software is provided as open source.
26 Any use, reproduction, or distribution of the software constitutes
27 recipient's acceptance of the terms of the accompanying license file.
28
29----------------------------------------------------------------------------*/
30
31#include <adolc/interfaces.h>
32#include <adolc/adalloc.h>
33#include <adolc/taping.h>
34#include "taping_p.h"
35#include <adolc/oplate.h>
36#include <adolc/externfcts.h>
37#include "externfcts_p.h"
38
39#include <math.h>
40
41#if defined(ADOLC_DEBUG)
42#include <string.h>
43#endif /* ADOLC_DEBUG */
44
45
46/****************************************************************************/
47/*                                                                   MACROS */
48#undef _ADOLC_VECTOR_
49#undef _HIGHER_ORDER_
50
51
52
53/*--------------------------------------------------------------------------*/
54#if defined(_ZOS_)
55#  define GENERATED_FILENAME "zos_forward"
56
57/*--------------------------------------------------------------------------*/
58#else
59#if defined(_FOS_)
60#define GENERATED_FILENAME "fos_forward"
61
62#define ARGUMENT(indexi,l,i) argument[indexi]
63#define TAYLORS(indexd,l,i)   taylors[indexd]
64
65/*--------------------------------------------------------------------------*/
66#else
67#if defined(_FOV_)
68#define GENERATED_FILENAME "fov_forward"
69
70#define _ADOLC_VECTOR_
71
72#if defined(_CHUNKED_)
73#define ARGUMENT(indexi,l,i) argument[indexi][l+offset]
74#define TAYLORS(indexd,l,i)   taylors[indexd][l+offset]
75#else
76#define ARGUMENT(indexi,l,i) argument[indexi][l]
77#define TAYLORS(indexd,l,i)   taylors[indexd][l]
78#endif
79
80/*--------------------------------------------------------------------------*/
81#else
82#if defined(_HOS_)
83#define GENERATED_FILENAME "hos_forward"
84
85#define _HIGHER_ORDER_
86
87#define ARGUMENT(indexi,l,i) argument[indexi][i]
88#define TAYLORS(indexd,l,i)   taylors[indexd][i]
89
90/*--------------------------------------------------------------------------*/
91#else
92#if defined(_HOV_)
93#define GENERATED_FILENAME "hov_forward"
94
95#define _ADOLC_VECTOR_
96#define _HIGHER_ORDER_
97
98#define ARGUMENT(indexi,l,i) argument[indexi][l][i]
99#define TAYLORS(indexd,l,i)   taylors[indexd][l][i]
100
101/*--------------------------------------------------------------------------*/
102#else
103#if defined(_HOV_WK_)
104#define GENERATED_FILENAME "hov_wk_forward"
105
106#define _ADOLC_VECTOR_
107#define _HIGHER_ORDER_
108
109#define ARGUMENT(indexi,l,i) argument[indexi][l][i]
110#define TAYLORS(indexd,l,i)   taylors[indexd][l][i]
111
112/*--------------------------------------------------------------------------*/
113#else
114#if defined(_INT_FOR_)
115#if defined(_TIGHT_)
116#define GENERATED_FILENAME "int_forward_t"
117#endif
118#if defined(_NTIGHT_)
119#define GENERATED_FILENAME "int_forward_s"
120#endif
121#define ARGUMENT(indexi,l,i) argument[indexi][l]
122#define TAYLORS(indexd,l,i)   taylors[indexd][l]
123/*--------------------------------------------------------------------------*/
124#else
125#if defined(_INDO_)
126void copy_index_domain(int res, int arg, locint **ind_dom);
127void merge_2_index_domains(int res, int arg, locint **ind_dom);
128void combine_2_index_domains(int res, int arg1, int arg2, locint **ind_dom);
129void merge_3_index_domains(int res, int arg1, int arg2, locint **ind_dom);
130
131#define NUMNNZ 20
132#define FMIN_ADOLC(x,y)  ((y<x)?y:x)
133
134#if defined(_INDOPRO_)
135#if defined(_TIGHT_)
136#define GENERATED_FILENAME "indopro_forward_t"
137#endif
138#if defined(_NTIGHT_)
139#define GENERATED_FILENAME "indopro_forward_s"
140#endif
141#endif
142#if defined(_NONLIND_)
143
144/*
145 * This is the type used for the list elements. The entry is either a counter
146 * (first element of the NID list) or the index of an independent variable.
147 */
148
149void extend_nonlinearity_domain_binary_step
150(int arg1, int arg2, locint **ind_dom, locint **nonl_dom);
151void extend_nonlinearity_domain_unary
152(int arg, locint **ind_dom, locint **nonl_dom);
153void extend_nonlinearity_domain_binary
154(int arg1, int arg2, locint **ind_dom, locint **nonl_dom);
155
156
157#if defined(_TIGHT_)
158#define GENERATED_FILENAME "nonl_ind_forward_t"
159#endif
160#if defined(_NTIGHT_)
161#define GENERATED_FILENAME "nonl_ind_forward_s"
162#endif
163#endif
164
165
166/*--------------------------------------------------------------------------*/
167#else
168#error Error ! Define [_ZOS_ | _FOS_ |\
169   _HOS_ | _FOV_ | _HOV_ | _HOV_WK_  | _INT_FOR_SAFE_ | _INT_FOR_TIGHT_ | _INDOPRO_ | _NONLIND_ ] [{_KEEP_}]
170#endif
171#endif
172#endif
173#endif
174#endif
175#endif
176#endif
177#endif
178
179/*--------------------------------------------------------------------------*/
180/*                                                               KEEP stuff */
181#if defined(_KEEP_)
182
183#if defined(_HOV_WK_) /* keep in this vector mode */
184#define IF_KEEP_TAYLOR_CLOSE \
185if (keep){\
186  fprintf(DIAG_OUT,"Succeeding reverse sweep will fail!\n");\
187  taylor_close(0);\
188}
189#define IF_KEEP_WRITE_TAYLOR(res,keep,k,p) \
190    { \
191        UPDATE_TAYLORWRITTEN(keep * k * p) \
192        if (keep) \
193        { \
194            ADOLC_WRITE_SCAYLOR(dp_T0[res]); \
195            if (keep > 1) \
196            write_taylors(res,(keep-1),k,p); \
197        } \
198    }
199#else
200#if defined(_ADOLC_VECTOR_) /* otherwise no keep */
201#define IF_KEEP_TAYLOR_CLOSE
202#define IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
203#else /* _ZOS_, _FOS_, _HOS_ */
204#define IF_KEEP_TAYLOR_CLOSE \
205if (keep){\
206  fprintf(DIAG_OUT,"Otherwise succeeding reverse sweep will fail!\n");\
207  taylor_close(0);\
208}
209#if defined(_ZOS_)
210#define IF_KEEP_WRITE_TAYLOR(res,keep,k,p) \
211    { \
212        UPDATE_TAYLORWRITTEN(keep) \
213        if (keep) \
214            ADOLC_WRITE_SCAYLOR(dp_T0[res]); \
215    }
216#else
217#if defined(_FOS_)
218#define IF_KEEP_WRITE_TAYLOR(res,keep,k,p) \
219    { \
220        UPDATE_TAYLORWRITTEN(keep) \
221        if (keep) \
222        { \
223            ADOLC_WRITE_SCAYLOR(dp_T0[res]); \
224            if (keep > 1) \
225                ADOLC_WRITE_SCAYLOR(dp_T[res]); \
226        } \
227    }
228#else
229#if defined(_HOS_)
230#define IF_KEEP_WRITE_TAYLOR(res,keep,k,p) \
231    { \
232        UPDATE_TAYLORWRITTEN(keep) \
233        if (keep) \
234        { \
235            ADOLC_WRITE_SCAYLOR(dp_T0[res]); \
236            if (keep > 1) \
237                write_taylor(res,keep-1); \
238        } \
239    }
240#endif
241#endif
242#endif
243#endif
244#endif
245
246#else  /* no _KEEP_ */
247#define IF_KEEP_TAYLOR_CLOSE
248#define IF_KEEP_WRITE_TAYLOR(res,keep,k,p)
249#endif
250
251/*--------------------------------------------------------------------------*/
252/*                                                      access to variables */
253#if !defined(_ZOS_)
254#if defined(_FOS_)
255#define TRES         *Tres
256#define TARG         *Targ
257#define TARG1        *Targ1
258#define TARG2        *Targ2
259#define TQO          *Tqo
260
261#define TRES_INC     *Tres
262#define TARG_INC     *Targ
263#define TARG1_INC    *Targ1
264#define TARG2_INC    *Targ2
265#define TQO_INC      *Tqo
266
267#define TRES_DEC     *Tres
268#define TARG_DEC     *Targ
269#define TARG1_DEC    *Targ1
270#define TARG2_DEC    *Targ2
271#define TQO_DEC      *Tqo
272
273#define TRES_FOINC   *Tres
274#define TARG_FOINC   *Targ
275#define TARG1_FOINC  *Targ1
276#define TARG2_FOINC  *Targ2
277#define TQO_FOINC    *Tqo
278
279#define TRES_FODEC   *Tres
280#define DEC_TRES_FO
281#define TARG_FODEC   *Targ
282#define TARG1_FODEC  *Targ1
283#define TARG2_FODEC  *Targ2
284#define TQO_FODEC    *Tqo
285
286#define ASSIGN_T(a,b)  a = &b;
287
288#else
289#if defined(_INT_FOR_)
290#define TRES         *Tres
291#define TARG         *Targ
292#define TARG1        *Targ1
293#define TARG2        *Targ2
294#define TQO          *Tqo
295
296#define TRES_INC     *Tres++
297#define TARG_INC     *Targ++
298#define TARG1_INC    *Targ1++
299#define TARG2_INC    *Targ2++
300#define TQO_INC      *Tqo++
301
302#define TRES_DEC     *Tres--
303#define TARG_DEC     *Targ--
304#define TARG1_DEC    *Targ1--
305#define TARG2_DEC    *Targ2--
306#define TQO_DEC      *Tqo--
307
308#define TRES_FOINC   *Tres++
309#define TARG_FOINC   *Targ++
310#define TARG1_FOINC  *Targ1++
311#define TARG2_FOINC  *Targ2++
312#define TQO_FOINC    *Tqo++
313
314#define TRES_FODEC   *Tres--
315#define TARG_FODEC   *Targ--
316#define TARG1_FODEC  *Targ1--
317#define TARG2_FODEC  *Targ2--
318#define TQO_FODEC    *Tqo--
319
320
321#define ASSIGN_T(a,b)  a = b;
322
323#else  /* _HOS_, _FOV_, _HOV_, _HOV_WK */
324#define TRES         *Tres
325#define TARG         *Targ
326#define TARG1        *Targ1
327#define TARG2        *Targ2
328#define TQO          *Tqo
329
330#define TRES_INC     *Tres++
331#define TARG_INC     *Targ++
332#define TARG1_INC    *Targ1++
333#define TARG2_INC    *Targ2++
334#define TQO_INC      *Tqo++
335
336#define TRES_DEC     *Tres--
337#define TARG_DEC     *Targ--
338#define TARG1_DEC    *Targ1--
339#define TARG2_DEC    *Targ2--
340#define TQO_DEC      *Tqo--
341
342#if defined(_FOV_)
343#define TRES_FOINC   *Tres++
344#define TARG_FOINC   *Targ++
345#define TARG1_FOINC  *Targ1++
346#define TARG2_FOINC  *Targ2++
347#define TQO_FOINC    *Tqo++
348
349#define TRES_FODEC   *Tres
350#define DEC_TRES_FO  Tres--;
351#define TARG_FODEC   *Targ--
352#define TARG1_FODEC  *Targ1--
353#define TARG2_FODEC  *Targ2--
354#define TQO_FODEC    *Tqo--
355#else /* _HOS_, _HOV_, _HOV_WK */
356#define TRES_FOINC   *Tres
357#define TARG_FOINC   *Targ
358#define TARG1_FOINC  *Targ1
359#define TARG2_FOINC  *Targ2
360#define TQO_FOINC    *Tqo
361
362#define TRES_FODEC   *Tres
363#define DEC_TRES_FO
364#define TARG_FODEC   *Targ
365#define TARG1_FODEC  *Targ1
366#define TARG2_FODEC  *Targ2
367#define TQO_FODEC    *Tqo
368#endif
369#endif
370
371#define ASSIGN_T(a,b)  a = b;
372#endif
373#endif
374
375
376/*--------------------------------------------------------------------------*/
377/*                                                               loop stuff */
378#if defined(_ADOLC_VECTOR_)
379#define FOR_0_LE_l_LT_p for (l=0; l<p; l++)
380#define FOR_p_GT_l_GE_0 for (l=p-1; l>=0; l--)
381#else
382#if defined(_INT_FOR_)
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#define FOR_0_LE_l_LT_p
387#define FOR_p_GT_l_GE_0
388#endif
389#endif
390
391#if defined(_HIGHER_ORDER_)
392#define FOR_0_LE_i_LT_k for (i=0; i<k; i++)
393#define FOR_k_GT_i_GE_0 for (i=k-1; i>=0; i--)
394#else
395#define FOR_0_LE_i_LT_k
396#define FOR_k_GT_i_GE_0
397#endif
398
399#if defined(_HOV_)
400#define FOR_0_LE_l_LT_pk for (l=0; l<pk; l++)
401#define INC_pk_1(T)      T += pk-1;
402#define VEC_INC(T,inc)   T += inc;
403#define HOV_INC(T,inc)   T += inc;
404#else
405#if defined(_HOV_WK_)
406#define FOR_0_LE_l_LT_pk for (l=0; l<pk; l++)
407#define INC_pk_1(T)      T += pk-1;
408#define VEC_INC(T,inc)   T += inc;
409#define HOV_INC(T,inc)   T += inc;
410#else
411#if defined(_FOV_)
412#define FOR_0_LE_l_LT_pk for (l=0; l<p; l++)
413#define INC_pk_1(T)      T += p-1;
414#define VEC_INC(T,inc)   T++;
415#define HOV_INC(T,inc)
416#else
417#if defined(_HOS_)
418#define FOR_0_LE_l_LT_pk for (l=0; l<k; l++)
419#define INC_pk_1(T)      T += k-1;
420#define VEC_INC(T,inc)
421#define HOV_INC(T,inc)
422#else
423#if defined(_INT_FOR_)
424#define FOR_0_LE_l_LT_pk for (l=0; l<p; l++)
425#define INC_pk_1(T)      T += p-1;
426#define VEC_INC(T,inc)   T++;
427#else
428#define FOR_0_LE_l_LT_pk
429#define INC_pk_1(T)
430#define VEC_INC(T,inc)
431#define HOV_INC(T,inc)
432#endif
433#endif
434#endif
435#endif
436#endif
437
438/*--------------------------------------------------------------------------*/
439/*                                                        higher order case */
440#if defined(_HIGHER_ORDER_)
441#define BREAK_FOR_I break;
442#else
443#define BREAK_FOR_I ;
444#endif
445
446/* END Macros */
447
448BEGIN_C_DECLS
449
450#if defined(_ZOS_)
451/****************************************************************************/
452/* Zero Order Scalar version of the forward mode.                           */
453/****************************************************************************/
454#if defined(_KEEP_)
455int  zos_forward(
456#else
457int  zos_forward_nk(
458#endif
459    short  tnum,              /* tape id */
460    int    depcheck,          /* consistency chk on # of deps */
461    int    indcheck,          /* consistency chk on # of indeps */
462#if defined(_KEEP_)
463    int    keep,              /* flag for reverse sweep */
464#endif
465    const double *basepoint,  /* independant variable values */
466    double       *valuepoint) /* dependent variable values */
467
468#else
469#if defined(_FOS_)
470/****************************************************************************/
471/* First Order Scalar version of the forward mode.                          */
472/****************************************************************************/
473#if defined(_KEEP_)
474int  fos_forward(
475#else
476int  fos_forward_nk(
477#endif
478    short  tnum,        /* tape id */
479    int    depcheck,    /* consistency chk on # of deps */
480    int    indcheck,    /* consistency chk on # of indeps */
481#if defined(_KEEP_)
482    int    keep,        /* flag for reverse sweep */
483#endif
484    const double *basepoint,  /* independent variable values */
485    double *argument,   /* Taylor coefficients (input) */
486    double *valuepoint, /* Taylor coefficients (output) */
487    double *taylors)    /* matrix of coefficient vectors */
488/* the order of the indices in argument and taylors is [var][taylor] */
489
490#else
491#if defined(_INT_FOR_)
492#if defined(_TIGHT_)
493/****************************************************************************/
494/* First Order Vector version of the forward mode for bit patterns, tight   */
495/****************************************************************************/
496int int_forward_tight(
497    short               tnum,     /* tape id                              */
498    int                 depcheck, /* consistency chk on # of dependents   */
499    int                 indcheck, /* consistency chk on # of independents */
500    int                 p,        /* # of taylor series, bit pattern      */
501    const double       *basepoint,  /* independent variable values   (in)*/
502    unsigned long int **argument,  /* Taylor coeff.                 (in)*/
503    double             *valuepoint, /* dependent variable values    (out)*/
504    unsigned long int **taylors)   /* matrix of coefficient vectors(out)*/
505
506/* int_forward_tight( tag, m, n, p, x[n], X[n][p], y[m], Y[m][p]),
507
508     nBV = number of Boolean Vectors to be packed
509                      (see Chapter Dependence Analysis, ADOL-C Documentation)
510     bits_per_long = 8*sizeof(unsigned long int)
511     p = nBV / bits_per_long + ( (nBV % bits_per_long) != 0 )
512
513     The order of the indices in argument and taylors is [var][taylor]
514
515     For the full Jacobian matrix set
516     p = indep / bits_per_long + ((indep % bits_per_long) != 0)
517     and pass a bit pattern version of the identity matrix as an argument   */
518
519
520#endif
521#if defined (_NTIGHT_)
522/****************************************************************************/
523/* First Order Vector version of the forward mode, bit pattern, safe        */
524/****************************************************************************/
525int int_forward_safe(
526    short             tnum,     /* tape id                              */
527    int               depcheck, /* consistency chk on # of dependents   */
528    int               indcheck, /* consistency chk on # of independents */
529    int               p,        /* # of taylor series, bit pattern      */
530    unsigned long int **argument, /* Taylor coeff.                  (in)*/
531    unsigned long int **taylors)  /* matrix of coefficient vectors (out)*/
532
533/* int_forward_safe( tag, m, n, p, X[n][p], Y[m][p]),
534
535nBV = number of Boolean Vectors to be packed
536(see Chapter Dependence Analysis, ADOL-C Documentation)
537bits_per_long = 8*sizeof(unsigned long int)
538p = nBV / bits_per_long + ( (nBV % bits_per_long) != 0 )
539
540The order of the indices in argument and taylors is [var][taylor]
541
542For the full Jacobian matrix set
543p = indep / bits_per_long + ((indep % bits_per_long) != 0)
544and pass a bit pattern version of the identity matrix as an argument    */
545#endif
546#else
547#if defined(_INDOPRO_)
548#if defined(_TIGHT_)
549/****************************************************************************/
550/* First Order Vector version of the forward mode for bit patterns, tight   */
551/****************************************************************************/
552int indopro_forward_tight(
553    short             tnum,        /* tape id                              */
554    int               depcheck,    /* consistency chk on # of dependents   */
555    int               indcheck,    /* consistency chk on # of independents */
556    const double     *basepoint,  /* independent variable values   (in)   */
557    unsigned int    **crs)        /* returned row index storage (out)     */
558
559/* indopro_forward_tight( tag, m, n, x[n], *crs[m]),
560
561  */
562
563
564#endif
565#if defined (_NTIGHT_)
566/****************************************************************************/
567/* First Order Vector version of the forward mode, bit pattern, safe        */
568/****************************************************************************/
569int indopro_forward_safe(
570    short             tnum,        /* tape id                              */
571    int               depcheck,    /* consistency chk on # of dependents   */
572    int               indcheck,    /* consistency chk on # of independents */
573    const double     *basepoint,   /* independent variable values   (in)   */
574    unsigned int    **crs)         /* returned row index storage (out)     */
575
576/* indopro_forward_safe( tag, m, n, x[n], *crs[m]),
577
578  */
579#endif
580#else
581#if defined(_NONLIND_)
582#if defined(_TIGHT_)
583/****************************************************************************/
584/* First Order Vector version of the forward mode for bit patterns, tight   */
585/****************************************************************************/
586int nonl_ind_forward_tight(
587    short             tnum,        /* tape id                              */
588    int               depcheck,    /* consistency chk on # of dependents   */
589    int               indcheck,    /* consistency chk on # of independents */
590    const double     *basepoint,  /* independent variable values   (in)   */
591    unsigned int     **crs)        /* returned row index storage (out)     */
592
593/* indopro_forward_tight( tag, m, n, x[n], *crs[m]),
594
595  */
596
597#endif
598#if defined (_NTIGHT_)
599/****************************************************************************/
600/* First Order Vector version of the forward mode, bit pattern, safe        */
601/****************************************************************************/
602int nonl_ind_forward_safe(
603    short             tnum,        /* tape id                              */
604    int               depcheck,    /* consistency chk on # of dependents   */
605    int               indcheck,    /* consistency chk on # of independents */
606    const double      *basepoint,  /* independent variable values   (in)   */
607    unsigned int    **crs)        /* returned row index storage (out)     */
608
609/* indopro_forward_safe( tag, m, n, x[n], *crs[m]),
610
611  */
612#endif
613#else
614#if defined(_FOV_)
615#if defined(_CHUNKED_)
616/****************************************************************************/
617/* First Order Vector version of the forward mode with p-offset in          */
618/* **argument and **taylors                                                 */
619/****************************************************************************/
620int  fov_offset_forward(
621    short  tnum,        /* tape id */
622    int    depcheck,    /* consistency chk on # of deps */
623    int    indcheck,    /* consistency chk on # of indeps */
624    int    p,           /* # of taylor series */
625    int    offset,      /* offset for assignments */
626    double *basepoint,  /* independent variable values */
627    double **argument,  /* Taylor coefficients (input) */
628    double *valuepoint, /* Taylor coefficients (output) */
629    double **taylors)   /* matrix of coifficient vectors */
630/* the order of the indices in argument and taylors is [var][taylor] */
631#else
632/****************************************************************************/
633/* First Order Vector version of the forward mode.                          */
634/****************************************************************************/
635int  fov_forward(
636    short         tnum,        /* tape id */
637    int           depcheck,    /* consistency chk on # of deps */
638    int           indcheck,    /* consistency chk on # of indeps */
639    int           p,           /* # of taylor series */
640    const double *basepoint,   /* independent variable values */
641    double      **argument,    /* Taylor coefficients (input) */
642    double       *valuepoint,  /* Taylor coefficients (output) */
643    double      **taylors)     /* matrix of coifficient vectors */
644/* the order of the indices in argument and taylors is [var][taylor] */
645#endif
646
647#else
648#if defined(_HOS_)
649/****************************************************************************/
650/* Higher Order Scalar version of the forward mode.                         */
651/****************************************************************************/
652#if defined(_KEEP_)
653int  hos_forward(
654#else
655int  hos_forward_nk(
656#endif
657    short  tnum,        /* tape id */
658    int    depcheck,    /* consistency chk on # of dependents */
659    int    indcheck,    /* consistency chk on # of independents */
660    int    gdegree,     /* highest derivative degree */
661#if defined(_KEEP_)
662    int    keep,        /* flag for reverse sweep */
663#endif
664    double *basepoint,  /* independent variable values */
665    double **argument,  /* independant variable values */
666    double *valuepoint, /* Taylor coefficients (output) */
667    double **taylors)   /* matrix of coifficient vectors */
668
669
670#else
671/****************************************************************************/
672/* Higher Order Vector version of the forward mode.                         */
673/****************************************************************************/
674#if defined(_KEEP_)
675int  hov_wk_forward(
676#else
677int  hov_forward(
678#endif
679    short  tnum,        /* tape id */
680    int    depcheck,    /* consistency chk on # of deps */
681    int    indcheck,    /* consistency chk on # of indeps */
682    int    gdegree,     /* highest derivative degree */
683#if defined(_KEEP_)
684    int    keep,        /* flag for reverse sweep */
685#endif
686    int    p,           /* # of taylor series */
687    double *basepoint,  /* independent variable values */
688    double ***argument, /* Taylor coefficients (input) */
689    double *valuepoint, /* Taylor coefficients (output) */
690    double ***taylors)  /* matrix of coifficient vectors */
691/* the order of the indices in argument and taylors is [var][taylor][deriv] */
692
693#endif
694#endif
695#endif
696#endif
697#endif
698#endif
699#endif
700{
701    /****************************************************************************/
702    /*                                                            ALL VARIABLES */
703    unsigned char operation;   /* operation code */
704    int ret_c =3;              /* return value */
705
706    locint size = 0;
707    locint res  = 0;
708    locint arg  = 0;
709    locint arg1 = 0;
710    locint arg2 = 0;
711
712    double coval = 0, *d = 0;
713
714    int indexi = 0,  indexd = 0;
715
716    /* loop indices */
717#if !defined (_ZOS_)
718#if !defined (_INT_FOR_)
719    int i;
720#if !defined (_INDO_)
721    int ii;
722#endif
723#endif
724#endif
725#if defined (_HIGHER_ORDER_)
726    int j, l=0;
727#endif
728    int ls;
729#if defined(_ADOLC_VECTOR_)
730#if !defined (_HIGHER_ORDER_)
731    int l=0;
732#endif
733#endif
734#if defined (_INT_FOR_)
735    int l=0;
736#endif
737#if defined (_INDO_)
738    int l=0;
739    int max_ind_dom;
740#if defined(_NONLIND_)
741    /* nonlinear interaction domains */
742    locint** nonl_dom;
743    locint*  temp;
744    locint*  temp1;
745#endif
746#endif
747
748    /* other necessary variables */
749#if !defined (_ZOS_)
750#if !defined (_INDO_)
751#if !defined (_INT_FOR_)
752    double r0=0.0, x, y, divs;
753    int even;
754#endif
755#endif
756#endif
757
758#if defined(_INT_FOR_)
759#ifdef _TIGHT_
760    double  *dp_T0;
761    double y, divs;
762#endif /* _TIGHT_ */
763
764    /* Taylor stuff */
765    unsigned long int  **up_T;
766
767    unsigned long int         *Tres, *Targ, *Targ1, *Targ2;
768#ifdef _TIGHT_
769    unsigned long int         *Tqo;
770    unsigned long int         *Targ1OP, *Targ2OP;
771#endif
772
773#define T0res  T0temp
774#else
775#if defined(_INDO_)
776#ifdef _TIGHT_
777    double  *dp_T0;
778    double  T0temp;
779    double divs;
780#endif /* _TIGHT_ */
781#define T0res  T0temp
782#define T0arg  T0temp
783
784    /* index domains */
785    locint** ind_dom;
786
787#else
788    double *dp_T0;
789#if !defined(_ZOS_)
790#if  defined(_FOS_)
791    double  *dp_T;
792# define T_TEMP Ttemp;
793# else
794    double *dp_Ttemp, **dpp_T;
795#endif
796    double         *Tres, *Targ, *Targ1, *Targ2, *Tqo;
797
798#if defined (_HIGHER_ORDER_)
799    double         *TresOP, *TresOP2, *zOP;
800    double *dp_z;
801#endif
802   double         *TargOP, *Targ1OP, *Targ2OP;
803   double         T0temp;
804#endif
805#define T0res  T0temp
806#define T0arg  T0temp
807#endif
808#endif
809
810#if defined(_HIGHER_ORDER_)
811    int k = gdegree;
812#endif
813
814#if defined(_KEEP_)
815    int taylbuf=0;
816#endif
817
818#if defined(_HOV_)
819    int pk = k*p;
820#else
821#if defined(_HOV_WK_)
822    int pk = k*p;
823#endif
824#endif
825
826    /* extern diff. function variables */
827#if defined(_EXTERN_)
828#  undef (_EXTERN_)
829#endif
830    /* ZOS_FORWARD */
831#if defined(_ZOS_)
832#   define _EXTERN_ 1
833#   define ADOLC_EXT_FCT_POINTER zos_forward
834#   define ADOLC_EXT_FCT_COMPLETE \
835    zos_forward(n, edfct->dp_x, m, edfct->dp_y)
836#   define ADOLC_EXT_LOOP
837#   define ADOLC_EXT_SUBSCRIPT
838#endif
839    /* FOS_FORWARD */
840#if defined(_FOS_)
841#   define _EXTERN_ 1
842#   define ADOLC_EXT_FCT_POINTER fos_forward
843#   define ADOLC_EXT_FCT_COMPLETE \
844    fos_forward(n, edfct->dp_x, edfct->dp_X, m, edfct->dp_y, edfct->dp_Y)
845#   define ADOLC_EXT_POINTER_X edfct->dp_X
846#   define ADOLC_EXT_POINTER_Y edfct->dp_Y
847#   define ADOLC_EXT_LOOP
848#   define ADOLC_EXT_SUBSCRIPT
849#endif
850    /* FOV_FORWARD */
851#if defined(_FOV_)
852#   define _EXTERN_ 1
853#   define ADOLC_EXT_FCT_POINTER fov_forward
854#   define ADOLC_EXT_FCT_COMPLETE \
855    fov_forward(n, edfct->dp_x, edfct->dpp_X, m, edfct->dp_y, edfct->dpp_Y)
856#   define ADOLC_EXT_POINTER_X edfct->dpp_X
857#   define ADOLC_EXT_POINTER_Y edfct->dpp_Y
858#   define ADOLC_EXT_LOOP for (loop2 = 0; loop2 < p; ++loop2)
859#   define ADOLC_EXT_SUBSCRIPT [loop2]
860#endif
861
862#if defined(_EXTERN_)
863    locint n, m;
864    ext_diff_fct *edfct;
865    int loop;
866#   if defined(_FOV_)
867        int loop2;
868#   endif
869    int ext_retc;
870#endif
871
872    ADOLC_OPENMP_THREAD_NUMBER;
873
874#if defined(ADOLC_DEBUG)
875    /****************************************************************************/
876    /*                                                           DEBUG MESSAGES */
877    fprintf(DIAG_OUT,"Call of %s(..) with tag: %d, n: %d, m %d,\n",
878            GENERATED_FILENAME, tnum, indcheck, depcheck);
879#if defined(_KEEP_)
880    fprintf(DIAG_OUT,"                    keep: %d\n", keep);
881#endif
882#if defined(_HIGHER_ORDER_)
883    fprintf(DIAG_OUT,"                    degree: %d\n",gdegree);
884#endif
885#if defined(_ADOLC_VECTOR_)
886    fprintf(DIAG_OUT,"                    p: %d\n\n",p);
887#endif
888
889#endif
890
891    /****************************************************************************/
892    /*                                                                    INITs */
893
894   /* Set up stuff for the tape */
895    ADOLC_OPENMP_GET_THREAD_NUMBER;
896
897    /* Initialize the Forward Sweep */
898
899    init_for_sweep(tnum);
900
901      if ( (depcheck != ADOLC_CURRENT_TAPE_INFOS.stats[NUM_DEPENDENTS]) ||
902            (indcheck != ADOLC_CURRENT_TAPE_INFOS.stats[NUM_INDEPENDENTS]) ) {
903        fprintf(DIAG_OUT,"ADOL-C error: forward sweep on tape %d  aborted!\n"
904                "Number of dependent(%u) and/or independent(%u) variables passed"
905                " to forward is\ninconsistent with number "
906                "recorded on tape (%d, %d) \n", tnum,
907                depcheck, indcheck,
908                ADOLC_CURRENT_TAPE_INFOS.stats[NUM_DEPENDENTS],
909                ADOLC_CURRENT_TAPE_INFOS.stats[NUM_INDEPENDENTS]);
910        exit (-1);
911    }
912
913
914    /****************************************************************************/
915    /*                                                        MEMORY ALLOCATION */
916    /* olvo 980626 has to be revised for common blocks */
917
918    /*--------------------------------------------------------------------------*/
919#if !defined(_NTIGHT_)
920    dp_T0 = myalloc1(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES]);
921    ADOLC_CURRENT_TAPE_INFOS.dp_T0 = dp_T0;
922#endif /* !_NTIGHT_ */
923#if defined(_ZOS_)                                                   /* ZOS */
924
925#if defined(_KEEP_)
926    if (keep>1) {
927        fprintf(DIAG_OUT,"\n ADOL-C error: zero order scalar forward cannot save"
928                " more\nthan zero order taylor coefficients!\n");
929        exit (-1);
930    }
931#endif
932#if defined(_KEEP_)
933    if (keep) {
934      taylbuf = ADOLC_CURRENT_TAPE_INFOS.stats[TAY_BUFFER_SIZE];
935
936        taylor_begin(taylbuf,&dp_T0,keep-1);
937    }
938#endif
939
940    /*--------------------------------------------------------------------------*/
941#else                                                                /* FOS */
942#if defined(_FOS_)
943#if defined(_KEEP_)
944    if (keep>2) {
945        fprintf(DIAG_OUT,"\n ADOL-C error: first order scalar forward cannot save"
946                " more  \nthan first order taylor coefficients!\n");
947        exit (-1);
948    }
949#endif
950    dp_T = myalloc1(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES]);
951# define TAYLOR_BUFFER dp_T
952#if defined(_KEEP_)
953    if (keep) {
954        taylbuf = ADOLC_CURRENT_TAPE_INFOS.stats[TAY_BUFFER_SIZE];
955        taylor_begin(taylbuf,&dp_T,keep-1);
956    }
957#endif
958
959    /*--------------------------------------------------------------------------*/
960#else                                                                /* INF_FOR */
961#if defined(_INT_FOR_)
962        up_T     = myalloc2_ulong(ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES],p);
963#define TAYLOR_BUFFER up_T
964
965    /*--------------------------------------------------------------------------*/
966#else                                                                /* INDOPRO */
967#if defined(_INDO_)
968    ind_dom = (locint **)  malloc(sizeof(locint*) * ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES]);
969
970    for(i=0;i<ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES];i++)
971    {
972        ind_dom[i] = (locint *)  malloc(sizeof(locint) * (NUMNNZ+2));
973        ind_dom[i][0] = 0;
974        ind_dom[i][1] = NUMNNZ;
975    }
976
977    max_ind_dom = ADOLC_CURRENT_TAPE_INFOS.stats[NUM_MAX_LIVES];
978#if defined(_NONLIND_)
979
980    nonl_dom = (locint**) malloc(sizeof(locint*) * indcheck);
981    for(i=0;i<indcheck;i++){
982          nonl_dom[i] = (locint*) malloc(sizeof(locint)*(NUMNNZ+2));
983          nonl_dom[i][0]=0;
984          nonl_dom[i][1]=NUMNNZ;
985       }
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#if defined(ADOLC_DEBUG)
3660    printf("\nTape contains:\n");
3661    for (v = 0; v < 256; ++v)
3662        if (countPerOperation[v] > 0)
3663            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]);
3664    printf("\n");
3665#endif /* ADOLC_DEBUG */
3666
3667#if defined(_KEEP_)
3668    if (keep) taylor_close(taylbuf);
3669#endif
3670
3671    /* clean up */
3672#if !defined (_NTIGHT_)
3673    free(dp_T0);
3674#endif /* !_NTIGHT_ */
3675#if !defined(_INDO_)
3676#if !defined(_ZOS_)
3677#   if defined(_FOS_)
3678    free(dp_T);
3679#   else
3680#if !defined (_INT_FOR_)
3681    myfree2(dpp_T);
3682    free(dp_Ttemp);
3683#endif /* !_NTIGHT_ */
3684#endif
3685#endif
3686#endif
3687#if defined(_HIGHER_ORDER_)
3688    free(dp_z);
3689#endif
3690
3691    end_sweep();
3692
3693
3694#if defined(_INDO_)
3695
3696    for(i=0;i<max_ind_dom;i++)
3697      {
3698        free(ind_dom[i]);
3699      }
3700    free(ind_dom);
3701
3702#if defined(_NONLIND_)
3703
3704    for( i=0; i < indcheck; i++) {
3705       crs[i] = (unsigned int*) malloc(sizeof(unsigned int) * (nonl_dom[i][0]+1));
3706       crs[i][0] = nonl_dom[i][0];
3707       for(l=1; l < crs[i][0]+1; l++)
3708          crs[i][l] = nonl_dom[i][l+1];
3709       free(nonl_dom[i]);
3710    }
3711    free(nonl_dom);
3712
3713#endif
3714#endif
3715    return ret_c;
3716}
3717
3718/****************************************************************************/
3719
3720#if defined(_INDOPRO_)
3721
3722/****************************************************************************/
3723/* set operations for propagation of index domains                          */
3724
3725/*--------------------------------------------------------------------------*/
3726/* operations on index domains                                              */
3727
3728#if defined(_TIGHT_)
3729void copy_index_domain(int res, int arg, locint **ind_dom) {
3730
3731   int i;
3732
3733   if (ind_dom[arg][0] > ind_dom[res][1]-2)
3734     {
3735       free(ind_dom[res]);
3736       ind_dom[res] = (locint *)  malloc(sizeof(locint) * 2*(ind_dom[arg][0])+1);
3737       ind_dom[res][1] = 2*ind_dom[arg][0];
3738     }
3739
3740
3741    for(i=2;i<ind_dom[arg][0]+2;i++)
3742       ind_dom[res][i] = ind_dom[arg][i];
3743    ind_dom[res][0] = ind_dom[arg][0];
3744}
3745
3746
3747void merge_2_index_domains(int res, int arg, locint **ind_dom) 
3748{
3749
3750  int num,num1,num2, i,j,k,l;
3751  locint *temp_array, *arg_ind_dom, *res_ind_dom;
3752
3753  if (ind_dom[res][0] == 0)
3754    copy_index_domain(res,arg,ind_dom);
3755  else
3756    {
3757      if (res != arg)
3758        {
3759          arg_ind_dom = ind_dom[arg];
3760          res_ind_dom = ind_dom[res];
3761
3762          num  = ind_dom[res][0];
3763          num1 = arg_ind_dom[0];
3764          num2 = ind_dom[res][1];
3765
3766          if (num2 < num1+num)
3767            num2 = num1+num;
3768         
3769          temp_array = (locint *)  malloc(sizeof(locint)* (num2+2));
3770          temp_array[1] = num2;
3771
3772          i = 2;
3773          j = 2;
3774          k = 2;
3775          num += 2;
3776          num1 += 2;
3777          while ((i< num) && (j < num1))
3778            {
3779              if (res_ind_dom[i] < arg_ind_dom[j])
3780                {
3781                  temp_array[k] = res_ind_dom[i];
3782                  i++; k++;
3783                }
3784              else
3785                {
3786                  if (res_ind_dom[i] == arg_ind_dom[j])
3787                    {
3788                      temp_array[k] = arg_ind_dom[j];
3789                      i++;j++;k++;
3790                    }
3791                  else
3792                    {
3793                      temp_array[k] = arg_ind_dom[j];
3794                      j++;k++;               
3795                    }
3796                }
3797            }
3798          for(l = i;l<num;l++)
3799            {
3800              temp_array[k] = res_ind_dom[l];
3801              k++;
3802            }
3803          for(l = j;l<num1;l++)
3804            {
3805              temp_array[k] = arg_ind_dom[l];
3806              k++;
3807            }
3808          temp_array[0] = k-2;
3809          free(ind_dom[res]);
3810          ind_dom[res]=temp_array;
3811        }
3812    }
3813
3814
3815}
3816
3817void combine_2_index_domains(int res, int arg1, int arg2, locint **ind_dom) {
3818
3819    if (res != arg1)
3820       copy_index_domain(res, arg1, ind_dom);
3821
3822    merge_2_index_domains(res, arg2, ind_dom);
3823}
3824
3825void merge_3_index_domains(int res, int arg1, int arg2, locint **ind_dom) {
3826    merge_2_index_domains(res, arg1, ind_dom);
3827    merge_2_index_domains(res, arg2, ind_dom);
3828}
3829
3830
3831
3832#endif
3833#endif
3834
3835
3836#if defined(_NONLIND_)
3837#if defined(_TIGHT_)
3838
3839void extend_nonlinearity_domain_binary_step
3840(int arg1, int arg2, locint **ind_dom, locint **nonl_dom) 
3841{
3842  int index,num,num1, num2, i,j,k,l,m;
3843  locint *temp_nonl, *index_nonl_dom, *arg1_ind_dom, *arg2_ind_dom;
3844
3845  num = ind_dom[arg2][0];
3846
3847  for(m=2;m<ind_dom[arg1][0]+2;m++) 
3848    {
3849      index = ind_dom[arg1][m];
3850      index_nonl_dom = nonl_dom[index];
3851
3852      if (index_nonl_dom[0] == 0)  /* empty list */
3853        {
3854          if ( index_nonl_dom[1] < num)
3855            {
3856              free(index_nonl_dom);
3857              index_nonl_dom = (locint*) malloc(sizeof(locint)*2*(num+1) );
3858              index_nonl_dom[1] = 2*num;
3859            }
3860          for(i=2;i<num+2;i++)      /* append index domain list of "arg" */
3861            index_nonl_dom[i] = ind_dom[arg2][i];
3862          index_nonl_dom[0] = num;
3863        } 
3864      else 
3865        { /* merge lists */
3866          num1 = index_nonl_dom[0];
3867          num2 = index_nonl_dom[1];
3868         
3869          if (num1+num > num2)
3870            num2 = num1+num;
3871         
3872          temp_nonl = (locint*) malloc(sizeof(locint)*(num2+2));
3873          temp_nonl[1] = num2;
3874         
3875          i = 2;
3876          k = 2;
3877          j = 2;
3878          num1 +=2;
3879          num2 = num+2;
3880          while ((i<num1) && (j < num2)){
3881            if (ind_dom[arg2][j] < index_nonl_dom[i]) /* < */ {
3882              temp_nonl[k] = ind_dom[arg2][j];
3883              j++; k++;
3884            } else {
3885              if (ind_dom[arg2][j] == index_nonl_dom[i])  /* == */ {
3886                temp_nonl[k] = ind_dom[arg2][j];
3887                j++; k++; i++;
3888              } else {
3889                temp_nonl[k] = index_nonl_dom[i];
3890                i++; k++;
3891              }
3892            }
3893          }
3894          for(l = j;l<num2;l++) {
3895            temp_nonl[k] = ind_dom[arg2][l];
3896            k++;
3897          }
3898          for(l = i;l<num1;l++) {
3899            temp_nonl[k] = index_nonl_dom[l];
3900            k++;
3901          }
3902          temp_nonl[0] = k-2; 
3903          free((char*) nonl_dom[index]);
3904          nonl_dom[index] = temp_nonl;
3905        }
3906    }
3907}
3908
3909void extend_nonlinearity_domain_unary
3910(int arg, locint **ind_dom, locint **nonl_dom) {
3911    extend_nonlinearity_domain_binary_step(arg, arg, ind_dom, nonl_dom);
3912}
3913
3914void extend_nonlinearity_domain_binary
3915(int arg1, int arg2, locint **ind_dom, locint **nonl_dom) {
3916    extend_nonlinearity_domain_binary_step(arg1, arg2, ind_dom, nonl_dom);
3917    extend_nonlinearity_domain_binary_step(arg2, arg1, ind_dom, nonl_dom);
3918}
3919
3920
3921#endif
3922#endif
3923END_C_DECLS
Note: See TracBrowser for help on using the repository browser.