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

Last change on this file since 42 was 42, checked in by awalther, 11 years ago

set svn keywords property

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