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

Last change on this file since 762 was 762, checked in by mbanovic, 8 months ago

Merged branch "medipacksupport" from "git" into "svn"

The following commits were merged:

commit 0d1b5eec2cca8afdeea3cdffa196efb6cfd60a53
Merge: 72d114b 33bfdb5
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Mon Nov 5 10:03:04 2018 +0000

Merge branch 'medipackSupport' into 'medipacksupport'

Medipack support

See merge request adol-c/adol-c!26

commit 33bfdb5a006c782489bfef1b651ca3bdbceefaf2
Merge: ac55eab cf82982
Author: Max Sagebaum <max.sagebaum@…>
Date: Tue Oct 30 11:19:31 2018 +0100

Merge branch 'medipackSupport' into temp

commit ac55eab9dd8cb8c84926ee56456076392a047c1a
Merge: 72d114b caaac60
Author: Max Sagebaum <max.sagebaum@…>
Date: Tue Oct 30 11:14:09 2018 +0100

Merge remote-tracking branch 'origin/master' into temp

commit cf82982421aaa7d83405ffa3d0c9b6ef88251d0c
Merge: 6aeca20 caaac60
Author: Max Sagebaum <max.sagebaum@…>
Date: Tue Oct 30 11:13:25 2018 +0100

Merge remote-tracking branch 'origin/master' into medipackSupport

commit 6aeca205c2448b4bbc915eb76153ebde19448573
Author: Max Sagebaum <max.sagebaum@…>
Date: Tue Oct 23 22:30:28 2018 +0200

Added suport for ZOS, FOS, FOV forward and reverse.

commit caaac60da4c61b370d106c68064d38c42a7cb6e3
Merge: cc2e0b3 70fc288
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Mon Oct 8 08:53:40 2018 +0000

Merge branch 'fix_adtl_hov_refcntr' into 'master'

Fix undefined reference to adtl_hov::refcounter::refcnt

See merge request adol-c/adol-c!24

commit 70fc288b9ab95b16d3179fcf239ee2208ae1a2c4
Author: Jean-Paul Pelteret <jppelteret@…>
Date: Mon Oct 1 20:53:03 2018 +0200

Fix undefined reference to adtl_hov::refcounter::refcnt

commit cc2e0b3154fb6e62580def4501c4cf3f3d8e32ef
Merge: d7400f5 7c7f24b
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Mon Oct 1 12:26:39 2018 +0000

Merge branch 'docu' into 'master'

Refactor tapeless to traceless

See merge request adol-c/adol-c!23

commit ca397962cde23bde80e03924893e09c84d8728bf
Merge: 9cbc432 d7400f5
Author: Max Sagebaum <max.sagebaum@…>
Date: Fri Sep 28 10:07:41 2018 +0200

Merge remote-tracking branch 'origin/master' into medipackSupport

commit 9cbc4324e0d3e19f97ba5c5474121f0189e60f83
Author: Max Sagebaum <max.sagebaum@…>
Date: Thu Sep 27 14:38:30 2018 +0200

Missing MeDiPack? initialization on trace_on.

commit 76c30290365830d2370a354af949f3bf42df3885
Author: Max Sagebaum <max.sagebaum@…>
Date: Thu Sep 27 09:55:42 2018 +0200

Null pointer fix for initialization.

commit 7c7f24b25479870d58ff19d78a6e394ca28ddb58
Author: mflehmig <martin.schroschk@…>
Date: Thu Sep 20 13:16:06 2018 +0200

Refactor tapeless to traceless

As far as I can see, the official wording is traceless forward mode.
Additonally, the latex label and refs changed to 'traceless'.

commit 72d114b7ac42b8ac493030cedd1df8c9746ba5d4
Author: Max Sagebaum <max.sagebaum@…>
Date: Thu Oct 19 09:25:19 2017 +0200

Added support for MeDiPack? library.

Enable it with the configure options:
--enable-medipack --with-medipack=<path to MeDiPack?>

Tutorial on a how to use will follow.

commit b4ca76279d28407f29901d40953d02a0c5c9140e
Author: Max Sagebaum <max.sagebaum@…>
Date: Mon May 7 14:45:13 2018 +0200

Added support for cbrt function.

commit bc7b8ca61865058fac097410fd94a44fba281131
Author: Max Sagebaum <max.sagebaum@…>
Date: Thu Mar 1 10:31:18 2018 +0100

Changes for new interface.

commit cd1e82778c8540221b24559d5097bf6d00597e19
Author: Max Sagebaum <max.sagebaum@…>
Date: Thu Nov 16 14:31:07 2017 +0100

Changes to new MeDiPack? interface for adjoint values.

commit 55bcb0ffd5a9496817bffac0bd2c9489ed8ce992
Author: Max Sagebaum <max.sagebaum@…>
Date: Thu Oct 19 09:25:19 2017 +0200

Added support for MeDiPack? library.

Enable it with the configure options:
--enable-medipack --with-medipack=<path to MeDiPack?>

Tutorial on a how to use will follow.

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