a571de7bb16226fd8f076ed626d0c4ebbbef71bc
[scilab.git] / scilab / modules / scicos / src / c / scicos.c
1 /*  Scicos
2 *
3 *  Copyright (C) INRIA - METALAU Project <scicos@inria.fr>
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 *
19 * See the file ./license.txt
20 */
21 /* 11-03-2005,  Masoud
22 * adding A-Jacobian
23 * istate =-1 case;
24 *
25 * xx/03/06, Alan
26 * enable matrix typed
27 * input/output regular ports
28 *
29 * xx/02/07, Alan
30 * added object paramaters/states
31 *
32 * 20/07/07,  Masoud
33 * CVODE (Sundials) replaced LSODAR
34 * IDA  (Sundials)  replaced DASKR
35 */
36 /*--------------------------------------------------------------------------*/
37 #include <stdlib.h>
38 #include <string.h>
39 #include <wchar.h>
40 #include <math.h>
41
42 /* Sundials includes */
43 #include <cvode/cvode.h>            /* prototypes for CVODES fcts. and consts. */
44 #include <cvode/cvode_dense.h>     /* prototype for CVDense */
45 #include <cvode/cvode_direct.h>    /* prototypes for various DlsMat operations */
46 #include <ida/ida.h>
47 #include <ida/ida_dense.h>
48 #include <ida/ida_direct.h>
49 #include <nvector/nvector_serial.h>   /* serial N_Vector types, fcts., and macros */
50 #include <sundials/sundials_dense.h>  /* prototypes for various DlsMat operations */
51 #include <sundials/sundials_direct.h> /* definitions of DlsMat and DENSE_ELEM */
52 #include <sundials/sundials_types.h>  /* definition of type realtype */
53 #include <sundials/sundials_math.h>
54 #include <kinsol/kinsol.h>
55 #include <kinsol/kinsol_dense.h>
56 #include <kinsol/kinsol_direct.h>
57 #include <sundials/sundials_extension.h> /* uses extension for scicos */
58 #include "ida_impl.h"
59
60 #include "machine.h" /* C2F */
61 #include "scicos-def.h"
62 #include "sciprint.h"
63 #include "scicos.h"
64 #include "import.h"
65 #include "scicos_internal.h"
66 #include "blocks.h"
67 #include "core_math.h"
68 #include "storeCommand.h"
69 #include "syncexec.h"
70 #include "realtime.h"
71 #include "sci_malloc.h"
72 #include "ezxml.h"
73
74 #include "sciblk2.h"
75 #include "sciblk4.h"
76 #include "dynlib_scicos.h"
77
78 #include "configvariable_interface.h" /* getEntryPointPosition() and getEntryPointFromPosition() */
79
80 #include "lsodar.h"           /* prototypes for lsodar fcts. and consts. */
81 #include "ddaskr.h"           /* prototypes for ddaskr fcts. and consts. */
82
83 #if defined(linux) && defined(__i386__)
84 #include "setPrecisionFPU.h"
85 #endif
86
87 #include "localization.h"
88 #include "charEncoding.h"
89 #include "common_structure.h"
90 #include "Sciwarning.h"
91
92 /*--------------------------------------------------------------------------*/
93 typedef struct
94 {
95     void *dae_mem;
96     N_Vector ewt;
97     double *rwork;
98     int *iwork;
99     double *gwork; /* just added for a very special use: a
100                                    space passing to grblkdakr for zero crossing surfaces
101                                    when updating mode variables during initialization */
102 } *UserData;
103
104 SCICOS_IMPEXP SCSPTR_struct C2F(scsptr);
105
106 enum Solver
107 {
108     LSodar_Dynamic = 0,
109     CVode_BDF_Newton,
110     CVode_BDF_Functional,
111     CVode_Adams_Newton,
112     CVode_Adams_Functional,
113     Dormand_Prince,
114     Runge_Kutta,
115     Implicit_Runge_Kutta,
116     IDA_BDF_Newton = 100,
117     DDaskr_BDF_Newton = 101,
118     DDaskr_BDF_GMRes = 102
119 };
120 /*--------------------------------------------------------------------------*/
121
122 #define freeall                                 \
123         if (*neq>0) ODEFree(&ode_mem);          \
124         if (*neq>0) N_VDestroy_Serial(y);               \
125         if ( ng>0 ) FREE(jroot);                        \
126         if ( ng>0 ) FREE(zcros);
127
128
129 /* TJacque allocated by sundials */
130 #define freeallx                                \
131         if (*neq>0) free(TJacque);      \
132         if (*neq>0) FREE(data->rwork);          \
133         if (( ng>0 )&& (*neq>0)) FREE(data->gwork);     \
134         if (*neq>0) N_VDestroy_Serial(data->ewt);       \
135         if (*neq>0) FREE(data);                 \
136         if (*neq>0) DAEFree(&dae_mem);  \
137         if (*neq>0) N_VDestroy_Serial(IDx);             \
138         if (*neq>0) N_VDestroy_Serial(yp);              \
139         if (*neq>0) N_VDestroy_Serial(yy);              \
140         if ( ng>0 ) FREE(jroot);                        \
141         if ( ng>0 ) FREE(zcros);                        \
142         if (nmod>0) FREE(Mode_save);
143
144 #define freeouttbptr                            \
145         FREE(outtbd);                                   \
146         FREE(outtbc);                                   \
147         FREE(outtbs);                                   \
148         FREE(outtbl);                                   \
149         FREE(outtbuc);                          \
150         FREE(outtbus);                          \
151         FREE(outtbul);
152
153 #define freekinsol                              \
154         FREE(Mode_save);                                \
155         N_VDestroy_Serial(y);                           \
156         N_VDestroy_Serial(fscale);                      \
157         N_VDestroy_Serial(yscale);                      \
158         KINFree(&kin_mem);
159
160
161 #define ONE   RCONST(1.0)
162 #define ZERO  RCONST(0.0)
163 #define T0    RCONST(0.0)
164
165 // Special values for elements of 'funtyp'
166 #define FORTRAN_GATEWAY   0
167 #define UNUSED1           1
168 #define UNUSED2           2
169 #define SCIFUNC_BLOCK     3
170 #define EXPRESSION_BLOCK  5
171 #define IFTHEL_BLOCK      11
172 #define ESELECT_BLOCK     12
173 #define DEBUG_BLOCK       99
174 #define OLD_SCI_BLOCK     10005
175 /*--------------------------------------------------------------------------*/
176 /* Table of constant values */
177 static int c__90 = 90;
178 static int c__91 = 91;
179 static int c__0 = 0;
180 static double c_b14 = 0.;
181 static int c__1 = 1;
182
183 int TCritWarning = 0;
184
185 static double *t0 = NULL, *tf = NULL;
186 static double *x = NULL, *tevts = NULL, *xd = NULL, *g = NULL;
187 static  double Atol = 0., rtol = 0., ttol = 0., deltat = 0., hmax = 0.;
188 static int *ierr = NULL;
189 static int *pointi = NULL;
190 static int *xptr = NULL, *modptr = NULL, *evtspt = NULL;
191 static int *funtyp = NULL, *inpptr = NULL, *outptr = NULL, *inplnk = NULL, *outlnk = NULL;
192 static int *clkptr = NULL, *ordptr = NULL, *ordclk = NULL, *cord = NULL, *iord = NULL, *oord = NULL,  *zord = NULL,  *critev = NULL,  *zcptr = NULL;
193 static int nblk = 0, nordptr = 0, nlnk = 0, nx = 0, ng = 0, ncord = 0, noord = 0, nzord = 0;
194 static int nordclk = 0, niord = 0, nmod = 0;
195 static int nelem = 0;
196 static int *mod = NULL;
197 static int *neq = NULL;
198 static int *xprop = NULL; /* xproperty */
199 static int debug_block = 0;
200 static int *iwa = NULL;
201 static int hot;
202
203 /* declaration of ptr for typed port */
204 static void **outtbptr = NULL; /*pointer array of object of outtb*/
205 static int *outtbsz = NULL;    /*size of object of outtb*/
206 static int *outtbtyp = NULL;   /*type of object of outtb*/
207
208 SCSREAL_COP *outtbdptr = NULL;     /*to store double of outtb*/
209 SCSINT8_COP *outtbcptr = NULL;     /*to store int8 of outtb*/
210 SCSINT16_COP *outtbsptr = NULL;    /*to store int16 of outtb*/
211 SCSINT32_COP *outtblptr = NULL;    /*to store int32 of outtb*/
212 SCSUINT8_COP *outtbucptr = NULL;   /*to store unsigned int8 of outtb */
213 SCSUINT16_COP *outtbusptr = NULL;  /*to store unsigned int16 of outtb */
214 SCSUINT32_COP *outtbulptr = NULL;  /*to store unsigned int32 of outtb */
215
216 static outtb_el *outtb_elem = NULL;
217
218 static scicos_block *Blocks = NULL;
219
220 /* pass to external variable for code generation */
221 /* reserved variable name */
222 int *block_error = NULL;
223 double scicos_time = 0.;
224 int phase = 0;
225 int Jacobian_Flag = 0;
226 // double CI = 0., CJ = 0.;  // doubles returned by Get_Jacobian_ci and Get_Jacobian_cj respectively
227 double  CJJ = 0.;            // returned by Get_Jacobian_parameter
228 double SQuround = 0.;
229 /* Jacobian*/
230 static int AJacobian_block = 0;
231
232
233 /* Variable declaration moved to scicos.c because it was in the scicos-def.h therefore
234 * multiple declaration of the variable and linkers were complaining about duplicate
235 * symbols
236 */
237 SCICOS_IMPEXP COSDEBUGCOUNTER_struct C2F(cosdebugcounter);
238 SCICOS_IMPEXP RTFACTOR_struct C2F(rtfactor);
239 SCICOS_IMPEXP SOLVER_struct C2F(cmsolver);
240 SCICOS_IMPEXP CURBLK_struct C2F(curblk);
241 SCICOS_IMPEXP COSDEBUG_struct C2F(cosdebug);
242 SCICOS_IMPEXP COSHLT_struct C2F(coshlt);
243 SCICOS_IMPEXP DBCOS_struct C2F(dbcos);
244 SCICOS_IMPEXP COSTOL_struct C2F(costol);
245 SCICOS_IMPEXP COSERR_struct coserr;
246 /*--------------------------------------------------------------------------*/
247 static void FREE_blocks();
248 static void cosini(double *told);
249 static void cosend(double *told);
250 static void cossim(double *told);
251 static void cossimdaskr(double *told);
252 static void doit(double *told);
253 static void idoit(double *told);
254 static void zdoit(double *told, double *xt, double *xtd, double *g);
255 static void cdoit(double *told);
256 static void ozdoit(double *told, double *xt, double *xtd, int *kiwa);
257 static void odoit(double *told, double *xt, double *xtd, double *residual);
258 static void ddoit(double *told);
259 static void edoit(double *told, int *kiwa);
260 static void reinitdoit(double *told);
261 static int CallKinsol(double *told);
262 static int simblk(realtype t, N_Vector yy, N_Vector yp, void *f_data);
263 static int grblkdaskr(realtype t, N_Vector yy, N_Vector yp, realtype *gout, void *g_data);
264 static int grblk(realtype t, N_Vector yy, realtype *gout, void *g_data);
265 static void simblklsodar(int * nequations, realtype * tOld, realtype * actual, realtype * res);
266 static void grblklsodar(int * nequations, realtype * tOld, realtype * actual, int * ngc, realtype * res);
267 static void simblkddaskr(realtype *tOld, realtype *actual, realtype *actualP, realtype *res, int *flag, double *dummy1, int *dummy2);
268 static void grblkddaskr(int *nequations, realtype *tOld, realtype *actual, int *ngc, realtype *res, double *dummy1, int *dummy2);
269 static void jacpsol(realtype *res, int *ires, int *nequations, realtype *tOld, realtype *actual, realtype *actualP,
270                     realtype *rewt, realtype *savr, realtype *wk, realtype *h, realtype *cj, realtype *wp,
271                     int *iwp, int *ier, double *dummy1, int *dummy2);
272 static void psol(int *nequations, realtype *tOld, realtype *actual, realtype *actualP,
273                  realtype *savr, realtype *wk, realtype *cj, realtype *wght, realtype *wp,
274                  int *iwp, realtype *b, realtype *eplin, int *ier, double *dummy1, int *dummy2);
275 static void addevs(double t, int *evtnb, int *ierr1);
276 static int synchro_g_nev(ScicosImport *scs_imp, double *g, int kf, int *ierr);
277 static void Multp(double *A, double *B, double *R, int ra, int rb, int ca, int cb);
278 static int read_id(ezxml_t *elements, char *id, double *value);
279 static int simblkdaskr(realtype tres, N_Vector yy, N_Vector yp, N_Vector resval, void *rdata);
280 static void SundialsErrHandler(int error_code, const char *module, const char *function, char *msg, void *user_data);
281 static int Jacobians(long int Neq, realtype tt, realtype cj, N_Vector yy,
282                      N_Vector yp, N_Vector resvec, DlsMat Jacque, void *jdata,
283                      N_Vector tempv1, N_Vector tempv2, N_Vector tempv3);
284 static void call_debug_scicos(scicos_block *block, scicos_flag *flag, int flagi, int deb_blk);
285 static int synchro_nev(ScicosImport *scs_imp, int kf, int *ierr);
286 /*--------------------------------------------------------------------------*/
287 extern int C2F(dset)(int *n, double *dx, double *dy, int *incy);
288 extern int C2F(dcopy)(int *, double *, int *, double *, int *);
289 extern int C2F(dgefa)(double *A, int *lead_dim_A, int *n, int *ipivots, int *info);
290 extern int C2F(dgesl)(double *A, int *lead_dim_A, int *n, int *ipivots, double *B, int *job);
291 extern int C2F(msgs)();
292 extern void C2F(clearscicosimport)();
293 /*--------------------------------------------------------------------------*/
294 void putevs(double *t, int *evtnb, int *ierr1);
295 void Jdoit(double *told, double *xt, double *xtd, double *residual, int *job);
296 int simblkKinsol(N_Vector yy, N_Vector resval, void *rdata);
297 /*--------------------------------------------------------------------------*/
298 int C2F(scicos)(double *x_in, int *xptr_in, double *z__,
299                 void **work, int *zptr, int *modptr_in,
300                 void **oz, int *ozsz, int *oztyp, int *ozptr,
301                 char **iz, int *izptr, char **uid, int *uidptr, double *t0_in,
302                 double *tf_in, double *tevts_in, int *evtspt_in,
303                 int *nevts, int *pointi_in, void **outtbptr_in,
304                 int *outtbsz_in, int *outtbtyp_in,
305                 outtb_el *outtb_elem_in, int *nelem1, int *nlnk1,
306                 void** funptr, int *funtyp_in, int *inpptr_in,
307                 int *outptr_in, int *inplnk_in, int *outlnk_in,
308                 double *rpar, int *rpptr, int *ipar, int *ipptr,
309                 void **opar, int *oparsz, int *opartyp, int *opptr,
310                 int *clkptr_in, int *ordptr_in, int *nordptr1,
311                 int *ordclk_in, int *cord_in, int *ncord1,
312                 int *iord_in, int *niord1, int *oord_in,
313                 int *noord1, int *zord_in, int *nzord1,
314                 int *critev_in, int *nblk1, int *ztyp,
315                 int *zcptr_in, int *subscr, int *nsubs,
316                 double *simpar, int *flag__, int *ierr_out)
317 {
318     int i1, kf, lprt, in, out, job = 1;
319
320
321     static int mxtb = 0, ierr0 = 0, kfun0 = 0, i = 0, j = 0, k = 0, jj = 0;
322     static int ni = 0, no = 0;
323     static int nz = 0, noz = 0, nopar = 0;
324     double *W = NULL;
325
326     // Set FPU Flag to Extended for scicos simulation
327     // in order to override Java setting it to Double.
328 #if defined(linux) && defined(__i386__)
329     setFPUToExtended();
330 #endif
331     /*     Copyright INRIA */
332     /* iz,izptr are used to pass block labels */
333     TCritWarning = 0;
334
335     t0 = t0_in;
336     tf = tf_in;
337     ierr = ierr_out;
338
339     /* Parameter adjustments */
340     pointi = pointi_in;
341     x = x_in;
342     xptr = xptr_in - 1;
343     modptr = modptr_in - 1;
344     --zptr;
345     //--izptr;
346     --ozptr;
347     evtspt = evtspt_in - 1;
348     tevts = tevts_in - 1;
349     outtbptr = outtbptr_in;
350     outtbsz = outtbsz_in;
351     outtbtyp = outtbtyp_in;
352     outtb_elem = outtb_elem_in;
353     funtyp = funtyp_in - 1;
354     inpptr = inpptr_in - 1;
355     outptr = outptr_in - 1;
356     inplnk = inplnk_in - 1;
357     outlnk = outlnk_in - 1;
358     --rpptr;
359     --ipptr;
360     --opptr;
361     clkptr = clkptr_in - 1;
362     ordptr = ordptr_in - 1;
363     ordclk = ordclk_in - 1;
364     cord = cord_in - 1;
365     iord = iord_in - 1;
366     oord = oord_in - 1;
367     zord = zord_in - 1;
368
369     critev = critev_in - 1;
370     --ztyp;
371     zcptr = zcptr_in - 1;
372     --simpar;
373
374     /* Function Body */
375     Atol = simpar[1];
376     rtol = simpar[2];
377     ttol = simpar[3];
378     deltat = simpar[4];
379     C2F(rtfactor).scale = simpar[5];
380     C2F(cmsolver).solver = (int) simpar[6];
381     hmax = simpar[7];
382
383     nordptr = *nordptr1;
384     nblk  = *nblk1;
385     ncord = *ncord1;
386     noord = *noord1;
387     nzord = *nzord1;
388     niord = *niord1;
389     nlnk  = *nlnk1;
390     nelem = *nelem1;
391     *ierr = 0;
392
393     nordclk = ordptr[nordptr] - 1;  /* number of rows in ordclk is ordptr(nclkp1)-1 */
394     ng      = zcptr[nblk + 1] - 1;  /* computes number of zero crossing surfaces */
395     nmod    = modptr[nblk + 1] - 1; /* computes number of modes */
396     nz      = zptr[nblk + 1] - 1;   /* number of discrete real states */
397     noz     = ozptr[nblk + 1] - 1;  /* number of discrete object states */
398     nopar   = opptr[nblk + 1] - 1;  /* number of object parameters */
399     nx      = xptr[nblk + 1] - 1;   /* number of object parameters */
400     neq     = &nx;
401
402     xd = &x[xptr[nblk + 1] - 1];
403
404     /* check for hard coded maxsize */
405     for (i = 1; i <= nblk; ++i)
406     {
407         if (funtyp[i] < 10000)
408         {
409             funtyp[i] %= 1000;
410         }
411         else
412         {
413             funtyp[i] = funtyp[i] % 1000 + 10000;
414         }
415         ni = inpptr[i + 1] - inpptr[i];
416         no = outptr[i + 1] - outptr[i];
417         if (funtyp[i] == 1)
418         {
419             if (ni + no > 11)
420             {
421                 /*     hard coded maxsize in callf.c */
422                 C2F(msgs)(&c__90, &c__0);
423                 C2F(curblk).kfun = i;
424                 *ierr = i + 1005;
425                 return 0;
426             }
427         }
428         else if (funtyp[i] == 2 || funtyp[i] == 3)
429         {
430             /*     hard coded maxsize in scicos.h */
431             if (ni + no > SZ_SIZE)
432             {
433                 C2F(msgs)(&c__90, &c__0);
434                 C2F(curblk).kfun = i;
435                 *ierr = i + 1005;
436                 return 0;
437             }
438         }
439         mxtb = 0;
440         if (funtyp[i] == 0)
441         {
442             if (ni > 1)
443             {
444                 for (j = 1; j <= ni; ++j)
445                 {
446                     k = inplnk[inpptr[i] - 1 + j];
447                     mxtb = mxtb + (outtbsz[k - 1] * outtbsz[(k - 1) + nlnk]);
448                 }
449             }
450             if (no > 1)
451             {
452                 for (j = 1; j <= no; ++j)
453                 {
454                     k = outlnk[outptr[i] - 1 + j];
455                     mxtb = mxtb + (outtbsz[k - 1] * outtbsz[(k - 1) + nlnk]);
456                 }
457             }
458             if (mxtb > TB_SIZE)
459             {
460                 C2F(msgs)(&c__91, &c__0);
461                 C2F(curblk).kfun = i;
462                 *ierr = i + 1005;
463                 return 0;
464             }
465         }
466     }
467
468     if (nx > 0) /* xprop */
469     {
470         if ((xprop = MALLOC(sizeof(int) * nx)) == NULL )
471         {
472             *ierr = 5;
473             return 0;
474         }
475     }
476     for (i = 0; i < nx; i++) /* initialize */
477     {
478         xprop[i] = 1;
479     }
480     if (nmod > 0) /* mod */
481     {
482         if ((mod = MALLOC(sizeof(int) * nmod)) == NULL )
483         {
484             *ierr = 5;
485             if (nx > 0)
486             {
487                 FREE(xprop);
488             }
489             return 0;
490         }
491     }
492     if (ng > 0) /* g becomes global */
493     {
494         if ((g = MALLOC(sizeof(double) * ng)) == NULL )
495         {
496             *ierr = 5;
497             if (nmod > 0)
498             {
499                 FREE(mod);
500             }
501             if (nx > 0)
502             {
503                 FREE(xprop);
504             }
505             return 0;
506         }
507     }
508
509     debug_block = -1; /* no debug block for start */
510     C2F(cosdebugcounter).counter = 0;
511
512     /** Create Block's array **/
513     if ((Blocks = MALLOC(sizeof(scicos_block) * nblk)) == NULL )
514     {
515         *ierr = 5;
516         if (nx > 0)
517         {
518             FREE(xprop);
519         }
520         if (nmod > 0)
521         {
522             FREE(mod);
523         }
524         if (ng > 0)
525         {
526             FREE(g);
527         }
528         return 0;
529     }
530
531     /** Setting blocks properties for each entry in Block's array **/
532
533     /* 1 : type and pointer on simulation function */
534     for (kf = 0; kf < nblk; ++kf)   /*for each block */
535     {
536         void* p = funptr[kf];
537         C2F(curblk).kfun = kf + 1;
538         Blocks[kf].type = funtyp[kf + 1];
539         if (Blocks[kf].type == IFTHEL_BLOCK)
540         {
541             funtyp[kf + 1] = -1;
542         }
543         else if (Blocks[kf].type == ESELECT_BLOCK)
544         {
545             funtyp[kf + 1] = -2;
546         }
547         else if (Blocks[kf].type < 0)
548         {
549             //macros
550             funtyp[kf + 1] *= -1; // Restore a positive 'funtyp' for later use
551             switch (-Blocks[kf].type)
552             {
553                 case FORTRAN_GATEWAY:
554                     Blocks[kf].funpt = (voidg)F2C(sciblk);
555                     break;
556                 case UNUSED1:
557                     sciprint(_("type 1 function not allowed for scilab blocks\n"));
558                     *ierr = 1000 + kf + 1;
559                     FREE_blocks();
560                     return 0;
561                 case UNUSED2:
562                     sciprint(_("type 2 function not allowed for scilab blocks\n"));
563                     *ierr = 1000 + kf + 1;
564                     FREE_blocks();
565                     return 0;
566                 case SCIFUNC_BLOCK:
567                     Blocks[kf].funpt = (voidg)sciblk2;
568                     Blocks[kf].type = 2;
569                     break;
570                 case EXPRESSION_BLOCK:
571                     Blocks[kf].funpt = (voidg)sciblk4;
572                     Blocks[kf].type = 4;
573                     break;
574                 case DEBUG_BLOCK: /* debugging block */
575                     Blocks[kf].funpt = (voidg)sciblk4;
576                     Blocks[kf].type *= -1;
577                     debug_block = kf;
578                     break;
579
580                 case OLD_SCI_BLOCK:
581                     Blocks[kf].funpt = (voidg)sciblk4;
582                     Blocks[kf].type = 10004;
583                     break;
584                 default:
585                     sciprint(_("Undefined Function type\n"));
586                     *ierr = 1000 + kf + 1;
587                     FREE_blocks();
588                     return 0;
589             }
590             Blocks[kf].scsptr = p; /* set scilab function adress for sciblk */
591         }
592         else
593         {
594             //linked functions (internal or external)
595             Blocks[kf].funpt = (voidf)p;
596             Blocks[kf].scsptr = NULL;   /* this is done for being able to test if a block
597                                         is a scilab block in the debugging phase when
598                                         sciblk4 is called */
599         }
600
601         /* 2 : Dimension properties */
602         Blocks[kf].ztyp = ztyp[kf + 1];
603         Blocks[kf].nx = xptr[kf + 2] - xptr[kf + 1]; /* continuous state dimension*/
604         Blocks[kf].ng = zcptr[kf + 2] - zcptr[kf + 1]; /* number of zero crossing surface*/
605         Blocks[kf].nz = zptr[kf + 2] - zptr[kf + 1]; /* number of double discrete state*/
606         Blocks[kf].noz = ozptr[kf + 2] - ozptr[kf + 1]; /* number of other discrete state*/
607         Blocks[kf].nrpar = rpptr[kf + 2] - rpptr[kf + 1]; /* size of double precision parameter vector*/
608         Blocks[kf].nipar = ipptr[kf + 2] - ipptr[kf + 1]; /* size of integer precision parameter vector*/
609         Blocks[kf].nopar = opptr[kf + 2] - opptr[kf + 1]; /* number of other parameters (matrix, data structure,..)*/
610         Blocks[kf].nin = inpptr[kf + 2] - inpptr[kf + 1]; /* number of input ports */
611         Blocks[kf].nout = outptr[kf + 2] - outptr[kf + 1]; /* number of output ports */
612
613         /* 3 : input port properties */
614         /* in insz, we store :
615         *  - insz[0..nin-1] : first dimension of input ports
616         *  - insz[nin..2*nin-1] : second dimension of input ports
617         *  - insz[2*nin..3*nin-1] : type of data of input ports
618         */
619         /* allocate size and pointer arrays (number of input ports)*/
620         Blocks[kf].insz = NULL;
621         Blocks[kf].inptr = NULL;
622         if (Blocks[kf].nin != 0)
623         {
624             if ((Blocks[kf].insz = MALLOC(Blocks[kf].nin * 3 * sizeof(int))) == NULL )
625             {
626                 FREE_blocks();
627                 *ierr = 5;
628                 return 0;
629             }
630             if ((Blocks[kf].inptr = MALLOC(Blocks[kf].nin * sizeof(double*))) == NULL )
631             {
632                 FREE_blocks();
633                 *ierr = 5;
634                 return 0;
635             }
636         }
637         for (in = 0; in < Blocks[kf].nin; in++)
638         {
639             lprt = inplnk[inpptr[kf + 1] + in];
640             Blocks[kf].inptr[in] = outtbptr[lprt - 1]; /* pointer on the data*/
641             Blocks[kf].insz[in] = outtbsz[lprt - 1]; /* row dimension of the input port*/
642             Blocks[kf].insz[Blocks[kf].nin + in] = outtbsz[(lprt - 1) + nlnk]; /* column dimension of the input port*/
643             Blocks[kf].insz[2 * Blocks[kf].nin + in] = outtbtyp[lprt - 1]; /*type of data of the input port*/
644         }
645         /* 4 : output port properties */
646         /* in outsz, we store :
647         *  - outsz[0..nout-1] : first dimension of output ports
648         *  - outsz[nout..2*nout-1] : second dimension of output ports
649         *  - outsz[2*nout..3*nout-1] : type of data of output ports
650         */
651         /* allocate size and pointer arrays (number of output ports)*/
652         Blocks[kf].outsz = NULL;
653         Blocks[kf].outptr = NULL;
654         if (Blocks[kf].nout != 0)
655         {
656             if ((Blocks[kf].outsz = MALLOC(Blocks[kf].nout * 3 * sizeof(int))) == NULL )
657             {
658                 FREE_blocks();
659                 *ierr = 5;
660                 return 0;
661             }
662             if ((Blocks[kf].outptr = MALLOC(Blocks[kf].nout * sizeof(double*))) == NULL )
663             {
664                 FREE_blocks();
665                 *ierr = 5;
666                 return 0;
667             }
668         }
669         /* set the values */
670         for (out = 0; out < Blocks[kf].nout; out++) /*for each output port */
671         {
672             lprt = outlnk[outptr[kf + 1] + out];
673             Blocks[kf].outptr[out] = outtbptr[lprt - 1]; /*pointer on data */
674             Blocks[kf].outsz[out] = outtbsz[lprt - 1]; /*row dimension of output port*/
675             Blocks[kf].outsz[Blocks[kf].nout + out] = outtbsz[(lprt - 1) + nlnk]; /*column dimension of output ports*/
676             Blocks[kf].outsz[2 * Blocks[kf].nout + out] = outtbtyp[lprt - 1]; /*type of data of output port */
677         }
678
679         /* 5 : event output port properties */
680         Blocks[kf].evout = NULL;
681         Blocks[kf].nevout = clkptr[kf + 2] - clkptr[kf + 1];
682         if (Blocks[kf].nevout != 0)
683         {
684             if ((Blocks[kf].evout = CALLOC(Blocks[kf].nevout, sizeof(double))) == NULL )
685             {
686                 FREE_blocks();
687                 *ierr = 5;
688                 return 0;
689             }
690         }
691
692         /* 6 : pointer on the begining of the double discrete state array ( z) */
693         Blocks[kf].z = &(z__[zptr[kf + 1] - 1]);
694
695         /* 7 : type, size and pointer on the other discrete states  data structures (oz) */
696         Blocks[kf].ozsz = NULL;
697         if (Blocks[kf].noz == 0)
698         {
699             Blocks[kf].ozptr = NULL;
700             Blocks[kf].oztyp = NULL;
701         }
702         else
703         {
704             Blocks[kf].ozptr = &(oz[ozptr[kf + 1] - 1]);
705             if ((Blocks[kf].ozsz = MALLOC(Blocks[kf].noz * 2 * sizeof(int))) == NULL )
706             {
707                 FREE_blocks();
708                 *ierr = 5;
709                 return 0;
710             }
711             for (i = 0; i < Blocks[kf].noz; i++)
712             {
713                 Blocks[kf].ozsz[i] = ozsz[(ozptr[kf + 1] - 1) + i];
714                 Blocks[kf].ozsz[i + Blocks[kf].noz] = ozsz[(ozptr[kf + 1] - 1 + noz) + i];
715             }
716             Blocks[kf].oztyp = &(oztyp[ozptr[kf + 1] - 1]);
717         }
718
719         /* 8 : pointer on the begining of the double parameter array ( rpar ) */
720         Blocks[kf].rpar = &(rpar[rpptr[kf + 1] - 1]);
721
722         /* 9 : pointer on the begining of the integer parameter array ( ipar ) */
723         Blocks[kf].ipar = &(ipar[ipptr[kf + 1] - 1]);
724
725         /* 10 : type, size and pointer on the other parameters  data structures (opar) */
726         Blocks[kf].oparsz = NULL;
727         if (Blocks[kf].nopar == 0)
728         {
729             Blocks[kf].oparptr = NULL;
730             Blocks[kf].opartyp = NULL;
731         }
732         else
733         {
734             Blocks[kf].oparptr = &(opar[opptr[kf + 1] - 1]);
735             if ((Blocks[kf].oparsz = MALLOC(Blocks[kf].nopar * 2 * sizeof(int))) == NULL)
736             {
737                 FREE_blocks();
738                 *ierr = 5;
739                 return 0;
740             }
741             for (i = 0; i < Blocks[kf].nopar; i++)
742             {
743                 Blocks[kf].oparsz[i] = oparsz[(opptr[kf + 1] - 1) + i];
744                 Blocks[kf].oparsz[i + Blocks[kf].nopar] = oparsz[(opptr[kf + 1] - 1 + nopar) + i];
745             }
746             Blocks[kf].opartyp = &(opartyp[opptr[kf + 1] - 1]);
747         }
748
749         /* 10 : pointer on the beginning of the residual array (res) */
750         Blocks[kf].res = NULL;
751         if (Blocks[kf].nx != 0)
752         {
753             if ((Blocks[kf].res = MALLOC(Blocks[kf].nx * sizeof(double))) == NULL)
754             {
755                 FREE_blocks();
756                 *ierr = 5;
757                 return 0;
758             }
759         }
760
761         /* 11 : block label (label) */
762         i1 = izptr[kf];
763         if ((Blocks[kf].label = MALLOC(sizeof(char) * (i1 + 1))) == NULL)
764         {
765             FREE_blocks();
766             *ierr = 5;
767             return 0;
768         }
769         Blocks[kf].label[i1] = '\0';
770         strcpy(Blocks[kf].label, iz[kf]);
771
772         /* block uid (uid) */
773         i1 = uidptr[kf];
774         if ((Blocks[kf].uid = MALLOC(sizeof(char) * (i1 + 1))) == NULL)
775         {
776             FREE_blocks();
777             *ierr = 5;
778             return 0;
779         }
780         Blocks[kf].uid[i1] = '\0';
781         strcpy(Blocks[kf].uid, uid[kf]);
782
783         /* 12 : block array of crossed surfaces (jroot) */
784         Blocks[kf].jroot = NULL;
785         if (Blocks[kf].ng > 0)
786         {
787             if ((Blocks[kf].jroot = CALLOC(Blocks[kf].ng, sizeof(int))) == NULL)
788             {
789                 FREE_blocks();
790                 *ierr = 5;
791                 return 0;
792             }
793         }
794
795         /* 13 : block work array  (work) */
796         Blocks[kf].work = (void **)(((double *)work) + kf);
797
798         /* 14 : block modes  array  (mode) */
799         Blocks[kf].nmode = modptr[kf + 2] - modptr[kf + 1];
800         if (Blocks[kf].nmode != 0)
801         {
802             Blocks[kf].mode = &(mod[modptr[kf + 1] - 1]);
803         }
804
805         /* 15 : block xprop  array  (xprop) */
806         Blocks[kf].xprop = NULL;
807         if (Blocks[kf].nx != 0)
808         {
809             Blocks[kf].xprop = &(xprop[xptr[kf + 1] - 1]);
810         }
811
812         /* 16 : pointer on the zero crossing surface computation function of the block (g) */
813         Blocks[kf].g = NULL;
814         if (Blocks[kf].ng != 0)
815         {
816             Blocks[kf].g = &(g[zcptr[kf + 1] - 1]);
817         }
818     }
819     /** all block properties are stored in the Blocks array **/
820
821     /* iwa */
822     iwa = NULL;
823     if ((*nevts) != 0)
824     {
825         if ((iwa = MALLOC(sizeof(int) * (*nevts))) == NULL)
826         {
827             FREE_blocks();
828             *ierr = 5;
829             return 0;
830         }
831     }
832
833     /* save ptr of scicos in import structure */
834     makescicosimport(x, &nx, &xptr[1], &zcptr[1], z__, &nz, &zptr[1],
835                      &noz, oz, ozsz, oztyp, &ozptr[1],
836                      g, &ng, mod, &nmod, &modptr[1], iz, &izptr[1],
837                      uid, uidptr,
838                      &inpptr[1], &inplnk[1], &outptr[1], &outlnk[1],
839                      outtbptr, outtbsz, outtbtyp,
840                      outtb_elem, &nelem,
841                      &nlnk, rpar, &rpptr[1], ipar, &ipptr[1],
842                      opar, oparsz, opartyp, &opptr[1],
843                      &nblk, subscr, nsubs,
844                      &tevts[1], &evtspt[1], nevts, pointi,
845                      &iord[1], &niord, &oord[1], &noord, &zord[1], &nzord,
846                      funptr, &funtyp[1], &ztyp[1],
847                      &cord[1], &ncord, &ordclk[1], &nordclk, &clkptr[1],
848                      &ordptr[1], &nordptr, &critev[1], iwa, Blocks,
849                      t0, tf, &Atol, &rtol, &ttol, &deltat, &hmax,
850                      xprop, xd);
851
852     if (*flag__ == 1)   /*start*/
853     {
854         /*      blocks initialization */
855         for (kf = 0; kf < nblk; ++kf)
856         {
857             *(Blocks[kf].work) = NULL;
858         }
859         cosini(t0);
860         if (*ierr != 0)
861         {
862             ierr0 = *ierr;
863             kfun0 = C2F(curblk).kfun;
864             cosend(t0);
865             *ierr = ierr0;
866             C2F(curblk).kfun = kfun0;
867         }
868
869     }
870     else if (*flag__ == 2)     /*run*/
871     {
872
873         /*     integration */
874         switch (C2F(cmsolver).solver)
875         {
876             case LSodar_Dynamic:
877             case CVode_BDF_Newton:
878             case CVode_BDF_Functional:
879             case CVode_Adams_Newton:
880             case CVode_Adams_Functional:
881             case Dormand_Prince:
882             case Runge_Kutta:
883             case Implicit_Runge_Kutta:
884                 cossim(t0);
885                 break;
886             case IDA_BDF_Newton:
887             case DDaskr_BDF_Newton:
888             case DDaskr_BDF_GMRes:
889                 cossimdaskr(t0);
890                 break;
891             default: // Unknown solver number
892                 *ierr = 1000;
893                 return 0;
894         }
895         if (*ierr != 0)
896         {
897             ierr0 = *ierr;
898             kfun0 = C2F(curblk).kfun;
899             cosend(t0);
900             *ierr = ierr0;
901             C2F(curblk).kfun = kfun0;
902         }
903
904     }
905     else if (*flag__ == 3)     /*finish*/
906     {
907         /*     blocks closing */
908         cosend(t0);
909     }
910     else if (*flag__ == 4)     /*linear*/
911     {
912         phase = 1;
913         idoit(t0);
914         if (*ierr == 0)
915         {
916             if ((W = MALLOC(sizeof(double) * (Max(nx, ng)))) == NULL )
917             {
918                 FREE(iwa);
919                 FREE_blocks();
920                 *ierr = 5;
921                 return 0;
922             }
923
924             /*---------instead of old simblk--------*/
925             /*  C2F(simblk)(&nx, t0, x, W);  */
926
927             if (ng > 0 && nmod > 0)
928             {
929                 zdoit(t0, x, x + nx, W); /* updating modes as a function of state values; this was necessary in iGUI*/
930             }
931             for (jj = 0; jj < nx; jj++)
932             {
933                 W[jj] = 0.0;
934             }
935             C2F(ierode).iero = 0;
936             *ierr = 0;
937             if (C2F(cmsolver).solver < 100)
938             {
939                 odoit(t0, x, W, W);
940             }
941             else
942             {
943                 odoit(t0, x, x + nx, W);
944             }
945             C2F(ierode).iero = *ierr;
946             /*-----------------------------------------*/
947             for (i = 0; i < nx; ++i)
948             {
949                 x[i] = W[i];
950             }
951             FREE(W);
952         }
953     }
954     else if (*flag__ == 5)     /* initial_KINSOL= "Kinsol" */
955     {
956         C2F(ierode).iero = 0;
957         *ierr = 0;
958         idoit(t0);
959         CallKinsol(t0);
960         *ierr = C2F(ierode).iero;
961     }
962
963
964     FREE(iwa);
965     FREE_blocks();
966
967     C2F(clearscicosimport)();
968     return 0;
969 } /* scicos_ */
970 /*--------------------------------------------------------------------------*/
971 /* check_flag */
972 static int check_flag(void *flagvalue, char *funcname, int opt)
973 {
974     int *errflag = NULL;
975
976     /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
977     if (opt == 0 && flagvalue == NULL)
978     {
979         sciprint(_("\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n"), funcname);
980         return (1);
981     }
982     /* Check if flag < 0 */
983     else if (opt == 1)
984     {
985         errflag = (int *) flagvalue;
986         if (*errflag < 0)
987         {
988             sciprint(_("\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n"),
989                      funcname, *errflag);
990             return (1);
991         }
992     }
993     /* Check if function returned NULL pointer - no memory allocated */
994     else if (opt == 2 && flagvalue == NULL)
995     {
996         sciprint(_("\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n"), funcname);
997         return (1);
998     }
999
1000     return (0);
1001 } /* check_flag */
1002
1003 /*--------------------------------------------------------------------------*/
1004 static void cosini(double *told)
1005 {
1006     static scicos_flag flag__ = 0;
1007     static int i = 0;
1008
1009     static int kfune = 0;
1010     static int jj = 0;
1011
1012     SCSREAL_COP *outtbd = NULL;    /*to save double of outtb*/
1013     SCSINT8_COP *outtbc = NULL;    /*to save int8 of outtb*/
1014     SCSINT16_COP *outtbs = NULL;   /*to save int16 of outtb*/
1015     SCSINT32_COP *outtbl = NULL;   /*to save int32 of outtb*/
1016     SCSUINT8_COP *outtbuc = NULL;  /*to save unsigned int8 of outtb*/
1017     SCSUINT16_COP *outtbus = NULL; /*to save unsigned int16 of outtb*/
1018     SCSUINT32_COP *outtbul = NULL; /*to save unsigned int32 of outtb*/
1019     int szouttbd = 0;  /*size of arrays*/
1020     int szouttbc = 0,  szouttbs = 0,  szouttbl = 0;
1021     int szouttbuc = 0, szouttbus = 0, szouttbul = 0;
1022     int curouttbd = 0; /*current position in arrays*/
1023     int curouttbc = 0,  curouttbs = 0,  curouttbl = 0;
1024     int curouttbuc = 0, curouttbus = 0, curouttbul = 0;
1025
1026     int ii = 0, kk = 0; /*local counters*/
1027     int sszz = 0;  /*local size of element of outtb*/
1028     /*Allocation of arrays for outtb*/
1029     for (ii = 0; ii < nlnk; ii++)
1030     {
1031         switch (outtbtyp[ii])
1032         {
1033             case SCSREAL_N    :
1034                 szouttbd += outtbsz[ii] * outtbsz[ii + nlnk]; /*double real matrix*/
1035                 outtbd = (SCSREAL_COP *) REALLOC (outtbd, szouttbd * sizeof(SCSREAL_COP));
1036                 memset(outtbd, 0, szouttbd * sizeof(SCSREAL_COP));
1037                 break;
1038
1039             case SCSCOMPLEX_N :
1040                 szouttbd += 2 * outtbsz[ii] * outtbsz[ii + nlnk]; /*double complex matrix*/
1041                 outtbd = (SCSCOMPLEX_COP *) REALLOC (outtbd, szouttbd * sizeof(SCSCOMPLEX_COP));
1042                 memset(outtbd, 0, szouttbd * sizeof(SCSCOMPLEX_COP));
1043                 break;
1044
1045             case SCSINT8_N    :
1046                 szouttbc += outtbsz[ii] * outtbsz[ii + nlnk]; /*int8*/
1047                 outtbc = (SCSINT8_COP *) REALLOC (outtbc, szouttbc * sizeof(SCSINT8_COP));
1048                 memset(outtbc, 0, szouttbc * sizeof(SCSINT8_COP));
1049                 break;
1050
1051             case SCSINT16_N   :
1052                 szouttbs += outtbsz[ii] * outtbsz[ii + nlnk]; /*int16*/
1053                 outtbs = (SCSINT16_COP *) REALLOC (outtbs, szouttbs * sizeof(SCSINT16_COP));
1054                 memset(outtbs, 0, szouttbs * sizeof(SCSINT16_COP));
1055                 break;
1056
1057             case SCSINT32_N   :
1058                 szouttbl += outtbsz[ii] * outtbsz[ii + nlnk]; /*int32*/
1059                 outtbl = (SCSINT32_COP *) REALLOC (outtbl, szouttbl * sizeof(SCSINT32_COP));
1060                 memset(outtbl, 0, szouttbl * sizeof(SCSINT32_COP));
1061                 break;
1062
1063             case SCSUINT8_N   :
1064                 szouttbuc += outtbsz[ii] * outtbsz[ii + nlnk]; /*uint8*/
1065                 outtbuc = (SCSUINT8_COP *) REALLOC (outtbuc, szouttbuc * sizeof(SCSUINT8_COP));
1066                 memset(outtbuc, 0, szouttbuc * sizeof(SCSUINT8_COP));
1067                 break;
1068
1069             case SCSUINT16_N  :
1070                 szouttbus += outtbsz[ii] * outtbsz[ii + nlnk]; /*uint16*/
1071                 outtbus = (SCSUINT16_COP *) REALLOC (outtbus, szouttbus * sizeof(SCSUINT16_COP));
1072                 memset(outtbus, 0, szouttbus * sizeof(SCSUINT16_COP));
1073                 break;
1074
1075             case SCSUINT32_N  :
1076                 szouttbul += outtbsz[ii] * outtbsz[ii + nlnk]; /*uint32*/
1077                 outtbul = (SCSUINT32_COP *) REALLOC (outtbul, szouttbul * sizeof(SCSUINT32_COP));
1078                 memset(outtbul, 0, szouttbul * sizeof(SCSUINT32_COP));
1079                 break;
1080
1081             default  : /* Add a message here */
1082                 break;
1083         }
1084     }
1085
1086     /* Jacobian*/
1087     AJacobian_block = 0;
1088
1089     /* Function Body */
1090     *ierr = 0;
1091
1092     /*     initialization (flag 4) */
1093     /*     loop on blocks */
1094     C2F(dset)(&ng, &c_b14, g, &c__1);
1095
1096     for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
1097     {
1098         flag__ = 4;
1099         if (Blocks[C2F(curblk).kfun - 1].nx > 0)
1100         {
1101             Blocks[C2F(curblk).kfun - 1].x  = &x[xptr[C2F(curblk).kfun] - 1];
1102             Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
1103         }
1104         Blocks[C2F(curblk).kfun - 1].nevprt = 0;
1105         if (funtyp[C2F(curblk).kfun] >= 0)   /* debug_block is not called here */
1106         {
1107             /*callf(told, xd, x, x,g,&flag__);*/
1108             Jacobian_Flag = 0;
1109             callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1110             if (flag__ < 0 && *ierr == 0)
1111             {
1112                 *ierr = 5 - flag__;
1113                 kfune = C2F(curblk).kfun;
1114             }
1115             if ((Jacobian_Flag == 1) && (AJacobian_block == 0))
1116             {
1117                 AJacobian_block = C2F(curblk).kfun;
1118             }
1119         }
1120     }
1121     if (*ierr != 0)
1122     {
1123         C2F(curblk).kfun = kfune;
1124         freeouttbptr;
1125         return;
1126     }
1127
1128     /*     initialization (flag 6) */
1129     flag__ = 6;
1130     for (jj = 1; jj <= ncord; ++jj)
1131     {
1132         C2F(curblk).kfun = cord[jj];
1133         Blocks[C2F(curblk).kfun - 1].nevprt = 0;
1134         if (funtyp[C2F(curblk).kfun] >= 0)
1135         {
1136             /*callf(told, xd, x, x,g,&flag__);*/
1137             callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1138             if (flag__ < 0)
1139             {
1140                 *ierr = 5 - flag__;
1141                 freeouttbptr;
1142                 return;
1143             }
1144         }
1145     }
1146
1147     /*     point-fix iterations */
1148     flag__ = 6;
1149     for (i = 1; i <= nblk + 1; ++i)   /*for each block*/
1150     {
1151         /*     loop on blocks */
1152         for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
1153         {
1154             Blocks[C2F(curblk).kfun - 1].nevprt = 0;
1155             if (funtyp[C2F(curblk).kfun] >= 0)
1156             {
1157                 /*callf(told, xd, x, x,g,&flag__);*/
1158                 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1159                 if (flag__ < 0)
1160                 {
1161                     *ierr = 5 - flag__;
1162                     freeouttbptr;
1163                     return;
1164                 }
1165             }
1166         }
1167
1168         flag__ = 6;
1169         for (jj = 1; jj <= ncord; ++jj)   /*for each continuous block*/
1170         {
1171             C2F(curblk).kfun = cord[jj];
1172             if (funtyp[C2F(curblk).kfun] >= 0)
1173             {
1174                 /*callf(told, xd, x, x,g,&flag__);*/
1175                 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1176                 if (flag__ < 0)
1177                 {
1178                     *ierr = 5 - flag__;
1179                     freeouttbptr;
1180                     return;
1181                 }
1182             }
1183         }
1184
1185         /*comparison between outtb and arrays*/
1186         curouttbd = 0;
1187         curouttbc = 0;
1188         curouttbs = 0;
1189         curouttbl = 0;
1190         curouttbuc = 0;
1191         curouttbus = 0;
1192         curouttbul = 0;
1193         for (jj = 0; jj < nlnk; jj++)
1194         {
1195             switch (outtbtyp[jj]) /*for each type of ports*/
1196             {
1197                 case SCSREAL_N    :
1198                     outtbdptr = (SCSREAL_COP *)outtbptr[jj]; /*double real matrix*/
1199                     sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1200                     for (kk = 0; kk < sszz; kk++)
1201                     {
1202                         int outtbdptr_isnan = outtbdptr[kk] != outtbdptr[kk];
1203                         int outtbd_isnan = (SCSREAL_COP)outtbd[curouttbd + kk] != (SCSREAL_COP)outtbd[curouttbd + kk];
1204
1205                         if (outtbdptr_isnan && outtbd_isnan)
1206                         {
1207                             continue;
1208                         }
1209                         if (outtbdptr[kk] != (SCSREAL_COP)outtbd[curouttbd + kk])
1210                         {
1211                             goto L30;
1212                         }
1213                     }
1214                     curouttbd += sszz;
1215                     break;
1216
1217                 case SCSCOMPLEX_N :
1218                     outtbdptr = (SCSCOMPLEX_COP *)outtbptr[jj]; /*double complex matrix*/
1219                     sszz = 2 * outtbsz[jj] * outtbsz[jj + nlnk];
1220                     for (kk = 0; kk < sszz; kk++)
1221                     {
1222                         int outtbdptr_isnan = outtbdptr[kk] != outtbdptr[kk];
1223                         int outtbd_isnan = (SCSCOMPLEX_COP)outtbd[curouttbd + kk] != (SCSCOMPLEX_COP)outtbd[curouttbd + kk];
1224
1225                         if (outtbdptr_isnan && outtbd_isnan)
1226                         {
1227                             continue;
1228                         }
1229                         if (outtbdptr[kk] != (SCSCOMPLEX_COP)outtbd[curouttbd + kk])
1230                         {
1231                             goto L30;
1232                         }
1233                     }
1234                     curouttbd += sszz;
1235                     break;
1236
1237                 case SCSINT8_N    :
1238                     outtbcptr = (SCSINT8_COP *)outtbptr[jj]; /*int8*/
1239                     sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1240                     for (kk = 0; kk < sszz; kk++)
1241                     {
1242                         if (outtbcptr[kk] != (SCSINT8_COP)outtbc[curouttbc + kk])
1243                         {
1244                             goto L30;
1245                         }
1246                     }
1247                     curouttbc += sszz;
1248                     break;
1249
1250                 case SCSINT16_N   :
1251                     outtbsptr = (SCSINT16_COP *)outtbptr[jj]; /*int16*/
1252                     sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1253                     for (kk = 0; kk < sszz; kk++)
1254                     {
1255                         if (outtbsptr[kk] != (SCSINT16_COP)outtbs[curouttbs + kk])
1256                         {
1257                             goto L30;
1258                         }
1259                     }
1260                     curouttbs += sszz;
1261                     break;
1262
1263                 case SCSINT32_N   :
1264                     outtblptr = (SCSINT32_COP *)outtbptr[jj]; /*int32*/
1265                     sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1266                     for (kk = 0; kk < sszz; kk++)
1267                     {
1268                         if (outtblptr[kk] != (SCSINT32_COP)outtbl[curouttbl + kk])
1269                         {
1270                             goto L30;
1271                         }
1272                     }
1273                     curouttbl += sszz;
1274                     break;
1275
1276                 case SCSUINT8_N   :
1277                     outtbucptr = (SCSUINT8_COP *)outtbptr[jj]; /*uint8*/
1278                     sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1279                     for (kk = 0; kk < sszz; kk++)
1280                     {
1281                         if (outtbucptr[kk] != (SCSUINT8_COP)outtbuc[curouttbuc + kk])
1282                         {
1283                             goto L30;
1284                         }
1285                     }
1286                     curouttbuc += sszz;
1287                     break;
1288
1289                 case SCSUINT16_N  :
1290                     outtbusptr = (SCSUINT16_COP *)outtbptr[jj]; /*uint16*/
1291                     sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1292                     for (kk = 0; kk < sszz; kk++)
1293                     {
1294                         if (outtbusptr[kk] != (SCSUINT16_COP)outtbus[curouttbus + kk])
1295                         {
1296                             goto L30;
1297                         }
1298                     }
1299                     curouttbus += sszz;
1300                     break;
1301
1302                 case SCSUINT32_N  :
1303                     outtbulptr = (SCSUINT32_COP *)outtbptr[jj]; /*uint32*/
1304                     sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1305                     for (kk = 0; kk < sszz; kk++)
1306                     {
1307                         if (outtbulptr[kk] != (SCSUINT32_COP)outtbul[curouttbul + kk])
1308                         {
1309                             goto L30;
1310                         }
1311                     }
1312                     curouttbul += sszz;
1313                     break;
1314
1315                 default  : /* Add a message here */
1316                     break;
1317             }
1318         }
1319         freeouttbptr;
1320         return;
1321
1322 L30:
1323         /*Save data of outtb in arrays*/
1324         curouttbd = 0;
1325         curouttbc = 0;
1326         curouttbs = 0;
1327         curouttbl = 0;
1328         curouttbuc = 0;
1329         curouttbus = 0;
1330         curouttbul = 0;
1331         for (ii = 0; ii < nlnk; ii++) /*for each link*/
1332         {
1333             switch (outtbtyp[ii])  /*switch to type of outtb object*/
1334             {
1335                 case SCSREAL_N    :
1336                     outtbdptr = (SCSREAL_COP *)outtbptr[ii]; /*double real matrix*/
1337                     sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1338                     C2F(dcopy)(&sszz, outtbdptr, &c__1, &outtbd[curouttbd], &c__1);
1339                     curouttbd += sszz;
1340                     break;
1341
1342                 case SCSCOMPLEX_N :
1343                     outtbdptr = (SCSCOMPLEX_COP *)outtbptr[ii]; /*double complex matrix*/
1344                     sszz = 2 * outtbsz[ii] * outtbsz[ii + nlnk];
1345                     C2F(dcopy)(&sszz, outtbdptr, &c__1, &outtbd[curouttbd], &c__1);
1346                     curouttbd += sszz;
1347                     break;
1348
1349                 case SCSINT8_N    :
1350                     outtbcptr = (SCSINT8_COP *)outtbptr[ii];  /*int8*/
1351                     sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1352                     for (kk = 0; kk < sszz; kk++)
1353                     {
1354                         outtbc[curouttbc + kk] = (SCSINT8_COP)outtbcptr[kk];
1355                     }
1356                     curouttbc += sszz;
1357                     break;
1358
1359                 case SCSINT16_N   :
1360                     outtbsptr = (SCSINT16_COP *)outtbptr[ii]; /*int16*/
1361                     sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1362                     for (kk = 0; kk < sszz; kk++)
1363                     {
1364                         outtbs[curouttbs + kk] = (SCSINT16_COP)outtbsptr[kk];
1365                     }
1366                     curouttbs += sszz;
1367                     break;
1368
1369                 case SCSINT32_N   :
1370                     outtblptr = (SCSINT32_COP *)outtbptr[ii];  /*int32*/
1371                     sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1372                     for (kk = 0; kk < sszz; kk++)
1373                     {
1374                         outtbl[curouttbl + kk] = (SCSINT32_COP)outtblptr[kk];
1375                     }
1376                     curouttbl += sszz;
1377                     break;
1378
1379                 case SCSUINT8_N   :
1380                     outtbucptr = (SCSUINT8_COP *)outtbptr[ii]; /*uint8*/
1381                     sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1382                     for (kk = 0; kk < sszz; kk++)
1383                     {
1384                         outtbuc[curouttbuc + kk] = (SCSUINT8_COP)outtbucptr[kk];
1385                     }
1386                     curouttbuc += sszz;
1387                     break;
1388
1389                 case SCSUINT16_N  :
1390                     outtbusptr = (SCSUINT16_COP *)outtbptr[ii]; /*uint16*/
1391                     sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1392                     for (kk = 0; kk < sszz; kk++)
1393                     {
1394                         outtbus[curouttbus + kk] = (SCSUINT16_COP)outtbusptr[kk];
1395                     }
1396                     curouttbus += sszz;
1397                     break;
1398
1399                 case SCSUINT32_N  :
1400                     outtbulptr = (SCSUINT32_COP *)outtbptr[ii]; /*uint32*/
1401                     sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1402                     for (kk = 0; kk < sszz; kk++)
1403                     {
1404                         outtbul[curouttbul + kk] = (SCSUINT32_COP)outtbulptr[kk];
1405                     }
1406                     curouttbul += sszz;
1407                     break;
1408
1409                 default  : /* Add a message here */
1410                     break;
1411             }
1412         }
1413     }
1414     *ierr = 20;
1415     freeouttbptr;
1416 } /* cosini_ */
1417
1418 /*--------------------------------------------------------------------------*/
1419 static void cossim(double *told)
1420 {
1421     /* System generated locals */
1422     int i3 = 0;
1423
1424     //** used for the [stop] button
1425     char* CommandToUnstack;
1426     static int CommandLength = 0;
1427     static int SeqSync = 0;
1428     static int one = 1;
1429
1430     /* Local variables */
1431     static scicos_flag flag__ = 0;
1432     static int ierr1 = 0;
1433     static int j = 0, k = 0;
1434     static double t = 0.;
1435     static int jj = 0;
1436     static double rhotmp = 0., tstop = 0.;
1437     static int kpo = 0, kev = 0;
1438     int Discrete_Jump = 0;
1439     int *jroot = NULL, *zcros = NULL;
1440     realtype reltol = 0., abstol = 0.;
1441     N_Vector y = NULL;
1442     void *ode_mem = NULL;
1443     int flag = 0, flagr = 0;
1444     int cnt = 0;
1445     /* Saving solver number */
1446     int solver = C2F(cmsolver).solver;
1447     /* Defining function pointers, for more readability */
1448     void(* ODEFree) (void**);
1449     int (* ODE) (void*, realtype, N_Vector, realtype*, int);
1450     int (* ODEReInit) (void*, realtype, N_Vector);
1451     int (* ODESetMaxStep) (void*, realtype);
1452     int (* ODESetStopTime) (void*, realtype);
1453     int (* ODEGetRootInfo) (void*, int*);
1454     int (* ODESStolerances) (void*, realtype, realtype);
1455     /* Generic flags for stop mode */
1456     int ODE_NORMAL   = 1;  /* ODE_NORMAL   = CV_NORMAL   = LS_NORMAL   = 1 */
1457     int ODE_ONE_STEP = 2;  /* ODE_ONE_STEP = CV_ONE_STEP = LS_ONE_STEP = 2 */
1458     switch (solver)
1459     {
1460         case LSodar_Dynamic:
1461             ODEFree = &LSodarFree;
1462             ODE = &LSodar;
1463             ODEReInit = &LSodarReInit;
1464             ODESetMaxStep = &LSodarSetMaxStep;
1465             ODESetStopTime = &LSodarSetStopTime;
1466             ODEGetRootInfo = &LSodarGetRootInfo;
1467             ODESStolerances = &LSodarSStolerances;
1468             break;
1469         case CVode_BDF_Newton:
1470         case CVode_BDF_Functional:
1471         case CVode_Adams_Newton:
1472         case CVode_Adams_Functional:
1473         case Dormand_Prince:
1474         case Runge_Kutta:
1475         case Implicit_Runge_Kutta:
1476             ODEFree = &CVodeFree;
1477             ODE = &CVode;
1478             ODEReInit = &CVodeReInit;
1479             ODESetMaxStep = &CVodeSetMaxStep;
1480             ODESetStopTime = &CVodeSetStopTime;
1481             ODEGetRootInfo = &CVodeGetRootInfo;
1482             ODESStolerances = &CVodeSStolerances;
1483             break;
1484         default: // Unknown solver number
1485             *ierr = 1000;
1486             return;
1487     }
1488
1489     jroot = NULL;
1490     if (ng > 0)
1491     {
1492         if ((jroot = MALLOC(sizeof(int) * ng)) == NULL )
1493         {
1494             *ierr = 10000;
1495             return;
1496         }
1497     }
1498
1499     for ( jj = 0 ; jj < ng ; jj++ )
1500     {
1501         jroot[jj] = 0 ;
1502     }
1503
1504     zcros = NULL;
1505     if (ng > 0)
1506     {
1507         if ((zcros = MALLOC(sizeof(int) * ng)) == NULL )
1508         {
1509             *ierr = 10000;
1510             if (ng > 0)
1511             {
1512                 FREE(jroot);
1513             }
1514             return;
1515         }
1516     }
1517
1518     reltol = (realtype) rtol;
1519     abstol = (realtype) Atol;  /* Ith(abstol,1) = realtype) Atol;*/
1520
1521     if (*neq > 0) /* Unfortunately CVODE does not work with NEQ==0 */
1522     {
1523         y = N_VNewEmpty_Serial(*neq);
1524         if (check_flag((void *)y, "N_VNewEmpty_Serial", 0))
1525         {
1526             *ierr = 10000;
1527             if (ng > 0)
1528             {
1529                 FREE(jroot);
1530             }
1531             if (ng > 0)
1532             {
1533                 FREE(zcros);
1534             }
1535             return;
1536         }
1537
1538         NV_DATA_S(y) = x;
1539
1540         ode_mem = NULL;
1541
1542         /* Set extension of Sundials for scicos */
1543         set_sundials_with_extension(TRUE);
1544
1545         switch (solver)
1546         {
1547             case LSodar_Dynamic:
1548                 ode_mem = LSodarCreate(neq, ng); /* Create the lsodar problem */
1549                 break;
1550             case CVode_BDF_Newton:
1551                 ode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
1552                 break;
1553             case CVode_BDF_Functional:
1554                 ode_mem = CVodeCreate(CV_BDF, CV_FUNCTIONAL);
1555                 break;
1556             case CVode_Adams_Newton:
1557                 ode_mem = CVodeCreate(CV_ADAMS, CV_NEWTON);
1558                 break;
1559             case CVode_Adams_Functional:
1560                 ode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);
1561                 break;
1562             case Dormand_Prince:
1563                 ode_mem = CVodeCreate(CV_DOPRI, CV_FUNCTIONAL);
1564                 break;
1565             case Runge_Kutta:
1566                 ode_mem = CVodeCreate(CV_ExpRK, CV_FUNCTIONAL);
1567                 break;
1568             case Implicit_Runge_Kutta:
1569                 ode_mem = CVodeCreate(CV_ImpRK, CV_FUNCTIONAL);
1570                 break;
1571         }
1572
1573         /*    ode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);*/
1574
1575         if (check_flag((void *)ode_mem, "CVodeCreate", 0))
1576         {
1577             *ierr = 10000;
1578             N_VDestroy_Serial(y);
1579             FREE(jroot);
1580             FREE(zcros);
1581             return;
1582         }
1583
1584         if (solver == LSodar_Dynamic)
1585         {
1586             flag = LSodarSetErrHandlerFn(ode_mem, SundialsErrHandler, NULL);
1587         }
1588         else
1589         {
1590             flag = CVodeSetErrHandlerFn(ode_mem, SundialsErrHandler, NULL);
1591         }
1592         if (check_flag(&flag, "CVodeSetErrHandlerFn", 1))
1593         {
1594             *ierr = 300 + (-flag);
1595             freeall
1596             return;
1597         }
1598
1599         if (solver == LSodar_Dynamic)
1600         {
1601             flag = LSodarInit(ode_mem, simblklsodar, T0, y);
1602         }
1603         else
1604         {
1605             flag = CVodeInit (ode_mem, simblk, T0, y);
1606         }
1607         if (check_flag(&flag, "CVodeInit", 1))
1608         {
1609             *ierr = 300 + (-flag);
1610             freeall
1611             return;
1612         }
1613
1614         flag = ODESStolerances(ode_mem, reltol, abstol);
1615         if (check_flag(&flag, "CVodeSStolerances", 1))
1616         {
1617             *ierr = 300 + (-flag);
1618             freeall
1619             return;
1620         }
1621
1622         if (solver == LSodar_Dynamic)
1623         {
1624             flag = LSodarRootInit(ode_mem, ng, grblklsodar);
1625         }
1626         else
1627         {
1628             flag = CVodeRootInit(ode_mem, ng, grblk);
1629         }
1630         if (check_flag(&flag, "CVodeRootInit", 1))
1631         {
1632             *ierr = 300 + (-flag);
1633             freeall
1634             return;
1635         }
1636
1637         /* Call CVDense to specify the CVDENSE dense linear solver, only for solvers needing CVode's Newton method */
1638         if (solver == CVode_BDF_Newton || solver == CVode_Adams_Newton)
1639         {
1640             flag = CVDense(ode_mem, *neq);
1641         }
1642         if (check_flag(&flag, "CVDense", 1))
1643         {
1644             *ierr = 300 + (-flag);
1645             freeall
1646             return;
1647         }
1648
1649         if (hmax > 0)
1650         {
1651             flag = ODESetMaxStep(ode_mem, (realtype) hmax);
1652             if (check_flag(&flag, "CVodeSetMaxStep", 1))
1653             {
1654                 *ierr = 300 + (-flag);
1655                 freeall;
1656                 return;
1657             }
1658         }
1659         /* Set the Jacobian routine to Jac (user-supplied)
1660         flag = CVDlsSetDenseJacFn(ode_mem, Jac);
1661         if (check_flag(&flag, "CVDlsSetDenseJacFn", 1)) return(1);  */
1662
1663     }/* testing if neq>0 */
1664
1665     /* Function Body */
1666     C2F(coshlt).halt = 0;
1667     *ierr = 0;
1668
1669     /*     initialization */
1670     C2F(realtimeinit)(told, &C2F(rtfactor).scale);
1671
1672     phase = 1;
1673     hot = 0;
1674
1675     jj = 0;
1676     for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
1677     {
1678         if (Blocks[C2F(curblk).kfun - 1].ng >= 1)
1679         {
1680             zcros[jj] = C2F(curblk).kfun;
1681             ++jj;
1682         }
1683     }
1684     /*     . ng >= jj required */
1685     if (jj != ng)
1686     {
1687         zcros[jj] = -1;
1688     }
1689     /*     initialization (propagation of constant blocks outputs) */
1690     idoit(told);
1691     if (*ierr != 0)
1692     {
1693         freeall;
1694         return;
1695     }
1696     /*--discrete zero crossings----dzero--------------------*/
1697     if (ng > 0) /* storing ZC signs just after a solver call*/
1698     {
1699         /*zdoit(told, g, x, x);*/
1700         zdoit(told, x, x, g);
1701         if (*ierr != 0)
1702         {
1703             freeall;
1704             return;
1705         }
1706         for (jj = 0; jj < ng; ++jj)
1707             if (g[jj] >= 0)
1708             {
1709                 jroot[jj] = 5;
1710             }
1711             else
1712             {
1713                 jroot[jj] = -5;
1714             }
1715     }
1716     /*--discrete zero crossings----dzero--------------------*/
1717
1718     /*     main loop on time */
1719     while (*told < *tf)
1720     {
1721         while (ismenu()) //** if the user has done something, do the actions
1722         {
1723             int ierr2 = 0;
1724             int iUnused;
1725             GetCommand(&CommandToUnstack, &SeqSync, &iUnused, NONE); //** get to the action
1726             CommandLength = (int)strlen(CommandToUnstack);
1727             //syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
1728             FREE(CommandToUnstack);
1729         }
1730         if (C2F(coshlt).halt != 0)
1731         {
1732             if (C2F(coshlt).halt == 2)
1733             {
1734                 *told = *tf;    /* end simulation */
1735             }
1736             C2F(coshlt).halt = 0;
1737             freeall;
1738             return;
1739         }
1740         if (*pointi == 0)
1741         {
1742             t = *tf;
1743         }
1744         else
1745         {
1746             t = tevts[*pointi];
1747         }
1748         if (fabs(t - *told) < ttol)
1749         {
1750             t = *told;
1751             /*     update output part */
1752         }
1753         if (*told > t)
1754         {
1755             /*     !  scheduling problem */
1756             *ierr = 1;
1757             freeall;
1758             return;
1759         }
1760         if (*told != t)
1761         {
1762             if (xptr[nblk + 1] == 1)
1763             {
1764                 /*     .     no continuous state */
1765                 if (*told + deltat + ttol > t)
1766                 {
1767                     *told = t;
1768                 }
1769                 else
1770                 {
1771                     *told += deltat;
1772                 }
1773                 /*     .     update outputs of 'c' type blocks with no continuous state */
1774                 if (*told >= *tf)
1775                 {
1776                     /*     .     we are at the end, update continuous part before leaving */
1777                     if (ncord > 0)
1778                     {
1779                         cdoit(told);
1780                         freeall;
1781                         return;
1782                     }
1783                 }
1784             }
1785             else
1786             {
1787                 /*     integrate */
1788                 rhotmp = *tf + ttol;
1789                 if (*pointi != 0)
1790                 {
1791                     kpo = *pointi;
1792 L20:
1793                     if (critev[kpo] == 1)
1794                     {
1795                         rhotmp = tevts[kpo];
1796                         goto L30;
1797                     }
1798                     kpo = evtspt[kpo];
1799                     if (kpo != 0)
1800                     {
1801                         goto L20;
1802                     }
1803 L30:
1804                     if (rhotmp < tstop)
1805                     {
1806                         hot = 0;
1807                     }
1808                 }
1809                 tstop = rhotmp;
1810                 t = Min(*told + deltat, Min(t, *tf + ttol));
1811
1812                 if (ng > 0 &&  hot == 0 && nmod > 0)
1813                 {
1814                     zdoit(told, x, x, g);
1815                     if (*ierr != 0)
1816                     {
1817                         freeall;
1818                         return;
1819                     }
1820                 }
1821
1822                 if (hot == 0) /* hot==0 : cold restart*/
1823                 {
1824                     flag = ODESetStopTime(ode_mem, (realtype)tstop);  /* Setting the stop time*/
1825                     if (check_flag(&flag, "CVodeSetStopTime", 1))
1826                     {
1827                         *ierr = 300 + (-flag);
1828                         freeall;
1829                         return;
1830                     }
1831
1832                     flag = ODEReInit(ode_mem, (realtype)(*told), y);
1833                     if (check_flag(&flag, "CVodeReInit", 1))
1834                     {
1835                         *ierr = 300 + (-flag);
1836                         freeall;
1837                         return;
1838                     }
1839                 }
1840
1841                 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1842                 {
1843                     sciprint(_("****SUNDIALS.Cvode from: %f to %f hot= %d  \n"), *told, t, hot);
1844                 }
1845
1846                 /*--discrete zero crossings----dzero--------------------*/
1847                 /*--check for Dzeros after Mode settings or ddoit()----*/
1848                 Discrete_Jump = 0;
1849
1850                 if (ng > 0 && hot == 0)
1851                 {
1852                     zdoit(told, x, x, g);
1853                     if (*ierr != 0)
1854                     {
1855                         freeall;
1856                         return;
1857                     }
1858                     for (jj = 0; jj < ng; ++jj)
1859                     {
1860                         if ((g[jj] >= 0.0) && (jroot[jj] == -5))
1861                         {
1862                             Discrete_Jump = 1;
1863                             jroot[jj] = 1;
1864                         }
1865                         else if ((g[jj] < 0.0) && (jroot[jj] == 5))
1866                         {
1867                             Discrete_Jump = 1;
1868                             jroot[jj] = -1;
1869                         }
1870                         else
1871                         {
1872                             jroot[jj] = 0;
1873                         }
1874                     }
1875                 }
1876                 /*--discrete zero crossings----dzero--------------------*/
1877
1878                 if (Discrete_Jump == 0) /* if there was a dzero, its event should be activated*/
1879                 {
1880                     phase = 2;
1881                     flag = ODE(ode_mem, t, y, told, ODE_NORMAL);
1882                     if (*ierr != 0)
1883                     {
1884                         freeall;
1885                         return;
1886                     }
1887                     phase = 1;
1888                 }
1889                 else
1890                 {
1891                     flag = CV_ROOT_RETURN; /* in order to handle discrete jumps */
1892                 }
1893
1894                 /*     .     update outputs of 'c' type  blocks if we are at the end*/
1895                 if (*told >= *tf)
1896                 {
1897                     if (ncord > 0)
1898                     {
1899                         cdoit(told);
1900                         freeall;
1901                         return;
1902                     }
1903                 }
1904
1905                 if (flag >= 0)
1906                 {
1907                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1908                     {
1909                         sciprint(_("****SUNDIALS.Cvode reached: %f\n"), *told);
1910                     }
1911                     hot = 1;
1912                     cnt = 0;
1913                 }
1914                 else if ( flag == CV_TOO_MUCH_WORK ||  flag == CV_CONV_FAILURE || flag == CV_ERR_FAILURE)
1915                 {
1916                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1917                     {
1918                         sciprint(_("****SUNDIALS.Cvode: too much work at time=%g (stiff region, change RTOL and ATOL)\n"), *told);
1919                     }
1920                     hot = 0;
1921                     cnt++;
1922                     if (cnt > 5)
1923                     {
1924                         *ierr = 300 + (-flag);
1925                         freeall;
1926                         return;
1927                     }
1928                 }
1929                 else
1930                 {
1931                     if (flag < 0)
1932                     {
1933                         *ierr = 300 + (-flag);    /* raising errors due to internal errors, otherwise error due to flagr*/
1934                     }
1935                     freeall;
1936                     return;
1937                 }
1938
1939                 if (flag == CV_ZERO_DETACH_RETURN)
1940                 {
1941                     hot = 0;
1942                 };  /* new feature of sundials, detects zero-detaching */
1943
1944                 if (flag == CV_ROOT_RETURN)
1945                 {
1946                     /*     .        at a least one root has been found */
1947                     hot = 0;
1948                     if (Discrete_Jump == 0)
1949                     {
1950                         flagr = ODEGetRootInfo(ode_mem, jroot);
1951                         if (check_flag(&flagr, "CVodeGetRootInfo", 1))
1952                         {
1953                             *ierr = 300 + (-flagr);
1954                             freeall;
1955                             return;
1956                         }
1957                     }
1958                     /*     .        at a least one root has been found */
1959                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1960                     {
1961                         sciprint(_("root found at t=: %f\n"), *told);
1962                     }
1963                     /*     .        update outputs affecting ztyp blocks ONLY FOR OLD BLOCKS */
1964                     zdoit(told, x, xd, g);
1965                     if (*ierr != 0)
1966                     {
1967                         freeall;
1968                         return;
1969                     }
1970                     for (jj = 0; jj < ng; ++jj)
1971                     {
1972                         C2F(curblk).kfun = zcros[ jj];
1973                         if (C2F(curblk).kfun == -1)
1974                         {
1975                             break;
1976                         }
1977                         kev = 0;
1978
1979                         for (j = zcptr[C2F(curblk).kfun] - 1 ;
1980                                 j < zcptr[C2F(curblk).kfun + 1] - 1 ; ++j)
1981                         {
1982                             if (jroot[j] != 0)
1983                             {
1984                                 kev = 1;
1985                                 break;
1986                             }
1987                         }
1988                         /*   */
1989                         if (kev != 0)
1990                         {
1991                             Blocks[C2F(curblk).kfun - 1].jroot = &jroot[zcptr[C2F(curblk).kfun] - 1];
1992                             if (funtyp[C2F(curblk).kfun] > 0)
1993                             {
1994
1995                                 if (Blocks[C2F(curblk).kfun - 1].nevout > 0)
1996                                 {
1997                                     flag__ = 3;
1998                                     if (Blocks[C2F(curblk).kfun - 1].nx > 0)
1999                                     {
2000                                         Blocks[C2F(curblk).kfun - 1].x  = &x[xptr[C2F(curblk).kfun] - 1];
2001                                         Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
2002                                     }
2003                                     /* call corresponding block to determine output event (kev) */
2004                                     Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
2005                                     /*callf(told, xd, x, x,g,&flag__);*/
2006                                     callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
2007                                     if (flag__ < 0)
2008                                     {
2009                                         *ierr = 5 - flag__;
2010                                         freeall;
2011                                         return;
2012                                     }
2013                                     /*     .              update event agenda */
2014                                     for (k = 0; k < Blocks[C2F(curblk).kfun - 1].nevout; ++k)
2015                                     {
2016                                         if (Blocks[C2F(curblk).kfun - 1].evout[k] >= 0.)
2017                                         {
2018                                             i3 = k + clkptr[C2F(curblk).kfun] ;
2019                                             addevs(Blocks[C2F(curblk).kfun - 1].evout[k] + (*told), &i3, &ierr1);
2020                                             if (ierr1 != 0)
2021                                             {
2022                                                 /*     .                       nevts too small */
2023                                                 *ierr = 3;
2024                                                 freeall;
2025                                                 return;
2026                                             }
2027                                         }
2028                                     }
2029                                 }
2030                                 /*     .              update state */
2031                                 if (Blocks[C2F(curblk).kfun - 1].nx > 0)
2032                                 {
2033                                     /*     .              call corresponding block to update state */
2034                                     flag__ = 2;
2035                                     Blocks[C2F(curblk).kfun - 1].x  = &x[xptr[C2F(curblk).kfun] - 1];
2036                                     Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
2037                                     Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
2038                                     /*callf(told, xd, x, x,g,&flag__);*/
2039                                     callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
2040
2041                                     if (flag__ < 0)
2042                                     {
2043                                         *ierr = 5 - flag__;
2044                                         freeall;
2045                                         return;
2046                                     }
2047                                 }
2048                             }
2049                         }
2050                     }
2051                 }
2052             }
2053             /*--discrete zero crossings----dzero--------------------*/
2054             if (ng > 0) /* storing ZC signs just after a sundials call*/
2055             {
2056                 zdoit(told, x, x, g);
2057                 if (*ierr != 0)
2058                 {
2059                     freeall;
2060                     return;
2061                 }
2062                 for (jj = 0; jj < ng; ++jj)
2063                 {
2064                     if (g[jj] >= 0)
2065                     {
2066                         jroot[jj] = 5;
2067                     }
2068                     else
2069                     {
2070                         jroot[jj] = -5;
2071                     }
2072                 }
2073             }
2074             /*--discrete zero crossings----dzero--------------------*/
2075
2076             C2F(realtime)(told);
2077         }
2078         else
2079         {
2080             /*     .  t==told */
2081             if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2082             {
2083                 sciprint(_("Event: %d activated at t=%f\n"), *pointi, *told);
2084                 for (kev = 0; kev < nblk; kev++)
2085                 {
2086                     if (Blocks[kev].nmode > 0)
2087                     {
2088                         sciprint(_("mode of block %d=%d, "), kev, Blocks[kev].mode[0]);
2089                     }
2090                 }
2091                 sciprint(_("**mod**\n"));
2092             }
2093
2094             ddoit(told);
2095             if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2096             {
2097                 sciprint(_("End of activation\n"));
2098             }
2099             if (*ierr != 0)
2100             {
2101                 freeall;
2102                 return;
2103             }
2104
2105         }
2106         /*     end of main loop on time */
2107     }
2108     freeall;
2109 } /* cossim_ */
2110
2111 /*--------------------------------------------------------------------------*/
2112 static void cossimdaskr(double *told)
2113 {
2114     /* System generated locals */
2115     int i3;
2116     //** used for the [stop] button
2117     char* CommandToUnstack;
2118     static int CommandLength = 0;
2119     static int SeqSync = 0;
2120     static int one = 1;
2121
2122     /* Local variables */
2123     static scicos_flag flag__ = 0;
2124     static int ierr1 = 0;
2125     static int j = 0, k = 0;
2126     static double t = 0.;
2127     static int jj = 0;
2128     static double rhotmp = 0., tstop = 0.;
2129     static int kpo = 0, kev = 0;
2130
2131     int *jroot = NULL, *zcros = NULL;
2132     int *Mode_save = NULL;
2133     int Mode_change = 0;
2134
2135     int flag = 0, flagr = 0;
2136     N_Vector   yy = NULL, yp = NULL;
2137     realtype reltol = 0., abstol = 0.;
2138     int Discrete_Jump = 0;
2139     N_Vector IDx = NULL;
2140     realtype *scicos_xproperty = NULL;
2141     DlsMat TJacque = NULL;
2142
2143     void *dae_mem = NULL;
2144     UserData data = NULL;
2145     IDAMem copy_IDA_mem = NULL;
2146     int maxnj = 0, maxnit = 0, maxnh = 0;
2147     /*-------------------- Analytical Jacobian memory allocation ----------*/
2148     int  Jn = 0, Jnx = 0, Jno = 0, Jni = 0, Jactaille = 0;
2149     double uround = 0.;
2150     int cnt = 0, N_iters = 0;
2151     /* Saving solver number */
2152     int solver = C2F(cmsolver).solver;
2153     /* Flags for initial values calculation */
2154     int DAE_YA_YDP_INIT = 1;
2155     int DAE_Y_INIT      = 2;
2156     /* Defining function pointers, for more readability*/
2157     void(* DAEFree) (void**);
2158     int (* DAESolve) (void*, realtype, realtype*, N_Vector, N_Vector, int);
2159     int (* DAEReInit) (void*, realtype, N_Vector, N_Vector);
2160     int (* DAESetId) (void*, N_Vector);
2161     int (* DAECalcIC) (void*, int, realtype);
2162     int (* DAESetMaxStep) (void*, realtype);
2163     int (* DAESetUserData) (void*, void*);
2164     int (* DAESetStopTime) (void*, realtype);
2165     int (* DAEGetRootInfo) (void*, int*);
2166     int (* DAESStolerances) (void*, realtype, realtype);
2167     int (* DAEGetConsistentIC) (void*, N_Vector, N_Vector);
2168     int (* DAESetMaxNumSteps) (void*, long int);
2169     int (* DAESetMaxNumJacsIC) (void*, int);
2170     int (* DAESetMaxNumItersIC) (void*, int);
2171     int (* DAESetMaxNumStepsIC) (void*, int);
2172     int (* DAESetLineSearchOffIC) (void*, int);
2173     /* For DAEs, the generic flags for stop mode depend on the used solver */
2174     int DAE_NORMAL = 0, DAE_ONE_STEP = 0;
2175     DAE_NORMAL   = (solver == IDA_BDF_Newton) ? 1 : 0;  /* IDA_NORMAL   = 1, DDAS_NORMAL   = 0 */
2176     DAE_ONE_STEP = (solver == IDA_BDF_Newton) ? 2 : 1;  /* IDA_ONE_STEP = 2, DDAS_ONE_STEP = 1 */
2177     switch (solver)
2178     {
2179         case IDA_BDF_Newton:
2180             DAEFree = &IDAFree;
2181             DAESolve = &IDASolve;
2182             DAESetId = &IDASetId;
2183             DAEReInit = &IDAReInit;
2184             DAECalcIC = &IDACalcIC;
2185             DAESetMaxStep = &IDASetMaxStep;
2186             DAESetUserData = &IDASetUserData;
2187             DAESetStopTime = &IDASetStopTime;
2188             DAEGetRootInfo = &IDAGetRootInfo;
2189             DAESStolerances = &IDASStolerances;
2190             DAESetMaxNumSteps = &IDASetMaxNumSteps;
2191             DAEGetConsistentIC = &IDAGetConsistentIC;
2192             DAESetMaxNumJacsIC = &IDASetMaxNumJacsIC;
2193             DAESetMaxNumItersIC = &IDASetMaxNumItersIC;
2194             DAESetMaxNumStepsIC = &IDASetMaxNumStepsIC;
2195             DAESetLineSearchOffIC = &IDASetLineSearchOffIC;
2196             break;
2197         case DDaskr_BDF_Newton:
2198         case DDaskr_BDF_GMRes:
2199             DAEFree = &DDaskrFree;
2200             DAESolve = &DDaskrSolve;
2201             DAESetId = &DDaskrSetId;
2202             DAEReInit = &DDaskrReInit;
2203             DAECalcIC = &DDaskrCalcIC;
2204             DAESetMaxStep = &DDaskrSetMaxStep;
2205             DAESetUserData = &DDaskrSetUserData;
2206             DAESetStopTime = &DDaskrSetStopTime;
2207             DAEGetRootInfo = &DDaskrGetRootInfo;
2208             DAESStolerances = &DDaskrSStolerances;
2209             DAESetMaxNumSteps = &DDaskrSetMaxNumSteps;
2210             DAEGetConsistentIC = &DDaskrGetConsistentIC;
2211             DAESetMaxNumJacsIC = &DDaskrSetMaxNumJacsIC;
2212             DAESetMaxNumItersIC = &DDaskrSetMaxNumItersIC;
2213             DAESetMaxNumStepsIC = &DDaskrSetMaxNumStepsIC;
2214             DAESetLineSearchOffIC = &DDaskrSetLineSearchOffIC;
2215             break;
2216         default: // Unknown solver number
2217             *ierr = 1000;
2218             return;
2219     }
2220
2221     /* Set extension of Sundials for scicos */
2222     set_sundials_with_extension(TRUE);
2223
2224     // CI=1.0;   /* for function Get_Jacobian_ci */
2225     jroot = NULL;
2226     if (ng > 0)
2227     {
2228         if ((jroot = MALLOC(sizeof(int) * ng)) == NULL )
2229         {
2230             *ierr = 10000;
2231             return;
2232         }
2233     }
2234     for ( jj = 0 ; jj < ng ; jj++ )
2235     {
2236         jroot[jj] = 0 ;
2237     }
2238
2239     zcros = NULL;
2240     if (ng > 0)
2241     {
2242         if ((zcros = MALLOC(sizeof(int) * ng)) == NULL )
2243         {
2244             *ierr = 10000;
2245             if (ng > 0)
2246             {
2247                 FREE(jroot);
2248             }
2249             return;
2250         }
2251     }
2252
2253     Mode_save = NULL;
2254     if (nmod > 0)
2255     {
2256         if ((Mode_save = MALLOC(sizeof(int) * nmod)) == NULL )
2257         {
2258             *ierr = 10000;
2259             if (ng > 0)
2260             {
2261                 FREE(jroot);
2262             }
2263             if (ng > 0)
2264             {
2265                 FREE(zcros);
2266             }
2267             return;
2268         }
2269     }
2270
2271     reltol = (realtype) rtol;
2272     abstol = (realtype) Atol;  /*  Ith(abstol,1) = (realtype) Atol;*/
2273
2274     if (*neq > 0)
2275     {
2276         yy = NULL;
2277         yy = N_VNewEmpty_Serial(*neq);
2278         if (check_flag((void *)yy, "N_VNew_Serial", 0))
2279         {
2280             if (ng > 0)
2281             {
2282                 FREE(jroot);
2283             }
2284             if (ng > 0)
2285             {
2286                 FREE(zcros);
2287             }
2288             if (nmod > 0)
2289             {
2290                 FREE(Mode_save);
2291             }
2292         }
2293         NV_DATA_S(yy) = x;
2294
2295         yp = NULL;
2296         yp = N_VNewEmpty_Serial(*neq);
2297         if (check_flag((void *)yp, "N_VNew_Serial", 0))
2298         {
2299             if (*neq > 0)
2300             {
2301                 N_VDestroy_Serial(yy);
2302             }
2303             if (ng > 0)
2304             {
2305                 FREE(jroot);
2306             }
2307             if (ng > 0)
2308             {
2309                 FREE(zcros);
2310             }
2311             if (nmod > 0)
2312             {
2313                 FREE(Mode_save);
2314             }
2315             return;
2316         }
2317         NV_DATA_S(yp) = xd;
2318
2319         IDx = NULL;
2320         IDx = N_VNew_Serial(*neq);
2321         if (check_flag((void *)IDx, "N_VNew_Serial", 0))
2322         {
2323             *ierr = 10000;
2324             if (*neq > 0)
2325             {
2326                 N_VDestroy_Serial(yp);
2327             }
2328             if (*neq > 0)
2329             {
2330                 N_VDestroy_Serial(yy);
2331             }
2332             if (ng > 0)
2333             {
2334                 FREE(jroot);
2335             }
2336             if (ng > 0)
2337             {
2338                 FREE(zcros);
2339             }
2340             if (nmod > 0)
2341             {
2342                 FREE(Mode_save);
2343             }
2344             return;
2345         }
2346
2347         /* Call the Create and Init functions to initialize DAE memory */
2348         dae_mem = NULL;
2349         if (solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes)
2350         {
2351             dae_mem = DDaskrCreate(neq, ng, solver);
2352         }
2353         else
2354         {
2355             dae_mem = IDACreate();
2356         }
2357         if (check_flag((void *)dae_mem, "IDACreate", 0))
2358         {
2359             if (*neq > 0)
2360             {
2361                 N_VDestroy_Serial(IDx);
2362             }
2363             if (*neq > 0)
2364             {
2365                 N_VDestroy_Serial(yp);
2366             }
2367             if (*neq > 0)
2368             {
2369                 N_VDestroy_Serial(yy);
2370             }
2371             if (ng > 0)
2372             {
2373                 FREE(jroot);
2374             }
2375             if (ng > 0)
2376             {
2377                 FREE(zcros);
2378             }
2379             if (nmod > 0)
2380             {
2381                 FREE(Mode_save);
2382             }
2383             return;
2384         }
2385         if (C2F(cmsolver).solver == 100)
2386         {
2387             copy_IDA_mem = (IDAMem) dae_mem;
2388         }
2389
2390         if (solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes)
2391         {
2392             flag = DDaskrSetErrHandlerFn(dae_mem, SundialsErrHandler, NULL);
2393         }
2394         else
2395         {
2396             flag = IDASetErrHandlerFn(dae_mem, SundialsErrHandler, NULL);
2397         }
2398         if (check_flag(&flag, "IDASetErrHandlerFn", 1))
2399         {
2400             *ierr = 200 + (-flag);
2401             if (*neq > 0)
2402             {
2403                 DAEFree(&dae_mem);
2404             }
2405             if (*neq > 0)
2406             {
2407                 N_VDestroy_Serial(IDx);
2408             }
2409             if (*neq > 0)
2410             {
2411                 N_VDestroy_Serial(yp);
2412             }
2413             if (*neq > 0)
2414             {
2415                 N_VDestroy_Serial(yy);
2416             }
2417             if (ng > 0)
2418             {
2419                 FREE(jroot);
2420             }
2421             if (ng > 0)
2422             {
2423                 FREE(zcros);
2424             }
2425             if (nmod > 0)
2426             {
2427                 FREE(Mode_save);
2428             }
2429             return;
2430         }
2431
2432         if (solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes)
2433         {
2434             flag = DDaskrInit(dae_mem, simblkddaskr, T0, yy, yp, jacpsol, psol);
2435         }
2436         else
2437         {
2438             flag = IDAInit(dae_mem, simblkdaskr, T0, yy, yp);
2439         }
2440         if (check_flag(&flag, "IDAInit", 1))
2441         {
2442             *ierr = 200 + (-flag);
2443             if (*neq > 0)
2444             {
2445                 DAEFree(&dae_mem);
2446             }
2447             if (*neq > 0)
2448             {
2449                 N_VDestroy_Serial(IDx);
2450             }
2451             if (*neq > 0)
2452             {
2453                 N_VDestroy_Serial(yp);
2454             }
2455             if (*neq > 0)
2456             {
2457                 N_VDestroy_Serial(yy);
2458             }
2459             if (ng > 0)
2460             {
2461                 FREE(jroot);
2462             }
2463             if (ng > 0)
2464             {
2465                 FREE(zcros);
2466             }
2467             if (nmod > 0)
2468             {
2469                 FREE(Mode_save);
2470             }
2471             return;
2472         }
2473
2474         flag = DAESStolerances(dae_mem, reltol, abstol);
2475         if (check_flag(&flag, "IDASStolerances", 1))
2476         {
2477             *ierr = 200 + (-flag);
2478             if (*neq > 0)
2479             {
2480                 DAEFree(&dae_mem);
2481             }
2482             if (*neq > 0)
2483             {
2484                 N_VDestroy_Serial(IDx);
2485             }
2486             if (*neq > 0)
2487             {
2488                 N_VDestroy_Serial(yp);
2489             }
2490             if (*neq > 0)
2491             {
2492                 N_VDestroy_Serial(yy);
2493             }
2494             if (ng > 0)
2495             {
2496                 FREE(jroot);
2497             }
2498             if (ng > 0)
2499             {
2500                 FREE(zcros);
2501             }
2502             if (nmod > 0)
2503             {
2504                 FREE(Mode_save);
2505             }
2506             return;
2507         }
2508
2509         if (solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes)
2510         {
2511             flag = DDaskrRootInit(dae_mem, ng, grblkddaskr);
2512         }
2513         else
2514         {
2515             flag = IDARootInit(dae_mem, ng, grblkdaskr);
2516         }
2517         if (check_flag(&flag, "IDARootInit", 1))
2518         {
2519             *ierr = 200 + (-flag);
2520             if (*neq > 0)
2521             {
2522                 DAEFree(&dae_mem);
2523             }
2524             if (*neq > 0)
2525             {
2526                 N_VDestroy_Serial(IDx);
2527             }
2528             if (*neq > 0)
2529             {
2530                 N_VDestroy_Serial(yp);
2531             }
2532             if (*neq > 0)
2533             {
2534                 N_VDestroy_Serial(yy);
2535             }
2536             if (ng > 0)
2537             {
2538                 FREE(jroot);
2539             }
2540             if (ng > 0)
2541             {
2542                 FREE(zcros);
2543             }
2544             if (nmod > 0)
2545             {
2546                 FREE(Mode_save);
2547             }
2548             return;
2549         }
2550
2551         if (solver == IDA_BDF_Newton)
2552         {
2553             flag = IDADense(dae_mem, *neq);
2554         }
2555         if (check_flag(&flag, "IDADense", 1))
2556         {
2557             *ierr = 200 + (-flag);
2558             if (*neq > 0)if (solver == IDA_BDF_Newton)
2559                 {
2560                     IDAFree(&dae_mem);
2561                 }
2562             if (*neq > 0)
2563             {
2564                 N_VDestroy_Serial(IDx);
2565             }
2566             if (*neq > 0)
2567             {
2568                 N_VDestroy_Serial(yp);
2569             }
2570             if (*neq > 0)
2571             {
2572                 N_VDestroy_Serial(yy);
2573             }
2574             if (ng > 0)
2575             {
2576                 FREE(jroot);
2577             }
2578             if (ng > 0)
2579             {
2580                 FREE(zcros);
2581             }
2582             if (nmod > 0)
2583             {
2584                 FREE(Mode_save);
2585             }
2586             return;
2587         }
2588
2589         data = NULL;
2590         if ((data = (UserData) MALLOC(sizeof(*data))) == NULL)
2591         {
2592             *ierr = 10000;
2593             if (*neq > 0)
2594             {
2595                 IDAFree(&dae_mem);
2596             }
2597             if (*neq > 0)
2598             {
2599                 N_VDestroy_Serial(IDx);
2600             }
2601             if (*neq > 0)
2602             {
2603                 N_VDestroy_Serial(yp);
2604             }
2605             if (*neq > 0)
2606             {
2607                 N_VDestroy_Serial(yy);
2608             }
2609             if (ng > 0)
2610             {
2611                 FREE(jroot);
2612             }
2613             if (ng > 0)
2614             {
2615                 FREE(zcros);
2616             }
2617             if (nmod > 0)
2618             {
2619                 FREE(Mode_save);
2620             }
2621             return;
2622         }
2623         data->dae_mem = dae_mem;
2624         data->ewt   = NULL;
2625         data->iwork = NULL;
2626         data->rwork = NULL;
2627         data->gwork = NULL;
2628
2629         data->ewt   = N_VNew_Serial(*neq);
2630         if (check_flag((void *)data->ewt, "N_VNew_Serial", 0))
2631         {
2632             *ierr = 200 + (-flag);
2633             if (*neq > 0)
2634             {
2635                 FREE(data);
2636             }
2637             if (*neq > 0)
2638             {
2639                 DAEFree(&dae_mem);
2640             }
2641             if (*neq > 0)
2642             {
2643                 N_VDestroy_Serial(IDx);
2644             }
2645             if (*neq > 0)
2646             {
2647                 N_VDestroy_Serial(yp);
2648             }
2649             if (*neq > 0)
2650             {
2651                 N_VDestroy_Serial(yy);
2652             }
2653             if (ng > 0)
2654             {
2655                 FREE(jroot);
2656             }
2657             if (ng > 0)
2658             {
2659                 FREE(zcros);
2660             }
2661             if (nmod > 0)
2662             {
2663                 FREE(Mode_save);
2664             }
2665             return;
2666         }
2667         if ( ng > 0 )
2668         {
2669             if ((data->gwork = (double *) MALLOC(ng * sizeof(double))) == NULL)
2670             {
2671                 if (*neq > 0)
2672                 {
2673                     N_VDestroy_Serial(data->ewt);
2674                 }
2675                 if (*neq > 0)
2676                 {
2677                     FREE(data);
2678                 }
2679                 if (*neq > 0)
2680                 {
2681                     DAEFree(&dae_mem);
2682                 }
2683                 if (*neq > 0)
2684                 {
2685                     N_VDestroy_Serial(IDx);
2686                 }
2687                 if (*neq > 0)
2688                 {
2689                     N_VDestroy_Serial(yp);
2690                 }
2691                 if (*neq > 0)
2692                 {
2693                     N_VDestroy_Serial(yy);
2694                 }
2695                 if (ng > 0)
2696                 {
2697                     FREE(jroot);
2698                 }
2699                 if (ng > 0)
2700                 {
2701                     FREE(zcros);
2702                 }
2703                 if (nmod > 0)
2704                 {
2705                     FREE(Mode_save);
2706                 }
2707                 return;
2708             }
2709         }
2710         /*Jacobian_Flag=0; */
2711         if (AJacobian_block > 0) /* set by the block with A-Jac in flag-4 using Set_Jacobian_flag(1); */
2712         {
2713             Jn = *neq;
2714             Jnx = Blocks[AJacobian_block - 1].nx;
2715             Jno = Blocks[AJacobian_block - 1].nout;
2716             Jni = Blocks[AJacobian_block - 1].nin;
2717         }
2718         else
2719         {
2720             Jn = *neq;
2721             Jnx = 0;
2722             Jno = 0;
2723             Jni = 0;
2724         }
2725         Jactaille = 3 * Jn + (Jn + Jni) * (Jn + Jno) + Jnx * (Jni + 2 * Jn + Jno) + (Jn - Jnx) * (2 * (Jn - Jnx) + Jno + Jni) + 2 * Jni * Jno;
2726
2727         if ((data->rwork = (double *) MALLOC(Jactaille * sizeof(double))) == NULL)
2728         {
2729             if ( ng > 0 )
2730             {
2731                 FREE(data->gwork);
2732             }
2733             if (*neq > 0)
2734             {
2735                 N_VDestroy_Serial(data->ewt);
2736             }
2737             if (*neq > 0)
2738             {
2739                 FREE(data);
2740             }
2741             if (*neq > 0)
2742             {
2743                 DAEFree(&dae_mem);
2744             }
2745             if (*neq > 0)
2746             {
2747                 N_VDestroy_Serial(IDx);
2748             }
2749             if (*neq > 0)
2750             {
2751                 N_VDestroy_Serial(yp);
2752             }
2753             if (*neq > 0)
2754             {
2755                 N_VDestroy_Serial(yy);
2756             }
2757             if (ng > 0)
2758             {
2759                 FREE(jroot);
2760             }
2761             if (ng > 0)
2762             {
2763                 FREE(zcros);
2764             }
2765             if (nmod > 0)
2766             {
2767                 FREE(Mode_save);
2768             }
2769             *ierr = 10000;
2770             return;
2771         }
2772
2773         if (solver == IDA_BDF_Newton)
2774         {
2775             flag = IDADlsSetDenseJacFn(dae_mem, Jacobians);
2776         }
2777         if (check_flag(&flag, "IDADlsSetDenseJacFn", 1))
2778         {
2779             *ierr = 200 + (-flag);
2780             freeallx
2781             return;
2782         }
2783
2784         TJacque = (DlsMat) NewDenseMat(*neq, *neq);
2785
2786         flag = DAESetUserData(dae_mem, data);
2787         if (check_flag(&flag, "IDASetUserData", 1))
2788         {
2789             *ierr = 200 + (-flag);
2790             freeallx
2791             return;
2792         }
2793
2794         if (hmax > 0)
2795         {
2796             flag = DAESetMaxStep(dae_mem, (realtype) hmax);
2797             if (check_flag(&flag, "IDASetMaxStep", 1))
2798             {
2799                 *ierr = 200 + (-flag);
2800                 freeallx
2801                 return;
2802             }
2803         }
2804
2805         maxnj = 100; /* setting the maximum number of Jacobian evaluations during a Newton step */
2806         flag = DAESetMaxNumJacsIC(dae_mem, maxnj);
2807         if (check_flag(&flag, "IDASetMaxNumJacsIC", 1))
2808         {
2809             *ierr = 200 + (-flag);
2810             freeallx
2811             return;
2812         }
2813
2814         maxnit = 10; /* setting the maximum number of Newton iterations in any attempt to solve CIC */
2815         if (C2F(cmsolver).solver == 102)
2816         {
2817             maxnit = 15;    /* By default, the Krylov max iterations should be 15 */
2818         }
2819         flag = DAESetMaxNumItersIC(dae_mem, maxnit);
2820         if (check_flag(&flag, "IDASetMaxNumItersIC", 1))
2821         {
2822             *ierr = 200 + (-flag);
2823             freeallx
2824             return;
2825         }
2826
2827         /* setting the maximum number of steps in an integration interval */
2828         maxnh = 2000;
2829         flag = DAESetMaxNumSteps(dae_mem, maxnh);
2830         if (check_flag(&flag, "IDASetMaxNumSteps", 1))
2831         {
2832             *ierr = 200 + (-flag);
2833             freeallx
2834             return;
2835         }
2836
2837     } /* testing if neq>0 */
2838
2839     uround = 1.0;
2840     do
2841     {
2842         uround = uround * 0.5;
2843     }
2844     while ( 1.0 + uround != 1.0);
2845     uround = uround * 2.0;
2846     SQuround = sqrt(uround);
2847     /* Function Body */
2848
2849     C2F(coshlt).halt = 0;
2850     *ierr = 0;
2851     /*     hot = .false. */
2852     phase = 1;
2853     hot = 0;
2854
2855     /*     initialization */
2856     C2F(realtimeinit)(told, &C2F(rtfactor).scale);
2857     /*     ATOL and RTOL are scalars */
2858
2859     jj = 0;
2860     for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
2861     {
2862         if (Blocks[C2F(curblk).kfun - 1].ng >= 1)
2863         {
2864             zcros[jj] = C2F(curblk).kfun;
2865             ++jj;
2866         }
2867     }
2868     /*     . ng >= jj required */
2869     if (jj != ng)
2870     {
2871         zcros[jj] = -1;
2872     }
2873     /*     initialization (propagation of constant blocks outputs) */
2874     idoit(told);
2875     if (*ierr != 0)
2876     {
2877         freeallx;
2878         return;
2879     }
2880
2881     /*--discrete zero crossings----dzero--------------------*/
2882     if (ng > 0) /* storing ZC signs just after a solver call*/
2883     {
2884         zdoit(told, x, x, g);
2885         if (*ierr != 0)
2886         {
2887             freeallx;
2888             return;
2889         }
2890         for (jj = 0; jj < ng; ++jj)
2891             if (g[jj] >= 0)
2892             {
2893                 jroot[jj] = 5;
2894             }
2895             else
2896             {
2897                 jroot[jj] = -5;
2898             }
2899     }
2900     /*     main loop on time */
2901     while (*told < *tf)
2902     {
2903         while (ismenu()) //** if the user has done something, do the actions
2904         {
2905             int ierr2 = 0;
2906             int iUnused;
2907             GetCommand(&CommandToUnstack, &SeqSync, &iUnused, NONE); //** get to the action
2908             CommandLength = (int)strlen(CommandToUnstack);
2909             //syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
2910             FREE(CommandToUnstack);
2911         }
2912         if (C2F(coshlt).halt != 0)
2913         {
2914             if (C2F(coshlt).halt == 2)
2915             {
2916                 *told = *tf;    /* end simulation */
2917             }
2918             C2F(coshlt).halt = 0;
2919             freeallx;
2920             return;
2921         }
2922         if (*pointi == 0)
2923         {
2924             t = *tf;
2925         }
2926         else
2927         {
2928             t = tevts[*pointi];
2929         }
2930         if (fabs(t - *told) < ttol)
2931         {
2932             t = *told;
2933             /*     update output part */
2934         }
2935         if (*told > t)
2936         {
2937             /*     !  scheduling problem */
2938             *ierr = 1;
2939             freeallx;
2940             return;
2941         }
2942         if (*told != t)
2943         {
2944             if (xptr[nblk + 1] == 1)
2945             {
2946                 /*     .     no continuous state */
2947                 if (*told + deltat + ttol > t)
2948                 {
2949                     *told = t;
2950                 }
2951                 else
2952                 {
2953                     *told += deltat;
2954                 }
2955                 /*     .     update outputs of 'c' type blocks with no continuous state */
2956                 if (*told >= *tf)
2957                 {
2958                     /*     .     we are at the end, update continuous part before leaving */
2959                     cdoit(told);
2960                     freeallx;
2961                     return;
2962                 }
2963             }
2964             else
2965             {
2966                 rhotmp = *tf + ttol;
2967                 if (*pointi != 0)
2968                 {
2969                     kpo = *pointi;
2970 L20:
2971                     if (critev[kpo] == 1)
2972                     {
2973                         rhotmp = tevts[kpo];
2974                         goto L30;
2975                     }
2976                     kpo = evtspt[kpo];
2977                     if (kpo != 0)
2978                     {
2979                         goto L20;
2980                     }
2981 L30:
2982                     if (rhotmp < tstop)
2983                     {
2984                         hot = 0;/* Cold-restart the solver if the new TSTOP isn't beyong the previous one*/
2985                     }
2986                 }
2987                 tstop = rhotmp;
2988                 t = Min(*told + deltat, Min(t, *tf + ttol));
2989
2990                 if (hot == 0)  /* CIC calculation when hot==0 */
2991                 {
2992
2993                     /* Setting the stop time*/
2994                     flag = DAESetStopTime(dae_mem, (realtype)tstop);
2995                     if (check_flag(&flag, "IDASetStopTime", 1))
2996                     {
2997                         *ierr = 200 + (-flag);
2998                         freeallx;
2999                         return;
3000                     }
3001
3002                     if (ng > 0 && nmod > 0)
3003                     {
3004                         phase = 1;
3005                         zdoit(told, x, xd, g);
3006                         if (*ierr != 0)
3007                         {
3008                             freeallx;
3009                             return;
3010                         }
3011                     }
3012
3013                     /*----------ID setting/checking------------*/
3014                     N_VConst(ONE, IDx); /* Initialize id to 1's. */
3015                     scicos_xproperty = NV_DATA_S(IDx);
3016                     reinitdoit(told);
3017                     if (*ierr > 0)
3018                     {
3019                         freeallx;
3020                         return;
3021                     }
3022                     for (jj = 0; jj < *neq; jj++)
3023                     {
3024                         if (xprop[jj] ==  1)
3025                         {
3026                             scicos_xproperty[jj] = ONE;
3027                         }
3028                         if (xprop[jj] == -1)
3029                         {
3030                             scicos_xproperty[jj] = ZERO;
3031                         }
3032                     }
3033                     /* CI=0.0;CJ=100.0; // for functions Get_Jacobian_ci and Get_Jacobian_cj
3034                     Jacobians(*neq, (realtype) (*told), yy, yp, bidon, (realtype) CJ, data, TJacque, tempv1, tempv2, tempv3);
3035                     for (jj=0;jj<*neq;jj++){
3036                     Jacque_col=DENSE_COL(TJacque,jj);
3037                     CI=ZERO;
3038                     for (kk=0;kk<*neq;kk++){
3039                     if ((Jacque_col[kk]-Jacque_col[kk]!=0)) {
3040                     CI=-ONE;
3041                     break;
3042                     }else{
3043                     if (Jacque_col[kk]!=0){
3044                     CI=ONE;
3045                     break;
3046                     }
3047                     }
3048                     }
3049                     if (CI>=ZERO){  scicos_xproperty[jj]=CI;}else{fprintf(stderr,"\nWarinng! Xproperties are not match for i=%d!",jj);}
3050                     } */
3051                     /* printf("\n"); for(jj=0;jj<*neq;jj++) { printf("x%d=%g ",jj,scicos_xproperty[jj]); }*/
3052
3053                     flag = DAESetId(dae_mem, IDx);
3054                     if (check_flag(&flag, "IDASetId", 1))
3055                     {
3056                         *ierr = 200 + (-flag);
3057                         freeallx
3058                         return;
3059                     }
3060                     // CI=1.0;  // for function Get_Jacobian_ci
3061                     /*--------------------------------------------*/
3062                     // maxnj=100; /* setting the maximum number of Jacobian evaluation during a Newton step */
3063                     // flag=DAESetMaxNumJacsIC(dae_mem, maxnj);
3064                     // if (check_flag(&flag, "IDASetMaxNumJacsIC", 1)) {
3065                     //   *ierr=200+(-flag);
3066                     //   freeallx;
3067                     //   return;
3068                     // };
3069                     // flag=DAESetLineSearchOffIC(dae_mem,FALSE);  /* (def=false)  */
3070                     // if (check_flag(&flag, "IDASetLineSearchOffIC", 1)) {
3071                     //   *ierr=200+(-flag);
3072                     //   freeallx;
3073                     //   return;
3074                     // };
3075                     // flag=DAESetMaxNumItersIC(dae_mem, 10);/* (def=10) setting the maximum number of Newton iterations in any attempt to solve CIC */
3076                     // if (check_flag(&flag, "IDASetMaxNumItersIC", 1)) {
3077                     //   *ierr=200+(-flag);
3078                     //   freeallx;
3079                     //   return;
3080                     // };
3081
3082                     N_iters = 4 + nmod * 4;
3083                     for (j = 0; j <= N_iters; j++)
3084                     {
3085                         /* counter to reevaluate the
3086                                                                         modes in  mode->CIC->mode->CIC-> loop
3087                                                                         do it once in the absence of mode (nmod=0) */
3088                         /* updating the modes through Flag==9, Phase==1 */
3089
3090                         /* Serge Steer 29/06/2009 */
3091                         while (ismenu()) //** if the user has done something, do the actions
3092                         {
3093                             int ierr2 = 0;
3094                             int iUnused;
3095                             GetCommand(&CommandToUnstack, &SeqSync, &iUnused, NONE); //** get to the action
3096                             CommandLength = (int)strlen(CommandToUnstack);
3097                             //syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
3098                             FREE(CommandToUnstack);
3099                         }
3100                         if (C2F(coshlt).halt != 0)
3101                         {
3102                             if (C2F(coshlt).halt == 2)
3103                             {
3104                                 *told = *tf;    /* end simulation */
3105                             }
3106                             C2F(coshlt).halt = 0;
3107                             freeallx;
3108                             return;
3109                         }
3110
3111                         /* yy->PH */
3112                         flag = DAEReInit(dae_mem, (realtype)(*told), yy, yp);
3113                         if (check_flag(&flag, "IDAReInit", 1))
3114                         {
3115                             *ierr = 200 + (-flag);
3116                             freeallx;
3117                             return;
3118                         }
3119
3120                         phase = 2; /* IDACalcIC: PHI-> yy0: if (ok) yy0_cic-> PHI*/
3121                         if (C2F(cmsolver).solver == 100)
3122                         {
3123                             copy_IDA_mem->ida_kk = 1;
3124                         }
3125                         // the initial conditons y0 and yp0 do not satisfy the DAE
3126                         flagr = DAECalcIC(dae_mem, DAE_YA_YDP_INIT, (realtype)(t));
3127                         phase = 1;
3128                         flag = DAEGetConsistentIC(dae_mem, yy, yp); /* PHI->YY */
3129                         if (*ierr > 5)    /* *ierr>5 => singularity in block */
3130                         {
3131                             freeallx;
3132                             return;
3133                         }
3134
3135                         if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3136                         {
3137                             if (flagr >= 0)
3138                             {
3139                                 sciprint(_("**** SUNDIALS.IDA successfully initialized *****\n") );
3140                             }
3141                             else
3142                             {
3143                                 sciprint(_("**** SUNDIALS.IDA failed to initialize ->try again *****\n") );
3144                             }
3145                         }
3146                         /*-------------------------------------*/
3147                         /* saving the previous modes*/
3148                         for (jj = 0; jj < nmod; ++jj)
3149                         {
3150                             Mode_save[jj] = mod[jj];
3151                         }
3152                         if (ng > 0 && nmod > 0)
3153                         {
3154                             phase = 1;
3155                             zdoit(told, x, xd, g);
3156                             if (*ierr != 0)
3157                             {
3158                                 freeallx;
3159                                 return;
3160                             }
3161                         }
3162                         /*------------------------------------*/
3163                         Mode_change = 0;
3164                         for (jj = 0; jj < nmod; ++jj)
3165                         {
3166                             if (Mode_save[jj] != mod[jj])
3167                             {
3168                                 Mode_change = 1;
3169                                 break;
3170                             }
3171                         }
3172                         if (Mode_change == 0)
3173                         {
3174                             if (flagr >= 0 )
3175                             {
3176                                 break;  /*   if (flagr>=0) break;  else{ *ierr=200+(-flagr); freeallx;  return; }*/
3177                             }
3178                             else if (j >= (int)( N_iters / 2))
3179                             {
3180                                 /* DAESetMaxNumStepsIC(mem,10); */     /* maxnh (def=5) */
3181                                 DAESetMaxNumJacsIC(dae_mem, 10);      /* maxnj 100 (def=4)*/
3182                                 /* DAESetMaxNumItersIC(mem,100000); */ /* maxnit in IDANewtonIC (def=10) */
3183                                 DAESetLineSearchOffIC(dae_mem, TRUE); /* (def=false)  */
3184                                 /* DAESetNonlinConvCoefIC(mem,1.01);*/ /* (def=0.01-0.33*/
3185                                 flag = DAESetMaxNumItersIC(dae_mem, 1000);
3186                                 if (check_flag(&flag, "IDASetMaxNumItersIC", 1))
3187                                 {
3188                                     *ierr = 200 + (-flag);
3189                                     freeallx;
3190                                     return;
3191                                 };
3192                             }
3193                         }
3194                     }/* mode-CIC  counter*/
3195                     if (Mode_change == 1)
3196                     {
3197                         /* In this case, we try again by relaxing all modes and calling DAECalcIC again
3198                         /Masoud */
3199                         phase = 1;
3200                         if (C2F(cmsolver).solver == 100)
3201                         {
3202                             copy_IDA_mem->ida_kk = 1;
3203                         }
3204                         flagr = DAECalcIC(dae_mem, DAE_YA_YDP_INIT, (realtype)(t));
3205                         phase = 1;
3206                         flag = DAEGetConsistentIC(dae_mem, yy, yp); /* PHI->YY */
3207                         if ((flagr < 0) || (*ierr > 5)) /* *ierr>5 => singularity in block */
3208                         {
3209                             *ierr = 23;
3210                             freeallx;
3211                             return;
3212                         }
3213                     }
3214                     /*-----If flagr<0 the initialization solver has not converged-----*/
3215                     if (flagr < 0)
3216                     {
3217                         *ierr = 237;
3218                         freeallx;
3219                         return;
3220                     }
3221
3222                 } /* CIC calculation when hot==0 */
3223
3224                 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3225                 {
3226                     sciprint(_("****daskr from: %f to %f hot= %d  \n"), *told, t, hot);
3227                 }
3228
3229                 /*--discrete zero crossings----dzero--------------------*/
3230                 /*--check for Dzeros after Mode settings or ddoit()----*/
3231                 Discrete_Jump = 0;
3232                 if (ng > 0 && hot == 0)
3233                 {
3234                     zdoit(told, x, xd, g);
3235                     if (*ierr != 0)
3236                     {
3237                         freeallx;
3238                         return;
3239                     }
3240                     for (jj = 0; jj < ng; ++jj)
3241                     {
3242                         if ((g[jj] >= 0.0) && ( jroot[jj] == -5))
3243                         {
3244                             Discrete_Jump = 1;
3245                             jroot[jj] = 1;
3246                         }
3247                         else if ((g[jj] < 0.0) && ( jroot[jj] == 5))
3248                         {
3249                             Discrete_Jump = 1;
3250                             jroot[jj] = -1;
3251                         }
3252                         else
3253                         {
3254                             jroot[jj] = 0;
3255                         }
3256                     }
3257                 }
3258
3259                 /*--discrete zero crossings----dzero--------------------*/
3260                 if (Discrete_Jump == 0) /* if there was a dzero, its event should be activated*/
3261                 {
3262                     phase = 2;
3263                     flagr = DAESolve(dae_mem, t, told, yy, yp, DAE_NORMAL);
3264                     phase = 1;
3265                     if (*ierr != 0)
3266                     {
3267                         freeallx;
3268                         return;
3269                     }
3270                 }
3271                 else
3272                 {
3273                     flagr = IDA_ROOT_RETURN; /* in order to handle discrete jumps */
3274                 }
3275                 if (flagr >= 0)
3276                 {
3277                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3278                     {
3279                         sciprint(_("****SUNDIALS.Ida reached: %f\n"), *told);
3280                     }
3281                     hot = 1;
3282                     cnt = 0;
3283                 }
3284                 else if ( flagr == IDA_TOO_MUCH_WORK ||  flagr == IDA_CONV_FAIL || flagr == IDA_ERR_FAIL)
3285                 {
3286                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3287                     {
3288                         sciprint(_("**** SUNDIALS.Ida: too much work at time=%g (stiff region, change RTOL and ATOL)\n"), *told);
3289                     }
3290                     hot = 0;
3291                     cnt++;
3292                     if (cnt > 5)
3293                     {
3294                         *ierr = 200 + (-flagr);
3295                         freeallx;
3296                         return;
3297                     }
3298                 }
3299                 else
3300                 {
3301                     if (flagr < 0)
3302                     {
3303                         *ierr = 200 + (-flagr);    /* raising errors due to internal errors, otherwise error due to flagr*/
3304                     }
3305                     freeallx;
3306                     return;
3307                 }
3308
3309                 /*     update outputs of 'c' type  blocks if we are at the end*/
3310                 if (*told >= *tf)
3311                 {
3312                     cdoit(told);
3313                     freeallx;
3314                     return;
3315                 }
3316
3317                 if (flagr == IDA_ZERO_DETACH_RETURN)
3318                 {
3319                     hot = 0;
3320                 }; /* new feature of sundials, detects unmasking */
3321                 if (flagr == IDA_ROOT_RETURN)
3322                 {
3323                     /*     .        at least one root has been found */
3324                     hot = 0;
3325                     if (Discrete_Jump == 0)
3326                     {
3327                         flagr = DAEGetRootInfo(dae_mem, jroot);
3328                         if (check_flag(&flagr, "IDAGetRootInfo", 1))
3329                         {
3330                             *ierr = 200 + (-flagr);
3331                             freeallx;
3332                             return;
3333                         }
3334                     }
3335
3336                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3337                     {
3338                         sciprint(_("root found at t=: %f\n"), *told);
3339                     }
3340                     /*     .        update outputs affecting ztyp blocks  ONLY FOR OLD BLOCKS*/
3341                     zdoit(told, x, xd, g);
3342                     if (*ierr != 0)
3343                     {
3344                         freeallx;
3345                         return;
3346                     }
3347                     for (jj = 0; jj < ng; ++jj)
3348                     {
3349                         C2F(curblk).kfun = zcros[jj];
3350                         if (C2F(curblk).kfun == -1)
3351                         {
3352                             break;
3353                         }
3354                         kev = 0;
3355                         for (j = zcptr[C2F(curblk).kfun] - 1 ;
3356                                 j < zcptr[C2F(curblk).kfun + 1] - 1 ; ++j)
3357                         {
3358                             if (jroot[j] != 0)
3359                             {
3360                                 kev = 1;
3361                                 break;
3362                             }
3363                         }
3364                         if (kev != 0)
3365                         {
3366                             Blocks[C2F(curblk).kfun - 1].jroot = &jroot[zcptr[C2F(curblk).kfun] - 1];
3367                             if (funtyp[C2F(curblk).kfun] > 0)
3368                             {
3369                                 if (Blocks[C2F(curblk).kfun - 1].nevout > 0)
3370                                 {
3371                                     flag__ = 3;
3372                                     if (Blocks[C2F(curblk).kfun - 1].nx > 0)
3373                                     {
3374                                         Blocks[C2F(curblk).kfun - 1].x  = &x[xptr[C2F(curblk).kfun] - 1];
3375                                         Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
3376                                     }
3377                                     /*     call corresponding block to determine output event (kev) */
3378                                     Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
3379                                     /*callf(told, xd, x, x,g,&flag__);*/
3380                                     callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
3381                                     if (flag__ < 0)
3382                                     {
3383                                         *ierr = 5 - flag__;
3384                                         freeallx;
3385                                         return;
3386                                     }
3387                                     /*     update event agenda */
3388                                     for (k = 0; k < Blocks[C2F(curblk).kfun - 1].nevout; ++k)
3389                                     {
3390                                         if (Blocks[C2F(curblk).kfun - 1].evout[k] >= 0)
3391                                         {
3392                                             i3 = k + clkptr[C2F(curblk).kfun] ;
3393                                             addevs(Blocks[C2F(curblk).kfun - 1].evout[k] + (*told), &i3, &ierr1);
3394                                             if (ierr1 != 0)
3395                                             {
3396                                                 /*     .                       nevts too small */
3397                                                 *ierr = 3;
3398                                                 freeallx;
3399                                                 return;
3400                                             }
3401                                         }
3402                                     }
3403                                 }
3404                                 /* update state */
3405                                 if ((Blocks[C2F(curblk).kfun - 1].nx > 0) || (*Blocks[C2F(curblk).kfun - 1].work != NULL) )
3406                                 {
3407                                     /* call corresponding block to update state */
3408                                     flag__ = 2;
3409                                     if (Blocks[C2F(curblk).kfun - 1].nx > 0)
3410                                     {
3411                                         Blocks[C2F(curblk).kfun - 1].x  = &x[xptr[C2F(curblk).kfun] - 1];
3412                                         Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
3413                                     }
3414                                     Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
3415
3416                                     Blocks[C2F(curblk).kfun - 1].xprop = &xprop[-1 + xptr[C2F(curblk).kfun]];
3417                                     /*callf(told, xd, x, x,g,&flag__);*/
3418                                     callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
3419
3420                                     if (flag__ < 0)
3421                                     {
3422                                         *ierr = 5 - flag__;
3423                                         freeallx;
3424                                         return;
3425                                     }
3426                                     for (j = 0; j < *neq; j++) /* Adjust xprop for IDx */
3427                                     {
3428                                         if (xprop[j] ==  1)
3429                                         {
3430                                             scicos_xproperty[j] = ONE;
3431                                         }
3432                                         if (xprop[j] == -1)
3433                                         {
3434                                             scicos_xproperty[j] = ZERO;
3435                                         }
3436                                     }
3437                                 }
3438                             }
3439                         }
3440                     }
3441                 }
3442                 /* Serge Steer 29/06/2009 */
3443                 while (ismenu()) //** if the user has done something, do the actions
3444                 {
3445                     int ierr2 = 0;
3446                     int iUnused;
3447                     GetCommand(&CommandToUnstack, &SeqSync, &iUnused, NONE); //** get to the action
3448                     CommandLength = (int)strlen(CommandToUnstack);
3449                     //syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
3450                     FREE(CommandToUnstack);
3451                 }
3452
3453                 if (C2F(coshlt).halt != 0)
3454                 {
3455                     if (C2F(coshlt).halt == 2)
3456                     {
3457                         *told = *tf;    /* end simulation */
3458                     }
3459                     C2F(coshlt).halt = 0;
3460                     freeallx;
3461                     return;
3462                 }
3463                 /* if(*pointi!=0){
3464                 t=tevts[*pointi];
3465                 if(*told<t-ttol){
3466                 cdoit(told);
3467                 goto L15;
3468                 }
3469                 }else{
3470                 if(*told<*tf){
3471                 cdoit(told);
3472                 goto L15;
3473                 }
3474                 }*/
3475
3476                 /*--discrete zero crossings----dzero--------------------*/
3477                 if (ng > 0) /* storing ZC signs just after a ddaskr call*/
3478                 {
3479                     zdoit(told, x, xd, g);
3480                     if (*ierr != 0)
3481                     {
3482                         freeallx;
3483                         return;
3484                     }
3485                     for (jj = 0; jj < ng; ++jj)
3486                     {
3487                         if (g[jj] >= 0)
3488                         {
3489                             jroot[jj] = 5;
3490                         }
3491                         else
3492                         {
3493                             jroot[jj] = -5;
3494                         }
3495                     }
3496                 }
3497                 /*--discrete zero crossings----dzero--------------------*/
3498             }
3499             C2F(realtime)(told);
3500         }
3501         else
3502         {
3503             /*     .  t==told */
3504             if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3505             {
3506                 sciprint(_("Event: %d activated at t=%f\n"), *pointi, *told);
3507             }
3508
3509             ddoit(told);
3510             if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3511             {
3512                 sciprint(_("End of activation\n"));
3513             }
3514             if (*ierr != 0)
3515             {
3516                 freeallx;
3517                 return;
3518             }
3519         }
3520         /*     end of main loop on time */
3521     }
3522     freeallx;
3523 } /* cossimdaskr_ */
3524 /*--------------------------------------------------------------------------*/
3525 /* Subroutine cosend */
3526 static void cosend(double *told)
3527 {
3528     /* Local variables */
3529     static scicos_flag flag__ = 0;
3530
3531     static int kfune = 0;
3532
3533     /* Function Body */
3534     *ierr = 0;
3535     /*     loop on blocks */
3536     for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
3537     {
3538         flag__ = 5;
3539         Blocks[C2F(curblk).kfun - 1].nevprt = 0;
3540         if (funtyp[C2F(curblk).kfun] >= 0)
3541         {
3542             if (Blocks[C2F(curblk).kfun - 1].nx > 0)
3543             {
3544                 Blocks[C2F(curblk).kfun - 1].x  = &x[xptr[C2F(curblk).kfun] - 1];
3545                 Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
3546             }
3547             /*callf(told, xd, x, x,x,&flag__);*/
3548             callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
3549             if (flag__ < 0 && *ierr == 0)
3550             {
3551                 *ierr = 5 - flag__;
3552                 kfune = C2F(curblk).kfun;
3553             }
3554         }
3555     }
3556     if (*ierr != 0)
3557     {
3558         C2F(curblk).kfun = kfune;
3559         return;
3560     }
3561 } /* cosend_ */
3562 /*--------------------------------------------------------------------------*/
3563 /* callf */
3564 void callf(double *t, scicos_block *block, scicos_flag *flag)
3565 {
3566     double* args[SZ_SIZE];
3567     int sz[SZ_SIZE];
3568     double intabl[TB_SIZE];
3569     double outabl[TB_SIZE];
3570
3571     int ii = 0, in = 0, out = 0, ki = 0, ko = 0, no = 0, ni = 0, k = 0, j = 0;
3572     int szi = 0, flagi = 0;
3573     double *ptr_d = NULL;
3574
3575     /* function pointers type def */
3576     voidf loc ;
3577     ScicosF0 loc0;
3578     ScicosF loc1;
3579     /*  ScicosFm1 loc3;*/
3580     ScicosF2 loc2;
3581     ScicosF2z loc2z;
3582     ScicosFi loci1;
3583     ScicosFi2 loci2;
3584     ScicosFi2z loci2z;
3585     ScicosF4 loc4;
3586
3587     int solver = C2F(cmsolver).solver;
3588     int cosd   = C2F(cosdebug).cosd;
3589     /*int kf     = C2F(curblk).kfun;*/
3590     scicos_time     = *t;
3591     block_error     = (int*) flag;
3592
3593     /* debug block is never called */
3594     /*if (kf==(debug_block+1)) return;*/
3595     if (block->type == DEBUG_BLOCK)
3596     {
3597         return;
3598     }
3599
3600     /* flag 7 implicit initialization */
3601     flagi = (int) * flag;
3602     /* change flag to zero if flagi==7 for explicit block */
3603     if (flagi == 7 && block->type < 10000)
3604     {
3605         *flag = 0;
3606     }
3607
3608     /* display information for debugging mode */
3609     if (cosd > 1)
3610     {
3611         if (cosd != 3)
3612         {
3613             sciprint(_("block %d [%s] is called "), C2F(curblk).kfun, block->uid);
3614             sciprint(_("with flag %d "), *flag);
3615             sciprint(_("at time %f \n"), *t);
3616         }
3617         if (debug_block > -1)
3618         {
3619             if (cosd != 3)
3620             {
3621                 sciprint(_("Entering the block \n"));
3622             }
3623             call_debug_scicos(block, flag, flagi, debug_block);
3624             if (*flag < 0)
3625             {
3626                 return;    /* error in debug block */
3627             }
3628         }
3629     }
3630
3631     C2F(scsptr).ptr = block->scsptr;
3632
3633     /* get pointer of the function */
3634     loc = block->funpt;
3635
3636     /* continuous state */
3637     if ((solver == IDA_BDF_Newton || solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes) && block->type < 10000 && *flag == 0)
3638     {
3639         ptr_d = block->xd;
3640         block->xd  = block->res;
3641     }
3642
3643     /* switch loop */
3644     //sciprint("callf type=%d flag=%d\n",block->type,flagi);
3645     switch (block->type)
3646     {
3647         /*******************/
3648         /* function type 0 */
3649         /*******************/
3650         case 0 :
3651         {
3652             /* This is for compatibility */
3653             /* jroot is returned in g for old type */
3654             if (block->nevprt < 0)
3655             {
3656                 for (j = 0; j < block->ng; ++j)
3657                 {
3658                     block->g[j] = (double)block->jroot[j];
3659                 }
3660             }
3661
3662             /* concatenated entries and concatened outputs */
3663             /* catenate inputs if necessary */
3664             ni = 0;
3665             if (block->nin > 1)
3666             {
3667                 ki = 0;
3668                 for (in = 0; in < block->nin; in++)
3669                 {
3670                     szi = block->insz[in] * block->insz[in + block->nin];
3671                     for (ii = 0; ii < szi; ii++)
3672                     {
3673                         intabl[ki++] = *((double *)(block->inptr[in]) + ii);
3674                     }
3675                     ni = ni + szi;
3676                 }
3677                 args[0] = &(intabl[0]);
3678             }
3679             else
3680             {
3681                 if (block->nin == 0)
3682                 {
3683                     args[0] = NULL;
3684                 }
3685                 else
3686                 {
3687                     args[0] = (double *)(block->inptr[0]);
3688                     ni = block->insz[0] * block->insz[1];
3689                 }
3690             }
3691
3692             /* catenate outputs if necessary */
3693             no = 0;
3694             if (block->nout > 1)
3695             {
3696                 ko = 0;
3697                 for (out = 0; out < block->nout; out++)
3698                 {
3699                     szi = block->outsz[out] * block->outsz[out + block->nout];
3700                     for (ii = 0; ii < szi; ii++)
3701                     {
3702                         outabl[ko++] = *((double *)(block->outptr[out]) + ii);
3703                     }
3704                     no = no + szi;
3705                 }
3706                 args[1] = &(outabl[0]);
3707             }
3708             else
3709             {
3710                 if (block->nout == 0)
3711                 {
3712                     args[1] = NULL;
3713                 }
3714                 else
3715                 {
3716                     args[1] = (double *)(block->outptr[0]);
3717                     no = block->outsz[0] * block->outsz[1];
3718                 }
3719             }
3720
3721             loc0 = (ScicosF0) loc;
3722
3723             (*loc0)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3724                     block->z, &block->nz,
3725                     block->evout, &block->nevout, block->rpar, &block->nrpar,
3726                     block->ipar, &block->nipar, (double *)args[0], &ni,
3727                     (double *)args[1], &no);
3728
3729             /* split output vector on each port if necessary */
3730             if (block->nout > 1)
3731             {
3732                 ko = 0;
3733                 for (out = 0; out < block->nout; out++)
3734                 {
3735                     szi = block->outsz[out] * block->outsz[out + block->nout];
3736                     for (ii = 0; ii < szi; ii++)
3737                     {
3738                         *((double *)(block->outptr[out]) + ii) = outabl[ko++];
3739                     }
3740                 }
3741             }
3742
3743             /* adjust values of output register */
3744             for (in = 0; in < block->nevout; ++in)
3745             {
3746                 block->evout[in] = block->evout[in] - *t;
3747             }
3748
3749             break;
3750         }
3751
3752         /*******************/
3753         /* function type 1 */
3754         /*******************/
3755         case 1 :
3756         {
3757             /* This is for compatibility */
3758             /* jroot is returned in g for old type */
3759             if (block->nevprt < 0)
3760             {
3761                 for (j = 0; j < block->ng; ++j)
3762                 {
3763                     block->g[j] = (double)block->jroot[j];
3764                 }
3765             }
3766
3767             /* one entry for each input or output */
3768             for (in = 0 ; in < block->nin ; in++)
3769             {
3770                 args[in] = block->inptr[in];
3771                 sz[in] = block->insz[in];
3772             }
3773             for (out = 0; out < block->nout; out++)
3774             {
3775                 args[in + out] = block->outptr[out];
3776                 sz[in + out] = block->outsz[out];
3777             }
3778             /* with zero crossing */
3779             if (block->ztyp > 0)
3780             {
3781                 args[block->nin + block->nout] = block->g;
3782                 sz[block->nin + block->nout] = block->ng;
3783             }
3784
3785             loc1 = (ScicosF) loc;
3786
3787             (*loc1)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3788                     block->z, &block->nz,
3789                     block->evout, &block->nevout, block->rpar, &block->nrpar,
3790                     block->ipar, &block->nipar,
3791                     (double *)args[0], &sz[0],
3792                     (double *)args[1], &sz[1], (double *)args[2], &sz[2],
3793                     (double *)args[3], &sz[3], (double *)args[4], &sz[4],
3794                     (double *)args[5], &sz[5], (double *)args[6], &sz[6],
3795                     (double *)args[7], &sz[7], (double *)args[8], &sz[8],
3796                     (double *)args[9], &sz[9], (double *)args[10], &sz[10],
3797                     (double *)args[11], &sz[11], (double *)args[12], &sz[12],
3798                     (double *)args[13], &sz[13], (double *)args[14], &sz[14],
3799                     (double *)args[15], &sz[15], (double *)args[16], &sz[16],
3800                     (double *)args[17], &sz[17]);
3801
3802             /* adjust values of output register */
3803             for (in = 0; in < block->nevout; ++in)
3804             {
3805                 block->evout[in] = block->evout[in] - *t;
3806             }
3807
3808             break;
3809         }
3810
3811         /*******************/
3812         /* function type 2 */
3813         /*******************/
3814         case 2 :
3815         {
3816             /* This is for compatibility */
3817             /* jroot is returned in g for old type */
3818             if (block->nevprt < 0)
3819             {
3820                 for (j = 0; j < block->ng; ++j)
3821                 {
3822                     block->g[j] = (double)block->jroot[j];
3823                 }
3824             }
3825
3826             /* no zero crossing */
3827             if (block->ztyp == 0)
3828             {
3829                 loc2 = (ScicosF2) loc;
3830                 (*loc2)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3831                         block->z, &block->nz,
3832                         block->evout, &block->nevout, block->rpar, &block->nrpar,
3833                         block->ipar, &block->nipar, (double **)block->inptr,
3834                         block->insz, &block->nin,
3835                         (double **)block->outptr, block->outsz, &block->nout,
3836                         block->scsptr);
3837             }
3838             /* with zero crossing */
3839             else
3840             {
3841                 loc2z = (ScicosF2z) loc;
3842                 (*loc2z)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3843                          block->z, &block->nz,
3844                          block->evout, &block->nevout, block->rpar, &block->nrpar,
3845                          block->ipar, &block->nipar, (double **)block->inptr,
3846                          block->insz, &block->nin,
3847                          (double **)block->outptr, block->outsz, &block->nout,
3848                          block->scsptr,
3849                          block->g, &block->ng);
3850             }
3851
3852             /* adjust values of output register */
3853             for (in = 0; in < block->nevout; ++in)
3854             {
3855                 block->evout[in] = block->evout[in] - *t;
3856             }
3857
3858             break;
3859         }
3860
3861         /*******************/
3862         /* function type 4 */
3863         /*******************/
3864         case 4 :
3865         {
3866             /* get pointer of the function type 4*/
3867             loc4 = (ScicosF4) loc;
3868
3869             (*loc4)(block, *flag);
3870
3871             break;
3872         }
3873
3874         /***********************/
3875         /* function type 10001 */
3876         /***********************/
3877         case 10001 :
3878         {
3879             /* This is for compatibility */
3880             /* jroot is returned in g for old type */
3881             if (block->nevprt < 0)
3882             {
3883                 for (j = 0; j < block->ng; ++j)
3884                 {
3885                     block->g[j] = (double)block->jroot[j];
3886                 }
3887             }
3888
3889             /* implicit block one entry for each input or output */
3890             for (in = 0 ; in < block->nin ; in++)
3891             {
3892                 args[in] = block->inptr[in];
3893                 sz[in] = block->insz[in];
3894             }
3895             for (out = 0; out < block->nout; out++)
3896             {
3897                 args[in + out] = block->outptr[out];
3898                 sz[in + out] = block->outsz[out];
3899             }
3900             /* with zero crossing */
3901             if (block->ztyp > 0)
3902             {
3903                 args[block->nin + block->nout] = block->g;
3904                 sz[block->nin + block->nout] = block->ng;
3905             }
3906
3907             loci1 = (ScicosFi) loc;
3908             (*loci1)(flag, &block->nevprt, t, block->res, block->xd, block->x,
3909                      &block->nx, block->z, &block->nz,
3910                      block->evout, &block->nevout, block->rpar, &block->nrpar,
3911                      block->ipar, &block->nipar,
3912                      (double *)args[0], &sz[0],
3913                      (double *)args[1], &sz[1], (double *)args[2], &sz[2],
3914                      (double *)args[3], &sz[3], (double *)args[4], &sz[4],
3915                      (double *)args[5], &sz[5], (double *)args[6], &sz[6],
3916                      (double *)args[7], &sz[7], (double *)args[8], &sz[8],
3917                      (double *)args[9], &sz[9], (double *)args[10], &sz[10],
3918                      (double *)args[11], &sz[11], (double *)args[12], &sz[12],
3919                      (double *)args[13], &sz[13], (double *)args[14], &sz[14],
3920                      (double *)args[15], &sz[15], (double *)args[16], &sz[16],
3921                      (double *)args[17], &sz[17]);
3922
3923             /* adjust values of output register */
3924             for (in = 0; in < block->nevout; ++in)
3925             {
3926                 block->evout[in] = block->evout[in] - *t;
3927             }
3928
3929             break;
3930         }
3931
3932         /***********************/
3933         /* function type 10002 */
3934         /***********************/
3935         case 10002 :
3936         {
3937             /* This is for compatibility */
3938             /* jroot is returned in g for old type */
3939             if (block->nevprt < 0)
3940             {
3941                 for (j = 0; j < block->ng; ++j)
3942                 {
3943                     block->g[j] = (double)block->jroot[j];
3944                 }
3945             }
3946
3947             /* implicit block, inputs and outputs given by a table of pointers */
3948             /* no zero crossing */
3949             if (block->ztyp == 0)
3950             {
3951                 loci2 = (ScicosFi2) loc;
3952                 (*loci2)(flag, &block->nevprt, t, block->res,
3953                          block->xd, block->x, &block->nx,
3954                          block->z, &block->nz,
3955                          block->evout, &block->nevout, block->rpar, &block->nrpar,
3956                          block->ipar, &block->nipar, (double **)block->inptr,
3957                          block->insz, &block->nin,
3958                          (double **)block->outptr, block->outsz, &block->nout);
3959             }
3960             /* with zero crossing */
3961             else
3962             {
3963                 loci2z = (ScicosFi2z) loc;
3964                 (*loci2z)(flag, &block->nevprt, t, block->res,
3965                           block->xd, block->x, &block->nx,
3966                           block->z, &block->nz,
3967                           block->evout, &block->nevout, block->rpar, &block->nrpar,
3968                           block->ipar, &block->nipar,
3969                           (double **)block->inptr, block->insz, &block->nin,
3970                           (double **)block->outptr, block->outsz, &block->nout,
3971                           block->g, &block->ng);
3972             }
3973
3974             /* adjust values of output register */
3975             for (in = 0; in < block->nevout; ++in)
3976             {
3977                 block->evout[in] = block->evout[in] - *t;
3978             }
3979
3980             break;
3981         }
3982
3983         /***********************/
3984         /* function type 10004 */
3985         /***********************/
3986         case 10004 :
3987         {
3988             /* get pointer of the function type 4*/
3989             loc4 = (ScicosF4) loc;
3990
3991             (*loc4)(block, *flag);
3992
3993             break;
3994         }
3995
3996         /***********/
3997         /* default */
3998         /***********/
3999         default :
4000         {
4001             sciprint(_("Undefined Function type\n"));
4002             *flag = -1000;
4003             return; /* exit */
4004         }
4005     }
4006     // sciprint("callf end  flag=%d\n",*flag);
4007     /* Implicit Solver & explicit block & flag==0 */
4008     /* adjust continuous state vector after call */
4009     if ((solver == IDA_BDF_Newton || solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes) && block->type < 10000 && *flag == 0)
4010     {
4011         block->xd  = ptr_d;
4012         if (flagi != 7)
4013         {
4014             for (k = 0; k < block->nx; k++)
4015             {
4016                 block->res[k] = block->res[k] - block->xd[k];
4017             }
4018         }
4019         else
4020         {
4021             for (k = 0; k < block->nx; k++)
4022             {
4023                 block->xd[k] = block->res[k];
4024             }
4025         }
4026     }
4027
4028     /* debug block */
4029     if (cosd > 1)
4030     {
4031         if (debug_block > -1)
4032         {
4033             if (*flag < 0)
4034             {
4035                 return;    /* error in block */
4036             }
4037             if (cosd != 3)
4038             {
4039                 sciprint(_("Leaving block %d \n"), C2F(curblk).kfun);
4040             }
4041             call_debug_scicos(block, flag, flagi, debug_block);
4042             /*call_debug_scicos(flag,kf,flagi,debug_block);*/
4043         }
4044     }
4045 } /* callf */
4046 /*--------------------------------------------------------------------------*/
4047 /* call_debug_scicos */
4048 static void call_debug_scicos(scicos_block *block, scicos_flag *flag, int flagi, int deb_blk)
4049 {
4050     voidf loc ;
4051     int solver = C2F(cmsolver).solver, k = 0;
4052     ScicosF4 loc4;
4053     double *ptr_d = NULL;
4054
4055     C2F(cosdebugcounter).counter = C2F(cosdebugcounter).counter + 1;
4056     C2F(scsptr).ptr = block->scsptr;
4057
4058     loc  = Blocks[deb_blk].funpt; /* GLOBAL */
4059     loc4 = (ScicosF4) loc;
4060
4061     /* continuous state */
4062     if ((solver == IDA_BDF_Newton || solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes) && block->type < 10000 && *flag == 0)
4063     {
4064         ptr_d = block->xd;
4065         block->xd  = block->res;
4066     }
4067
4068     // Temporarily replacing the block's computational function with DEBUG_BLOCK's so that sciblk4 will call %debug_scicos()
4069     block->scsptr = Blocks[deb_blk].scsptr;
4070     (*loc4)(block, *flag);
4071     block->scsptr = C2F(scsptr).ptr;
4072
4073     /* Implicit Solver & explicit block & flag==0 */
4074     /* adjust continuous state vector after call */
4075     if ((solver == IDA_BDF_Newton || solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes) && block->type < 10000 && *flag == 0)
4076     {
4077         block->xd  = ptr_d;
4078         if (flagi != 7)
4079         {
4080             for (k = 0; k < block->nx; k++)
4081             {
4082                 block->res[k] = block->res[k] - block->xd[k];
4083             }
4084         }
4085         else
4086         {
4087             for (k = 0; k < block->nx; k++)
4088             {
4089                 block->xd[k] = block->res[k];
4090             }
4091         }
4092     }
4093
4094     if (*flag < 0)
4095     {
4096         sciprint(_("Error in the Debug block \n"));
4097     }
4098 } /* call_debug_scicos */
4099 /*--------------------------------------------------------------------------*/
4100 /* simblk */
4101 static int simblk(realtype t, N_Vector yy, N_Vector yp, void *f_data)
4102 {
4103     double tx = 0., *x = NULL, *xd = NULL;
4104     int i = 0, nantest = 0;
4105
4106     tx = (double) t;
4107     x =  (double *) NV_DATA_S(yy);
4108     xd = (double *) NV_DATA_S(yp);
4109
4110     for (i = 0; i < *neq; i++)
4111     {
4112         xd[i] = 0;    /* à la place de "C2F(dset)(neq, &c_b14,xcdot , &c__1);"*/
4113     }
4114     C2F(ierode).iero = 0;
4115     *ierr = 0;
4116     odoit(&tx, x, xd, xd);
4117     C2F(ierode).iero = *ierr;
4118
4119     if (*ierr == 0)
4120     {
4121         nantest = 0;
4122         for (i = 0; i < *neq; i++) /* NaN checking */
4123         {
4124             if ((xd[i] - xd[i] != 0))
4125             {
4126                 Sciwarning(_("\nWarning: The computing function #%d returns a NaN/Inf"), i);
4127                 nantest = 1;
4128                 break;
4129             }
4130         }
4131         if (nantest == 1)
4132         {
4133             return 349;    /* recoverable error; */
4134         }
4135     }
4136
4137     return (abs(*ierr)); /* ierr>0 recoverable error; ierr>0 unrecoverable error; ierr=0: ok*/
4138
4139 } /* simblk */
4140 /*--------------------------------------------------------------------------*/
4141 /* grblk */
4142 static int grblk(realtype t, N_Vector yy, realtype *gout, void *g_data)
4143 {
4144     double tx = 0., *x = NULL;
4145     int jj = 0, nantest = 0;
4146
4147     tx = (double) t;
4148     x =  (double *) NV_DATA_S(yy);
4149
4150     C2F(ierode).iero = 0;
4151     *ierr = 0;
4152
4153     zdoit(&tx, x, x, (double*) gout);
4154
4155     if (*ierr == 0)
4156     {
4157         nantest = 0;
4158         for (jj = 0; jj < ng; jj++)
4159             if (gout[jj] - gout[jj] != 0)
4160             {
4161                 Sciwarning(_("\nWarning: The zero_crossing function #%d returns a NaN/Inf"), jj);
4162                 nantest = 1;
4163                 break;
4164             } /* NaN checking */
4165         if (nantest == 1)
4166         {
4167             return 350;    /* recoverable error; */
4168         }
4169     }
4170     C2F(ierode).iero = *ierr;
4171
4172     return 0;
4173 } /* grblk */
4174 /*--------------------------------------------------------------------------*/
4175 /* simblklsodar */
4176 static void simblklsodar(int * nequations, realtype * tOld, realtype * actual, realtype * res)
4177 {
4178     double tx = 0.;
4179     int i = 0;
4180
4181     tx = (double) * tOld;
4182
4183     for (i = 0; i < *nequations; ++i)
4184     {
4185         res[i] = 0;    /* à la place de "C2F(dset)(neq, &c_b14,xcdot , &c__1);"*/
4186     }
4187     C2F(ierode).iero = 0;
4188     *ierr = 0;
4189     odoit(&tx, actual, res, res);
4190     C2F(ierode).iero = *ierr;
4191
4192     if (*ierr == 0)
4193     {
4194         for (i = 0; i < *nequations; i++) /* NaN checking */
4195         {
4196             if ((res[i] - res[i] != 0))
4197             {
4198                 Sciwarning(_("\nWarning: The computing function #%d returns a NaN/Inf"), i);
4199             }
4200         }
4201     }
4202 } /* simblklsodar */
4203 /*--------------------------------------------------------------------------*/
4204 /* grblklsodar */
4205 static void grblklsodar(int * nequations, realtype * tOld, realtype * actual, int * ngc, realtype * res)
4206 {
4207     double tx = 0.;
4208     int jj = 0;
4209
4210     tx = (double) * tOld;
4211
4212     C2F(ierode).iero = 0;
4213     *ierr = 0;
4214
4215     zdoit(&tx, actual, actual, res);
4216
4217     if (*ierr == 0)
4218     {
4219         for (jj = 0; jj < *ngc; jj++)
4220         {
4221             if (res[jj] - res[jj] != 0)
4222             {
4223                 Sciwarning(_("\nWarning: The zero_crossing function #%d returns a NaN/Inf"), jj);
4224             } /* NaN checking */
4225         }
4226     }
4227 } /* grblklsodar */
4228 /*--------------------------------------------------------------------------*/
4229 /* simblkdaskr */
4230 static int simblkdaskr(realtype tres, N_Vector yy, N_Vector yp, N_Vector resval, void *rdata)
4231 {
4232     double tx = 0.;
4233     double *xc = NULL, *xcdot = NULL, *residual = NULL;
4234     realtype alpha = 0.;
4235
4236     UserData data;
4237
4238     realtype hh = 0.;
4239     int qlast = 0;
4240     int jj = 0, flag = 0, nantest = 0;
4241
4242     data = (UserData) rdata;
4243
4244     if (get_phase_simulation() == 1)
4245     {
4246         /* Just to update mode in a very special case, i.e., when initialization using modes fails.
4247         in this case, we relax all modes and try again one more time.
4248         */
4249         zdoit(&tx, NV_DATA_S(yy), NV_DATA_S(yp), NULL);
4250     }
4251
4252     hh = ZERO;
4253     flag = IDAGetCurrentStep(data->dae_mem, &hh);
4254     if (flag < 0)
4255     {
4256         *ierr = 200 + (-flag);
4257         return (*ierr);
4258     };
4259
4260     qlast = 0;
4261     flag = IDAGetCurrentOrder(data->dae_mem, &qlast);
4262     if (flag < 0)
4263     {
4264         *ierr = 200 + (-flag);
4265         return (*ierr);
4266     };
4267
4268     alpha = ZERO;
4269     for (jj = 0; jj < qlast; jj++)
4270     {
4271         alpha = alpha - ONE / (jj + 1);
4272     }
4273     if (hh != 0)
4274         // CJ=-alpha/hh;  // For function Get_Jacobian_cj
4275     {
4276         CJJ = -alpha / hh;
4277     }
4278     else
4279     {
4280         *ierr = 217;
4281         return (*ierr);
4282     }
4283     xc = (double *)  NV_DATA_S(yy);
4284     xcdot = (double *) NV_DATA_S(yp);
4285     residual = (double *) NV_DATA_S(resval);
4286     tx = (double) tres;
4287
4288     C2F(dcopy)(neq, xcdot, &c__1, residual, &c__1);
4289     *ierr = 0;
4290     C2F(ierode).iero = 0;
4291     odoit(&tx, xc, xcdot, residual);
4292
4293     C2F(ierode).iero = *ierr;
4294
4295     if (*ierr == 0)
4296     {
4297         nantest = 0;
4298         for (jj = 0; jj < *neq; jj++)
4299             if (residual[jj] - residual[jj] != 0) /* NaN checking */
4300             {
4301                 //sciprint(_("\nWarning: The residual function #%d returns a NaN"),jj);
4302                 nantest = 1;
4303                 break;
4304             }
4305         if (nantest == 1)
4306         {
4307             return 257;    /* recoverable error; */
4308         }
4309     }
4310
4311     return (abs(*ierr)); /* ierr>0 recoverable error; ierr>0 unrecoverable error; ierr=0: ok*/
4312 }/* simblkdaskr */
4313 /*--------------------------------------------------------------------------*/
4314 /* grblkdaskr */
4315 static int grblkdaskr(realtype t, N_Vector yy, N_Vector yp, realtype *gout, void *g_data)
4316 {
4317     double tx = 0.;
4318     int jj = 0, nantest = 0;
4319
4320     tx = (double) t;
4321
4322     *ierr = 0;
4323     C2F(ierode).iero = 0;
4324     zdoit(&tx, NV_DATA_S(yy), NV_DATA_S(yp), (double *) gout);
4325     if (*ierr == 0)
4326     {
4327         nantest = 0; /* NaN checking */
4328         for (jj = 0; jj < ng; jj++)
4329         {
4330             if (gout[jj] - gout[jj] != 0)
4331             {
4332                 Sciwarning(_("\nWarning: The zero-crossing function #%d returns a NaN"), jj);
4333                 nantest = 1;
4334                 break;
4335             }
4336         }
4337         if (nantest == 1)
4338         {
4339             return 258; /* recoverable error; */
4340         }
4341     }
4342     C2F(ierode).iero = *ierr;
4343     return (*ierr);
4344 }/* grblkdaskr */
4345 /*--------------------------------------------------------------------------*/
4346 /* simblkddaskr */
4347 static void simblkddaskr(realtype *tOld, realtype *actual, realtype *actualP, realtype *res, int *flag, double *dummy1, int *dummy2)
4348 {
4349     double tx = 0.;
4350
4351     int jj = 0;
4352
4353     if (get_phase_simulation() == 1)
4354     {
4355         /* Just to update mode in a very special case, i.e., when initialization using modes fails.
4356         in this case, we relax all modes and try again one more time.
4357         */
4358         zdoit(&tx, actual, actualP, NULL);
4359     }
4360
4361     CJJ = 6;
4362
4363     tx = (double) * tOld;
4364     *flag = 0;
4365
4366     C2F(dcopy)(neq, actualP, &c__1, res, &c__1);
4367     *ierr = 0;
4368     C2F(ierode).iero = 0;
4369     odoit(&tx, actual, actualP, res);
4370     C2F(ierode).iero = *ierr;
4371
4372     if (*ierr == 0)
4373     {
4374         for (jj = 0; jj < *neq; jj++)
4375             if (res[jj] - res[jj] != 0) /* NaN checking */
4376             {
4377                 Sciwarning(_("\nWarning: The residual function #%d returns a NaN"), jj);
4378                 *flag = -1; /* recoverable error; */
4379                 return;
4380             }
4381     }
4382     else
4383     {
4384         *flag = -2;
4385         return;
4386     }
4387     /* *flag=-1 recoverable error; *flag=-2 unrecoverable error; *flag=0: ok*/
4388
4389 }/* simblkddaskr */
4390 /*--------------------------------------------------------------------------*/
4391 /* grblkddaskr */
4392 static void grblkddaskr(int *nequations, realtype *tOld, realtype *actual, int *ngc, realtype *res, double *dummy1, int *dummy2)
4393 {
4394     double tx = 0.;
4395     int jj = 0;
4396
4397     tx = (double) * tOld;
4398
4399     *ierr = 0;
4400     C2F(ierode).iero = 0;
4401     zdoit(&tx, actual, actual, res);
4402     C2F(ierode).iero = *ierr;
4403
4404     if (*ierr == 0)
4405     {
4406         /* NaN checking */
4407         for (jj = 0; jj < *ngc; jj++)
4408         {
4409             if (res[jj] - res[jj] != 0)
4410             {
4411                 Sciwarning(_("\nWarning: The zero-crossing function #%d returns a NaN"), jj);
4412                 return;
4413             }
4414         }
4415     }
4416     else
4417     {
4418         sciprint(_("\nError: Problem in the evaluation of a root function"));
4419         return;
4420     }
4421
4422 }/* grblkddaskr */
4423 /*--------------------------------------------------------------------------*/
4424 /* jacpsol */
4425 static void jacpsol(realtype *res, int *ires, int *neq, realtype *tOld, realtype *actual, realtype *actualP,
4426                     realtype *rewt, realtype *savr, realtype *wk, realtype *h, realtype *cj, realtype *wp, int *iwp,
4427                     int *ier, double *dummy1, int *dummy2)
4428 {
4429     /* Here, we compute the system preconditioner matrix P, which is actually the jacobian matrix,
4430        so P(i,j) = dres(i)/dactual(j) + cj*dres(i)/dactualP(j), and we LU-decompose it. */
4431     int i = 0, j = 0, nrow = 0, info = 0;
4432     realtype tx = 0, del = 0, delinv = 0, ysave = 0, ypsave = 0;
4433     realtype * e = NULL;
4434
4435     tx = *tOld;
4436
4437     /* Work array used to evaluate res(*tOld, actual + small_increment, actualP + small_increment).
4438        savr already contains res(*tOld, actual, actualP). */
4439     e = (realtype *) calloc(*neq, sizeof(realtype));
4440
4441     for (i = 0; i < *neq; ++i)
4442     {
4443         del =  max (SQuround * max(fabs(actual[i]), fabs(*h * actualP[i])), 1. / rewt[i]);
4444         del *= (*h * actualP[i] >= 0) ? 1 : -1;
4445         del =  (actual[i] + del) - actual[i];
4446         ysave  = actual[i];
4447         ypsave = actualP[i];
4448         actual[i]  += del;
4449         actualP[i] += *cj * del;
4450         *ierr = 0;
4451         C2F(ierode).iero = 0;
4452         simblkddaskr(tOld, actual, actualP, e, ires, dummy1, dummy2);
4453         C2F(ierode).iero = *ierr;
4454         if (*ires == 0)
4455         {
4456             delinv = 1. / del;
4457             for (j = 0; j < *neq; ++j)
4458             {
4459                 wp[nrow + j] = (e[j] - savr[j]) * delinv;
4460                 /* NaN test */
4461                 if (wp[nrow + j] - wp[nrow + j] != 0)
4462                 {
4463                     Sciwarning(_("\nWarning: The preconditioner evaluation function returns a NaN at index #%d."), nrow + j);
4464                     *ier = -1;
4465                 }
4466             }
4467             nrow       += *neq;
4468             actual[i]  =  ysave;
4469             actualP[i] =  ypsave;
4470         }
4471         else
4472         {
4473             sciprint(_("\nError: The preconditioner evaluation function failed."));
4474             *ier = -1;
4475             free(e);
4476             return;
4477         }
4478     }
4479
4480     /* Proceed to LU factorization of P. */
4481     C2F(dgefa) (wp, neq, neq, iwp, &info);
4482     if (info != 0)
4483     {
4484         sciprint(_("\nError: Failed to factor the preconditioner."));
4485         *ier = -1;
4486     }
4487
4488     free(e);
4489 }/* jacpsol */
4490 /*--------------------------------------------------------------------------*/
4491 /* psol */
4492 static void psol(int *neq, realtype *tOld, realtype *actual, realtype *actualP,
4493                  realtype *savr, realtype *wk, realtype *cj, realtype *wght, realtype *wp,
4494                  int *iwp, realtype *b, realtype *eplin, int *ier, double *dummy1, int *dummy2)
4495 {
4496     /* This function "applies" the inverse of the preconditioner to 'b' (computes P^-1*b).
4497        It is done by solving P*x = b using the linpack routine 'dgesl'. */
4498     int i = 0, job = 0;
4499
4500     C2F(dgesl) (wp, neq, neq, iwp, b, &job);
4501
4502     /* NaN test */
4503     for (i = 0; i < *neq; ++i)
4504     {
4505         if (b[i] - b[i] != 0)
4506         {
4507             Sciwarning(_("\nWarning: The preconditioner application function returns a NaN at index #%d."), i);
4508             /* Indicate a recoverable error, meaning that the step will be retried with the same step size
4509                but with a call to 'jacpsol' to update necessary data, unless the Jacobian data is current,
4510                in which case the step will be retried with a smaller step size. */
4511             *ier = 1;
4512         }
4513     }
4514 }/* psol */
4515 /*--------------------------------------------------------------------------*/
4516 /* Subroutine addevs */
4517 static void addevs(double t, int *evtnb, int *ierr1)
4518 {
4519     static int i = 0, j = 0;
4520
4521     /* Function Body */
4522     *ierr1 = 0;
4523     if (evtspt[*evtnb] != -1)
4524     {
4525         if ((evtspt[*evtnb] == 0) && (*pointi == *evtnb))
4526         {
4527             tevts[*evtnb] = t;
4528             return;
4529         }
4530         else
4531         {
4532             if (*pointi == *evtnb)
4533             {
4534                 *pointi = evtspt[*evtnb]; /* remove from chain */
4535             }
4536             else
4537             {
4538                 i = *pointi;
4539                 while (*evtnb != evtspt[i])
4540                 {
4541                     i = evtspt[i];
4542                 }
4543                 evtspt[i] = evtspt[*evtnb]; /* remove old evtnb from chain */
4544                 if (TCritWarning == 0)
4545                 {
4546                     sciprint(_("\n Warning: an event is reprogrammed at t=%g by removing another"), t );
4547                     sciprint(_("\n         (already programmed) event. There may be an error in"));
4548                     Sciwarning(_("\n         your model. Please check your model\n"));
4549                     TCritWarning = 1;
4550                 }
4551                 do_cold_restart(); /* the erased event could be a critical
4552                                                                    event, so do_cold_restart is added to
4553                                                                    refresh the critical event table */
4554             }
4555             evtspt[*evtnb] = 0;
4556             tevts[*evtnb] = t;
4557         }
4558     }
4559     else
4560     {
4561         evtspt[*evtnb] = 0;
4562         tevts[*evtnb] = t;
4563     }
4564     if (*pointi == 0)
4565     {
4566         *pointi = *evtnb;
4567         return;
4568     }
4569     if (t < tevts[*pointi])
4570     {
4571         evtspt[*evtnb] = *pointi;
4572         *pointi = *evtnb;
4573         return;
4574     }
4575     i = *pointi;
4576
4577 L100:
4578     if (evtspt[i] == 0)
4579     {
4580         evtspt[i] = *evtnb;
4581         return;
4582     }
4583     if (t >= tevts[evtspt[i]])
4584     {
4585         j = evtspt[i];
4586         if (evtspt[j] == 0)
4587         {
4588             evtspt[j] = *evtnb;
4589             return;
4590         }
4591         i = j;
4592         goto L100;
4593     }
4594     else
4595     {
4596         evtspt[*evtnb] = evtspt[i];
4597         evtspt[i] = *evtnb;
4598     }
4599 } /* addevs */
4600 /*--------------------------------------------------------------------------*/
4601 /* Subroutine putevs */
4602 void putevs(double *t, int *evtnb, int *ierr1)
4603 {
4604     /* Function Body */
4605     *ierr1 = 0;
4606     if (evtspt[*evtnb] != -1)
4607     {
4608         *ierr1 = 1;
4609         return;
4610     }
4611     else
4612     {
4613         evtspt[*evtnb] = 0;
4614         tevts[*evtnb] = *t;
4615     }
4616     if (*pointi == 0)
4617     {
4618         *pointi = *evtnb;
4619         return;
4620     }
4621     evtspt[*evtnb] = *pointi;
4622     *pointi = *evtnb;
4623 } /* putevs */
4624 /*--------------------------------------------------------------------------*/
4625 /* Subroutine idoit */
4626 static void idoit(double *told)
4627 {
4628     /* initialisation (propagation of constant blocks outputs) */
4629     /*     Copyright INRIA */
4630
4631     int i2 = 0;
4632     scicos_flag flag = 0;
4633     int i = 0, j = 0;
4634     int ierr1 = 0;
4635
4636     ScicosImport *scs_imp = NULL;
4637     int *kf = NULL;
4638
4639     scs_imp = getscicosimportptr();
4640
4641     flag = 1;
4642     for (j = 0; j < * (scs_imp->niord); j++)
4643     {
4644         kf = &scs_imp->iord[j];
4645         C2F(curblk).kfun = *kf; /* */
4646         if (scs_imp->funtyp[*kf - 1] > -1)
4647         {
4648             /* continuous state */
4649             if (scs_imp->blocks[*kf - 1].nx > 0)
4650             {
4651                 scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4652                 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4653             }
4654             scs_imp->blocks[*kf - 1].nevprt = scs_imp->iord[j + * (scs_imp->niord)];
4655             /*callf(told, xd, x, x,x,&flag);*/
4656             callf(told, &scs_imp->blocks[*kf - 1], &flag);
4657             if (flag < 0)
4658             {
4659                 *ierr = 5 - flag;
4660                 return;
4661             }
4662         }
4663         if (scs_imp->blocks[*kf - 1].nevout > 0)
4664         {
4665             if (scs_imp->funtyp[*kf - 1] < 0)
4666             {
4667                 i = synchro_nev(scs_imp, *kf, ierr);
4668                 if (*ierr != 0)
4669                 {
4670                     return;
4671                 }
4672                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4673                 putevs(told, &i2, &ierr1);
4674                 if (ierr1 != 0)
4675                 {
4676                     /* event conflict */
4677                     *ierr = 3;
4678                     return;
4679                 }
4680                 doit(told);
4681                 if (*ierr != 0)
4682                 {
4683                     return;
4684                 }
4685             }
4686         }
4687     }
4688 } /* idoit_ */
4689 /*--------------------------------------------------------------------------*/
4690 /* Subroutine doit */
4691 static void doit(double *told)
4692 {
4693     /* propagation of blocks outputs on discrete activations */
4694     /*     Copyright INRIA */
4695
4696     int i = 0, i2 = 0;
4697     scicos_flag flag = 0;
4698     int nord = 0;
4699     int ierr1 = 0;
4700     int ii = 0, kever = 0;
4701
4702     ScicosImport *scs_imp = NULL;
4703     int *kf = NULL;
4704
4705     scs_imp = getscicosimportptr();
4706
4707     kever = *pointi;
4708     *pointi = evtspt[kever];
4709     evtspt[kever] = -1;
4710
4711     nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever - 1];
4712     if (nord == 0)
4713     {
4714         return;
4715     }
4716
4717     for (ii = scs_imp->ordptr[kever - 1]; ii <= scs_imp->ordptr[kever] - 1 ; ii++)
4718     {
4719         kf = &scs_imp->ordclk[ii - 1];
4720         C2F(curblk).kfun = *kf;
4721         if (scs_imp->funtyp[*kf - 1] > -1)
4722         {
4723             /* continuous state */
4724             if (scs_imp->blocks[*kf - 1].nx > 0)
4725             {
4726                 scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4727                 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4728             }
4729             scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
4730             flag = 1;
4731             /*callf(told, xd, x, x,x,&flag);*/
4732             callf(told, &scs_imp->blocks[*kf - 1], &flag);
4733             if (flag < 0)
4734             {
4735                 *ierr = 5 - flag;
4736                 return;
4737             }
4738         }
4739
4740         /* Initialize tvec */
4741         if (scs_imp->blocks[*kf - 1].nevout > 0)
4742         {
4743             if (scs_imp->funtyp[*kf - 1] < 0)
4744             {
4745                 i = synchro_nev(scs_imp, *kf, ierr);
4746                 if (*ierr != 0)
4747                 {
4748                     return;
4749                 }
4750                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4751                 putevs(told, &i2, &ierr1);
4752                 if (ierr1 != 0)
4753                 {
4754                     /* event conflict */
4755                     *ierr = 3;
4756                     return;
4757                 }
4758                 doit(told);
4759                 if (*ierr != 0)
4760                 {
4761                     return;
4762                 }
4763             }
4764         }
4765     }
4766 } /* doit_ */
4767 /*--------------------------------------------------------------------------*/
4768 /* Subroutine cdoit */
4769 static void cdoit(double *told)
4770 {
4771     /* propagation of continuous blocks outputs */
4772     /*     Copyright INRIA */
4773
4774     int i2 = 0;
4775     scicos_flag flag = 0;
4776     int ierr1 = 0;
4777     int i = 0, j = 0;
4778
4779     ScicosImport *scs_imp = NULL;
4780     int *kf = NULL;
4781
4782     scs_imp = getscicosimportptr();
4783
4784     /* Function Body */
4785     for (j = 0; j < * (scs_imp->ncord); j++)
4786     {
4787         kf = &scs_imp->cord[j];
4788         C2F(curblk).kfun = *kf;
4789         /* continuous state */
4790         if (scs_imp->blocks[*kf - 1].nx > 0)
4791         {
4792             scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4793             scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4794         }
4795         scs_imp->blocks[*kf - 1].nevprt = scs_imp->cord[j + * (scs_imp->ncord)];
4796         if (scs_imp->funtyp[*kf - 1] > -1)
4797         {
4798             flag = 1;
4799             /*callf(told, xd, x, x,x,&flag);*/
4800             callf(told, &scs_imp->blocks[*kf - 1], &flag);
4801             if (flag < 0)
4802             {
4803                 *ierr = 5 - flag;
4804                 return;
4805             }
4806         }
4807
4808         /* Initialize tvec */
4809         if (scs_imp->blocks[*kf - 1].nevout > 0)
4810         {
4811             if (scs_imp->funtyp[*kf - 1] < 0)
4812             {
4813                 i = synchro_nev(scs_imp, *kf, ierr);
4814                 if (*ierr != 0)
4815                 {
4816                     return;
4817                 }
4818                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4819                 putevs(told, &i2, &ierr1);
4820                 if (ierr1 != 0)
4821                 {
4822                     /* event conflict */
4823                     *ierr = 3;
4824                     return;
4825                 }
4826                 doit(told);
4827                 if (*ierr != 0)
4828                 {
4829                     return;
4830                 }
4831             }
4832         }
4833     }
4834 } /* cdoit_ */
4835 /*--------------------------------------------------------------------------*/
4836 /* Subroutine ddoit */
4837 static void ddoit(double *told)
4838 {
4839     /* update states & event out on discrete activations */
4840     /*     Copyright INRIA */
4841
4842     int i2 = 0, j = 0;
4843     scicos_flag flag = 0;
4844     int kiwa = 0;
4845     int i = 0, i3 = 0, ierr1 = 0;
4846     int ii = 0, keve = 0;
4847
4848     ScicosImport *scs_imp = NULL;
4849     int *kf = NULL;
4850
4851     scs_imp = getscicosimportptr();
4852
4853     /* Function Body */
4854     kiwa = 0;
4855     edoit(told, &kiwa);
4856     if (*ierr != 0)
4857     {
4858         return;
4859     }
4860
4861     /* update continuous and discrete states on event */
4862     if (kiwa == 0)
4863     {
4864         return;
4865     }
4866     for (i = 0; i < kiwa; i++)
4867     {
4868         keve = iwa[i];
4869         if (critev[keve] != 0)
4870         {
4871             hot = 0;
4872         }
4873         i2 = scs_imp->ordptr[keve] - 1;
4874         for (ii = scs_imp->ordptr[keve - 1]; ii <= i2; ii++)
4875         {
4876             kf = &scs_imp->ordclk[ii - 1];
4877             C2F(curblk).kfun = *kf;
4878             /* continuous state */
4879             if (scs_imp->blocks[*kf - 1].nx > 0)
4880             {
4881                 scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4882                 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4883             }
4884
4885             scs_imp->blocks[*kf - 1].nevprt = scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1];
4886
4887             if (scs_imp->blocks[*kf - 1].nevout > 0)
4888             {
4889                 if (scs_imp->funtyp[*kf - 1] >= 0)
4890                 {
4891                     /* initialize evout */
4892                     for (j = 0; j < scs_imp->blocks[*kf - 1].nevout; j++)
4893                     {
4894                         scs_imp->blocks[*kf - 1].evout[j] = -1;
4895                     }
4896                     flag = 3;
4897
4898                     if (scs_imp->blocks[*kf - 1].nevprt > 0) /* if event has continuous origin don't call*/
4899                     {
4900                         /*callf(told, xd, x, x ,x,&flag);*/
4901                         callf(told, &scs_imp->blocks[*kf - 1], &flag);
4902                         if (flag < 0)
4903                         {
4904                             *ierr = 5 - flag;
4905                             return;
4906                         }
4907                     }
4908
4909                     for (j = 0; j < scs_imp->blocks[*kf - 1].nevout; j++)
4910                     {
4911                         if (scs_imp->blocks[*kf - 1].evout[j] >= 0.)
4912                         {
4913                             i3 = j + scs_imp->clkptr[*kf - 1] ;
4914                             addevs(scs_imp->blocks[*kf - 1].evout[j] + (*told), &i3, &ierr1);
4915                             if (ierr1 != 0)
4916                             {
4917                                 /* event conflict */
4918                                 *ierr = 3;
4919                                 return;
4920                             }
4921                         }
4922                     }
4923                 }
4924             }
4925
4926             if (scs_imp->blocks[*kf - 1].nevprt > 0)
4927             {
4928                 if (scs_imp->blocks[*kf - 1].nx + scs_imp->blocks[*kf - 1].nz + scs_imp->blocks[*kf - 1].noz > 0 || \
4929                         *scs_imp->blocks[*kf - 1].work != NULL)
4930                 {
4931                     /*  if a hidden state exists, must also call (for new scope eg)  */
4932                     /*  to avoid calling non-real activations */
4933                     flag = 2;
4934                     /*callf(told, xd, x, x,x,&flag);*/
4935                     callf(told, &scs_imp->blocks[*kf - 1], &flag);
4936                     if (flag < 0)
4937                     {
4938                         *ierr = 5 - flag;
4939                         return;
4940                     }
4941                 }
4942             }
4943             else
4944             {
4945                 if (*scs_imp->blocks[*kf - 1].work != NULL)
4946                 {
4947                     flag = 2;
4948                     scs_imp->blocks[*kf - 1].nevprt = 0; /* in case some hidden continuous blocks need updating */
4949                     /*callf(told, xd, x, x,x,&flag);*/
4950                     callf(told, &scs_imp->blocks[*kf - 1], &flag);
4951                     if (flag < 0)
4952                     {
4953                         *ierr = 5 - flag;
4954                         return;
4955                     }
4956                 }
4957             }
4958         }
4959     }
4960 } /* ddoit_ */
4961 /*--------------------------------------------------------------------------*/
4962 /* Subroutine edoit */
4963 static void edoit(double *told, int *kiwa)
4964 {
4965     /* update blocks output on discrete activations */
4966     /*     Copyright INRIA */
4967
4968     int i2 = 0;
4969     scicos_flag flag = 0;
4970     int ierr1 = 0, i = 0;
4971     int kever = 0, ii = 0;
4972
4973     ScicosImport *scs_imp = NULL;
4974     int *kf = NULL;
4975     int nord = 0;
4976
4977     scs_imp = getscicosimportptr();
4978
4979     /* Function Body */
4980     kever = *pointi;
4981
4982     *pointi = evtspt[kever];
4983     evtspt[kever] = -1;
4984
4985     nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever - 1];
4986     if (nord == 0)
4987     {
4988         return;
4989     }
4990     iwa[*kiwa] = kever;
4991     ++(*kiwa);
4992     for (ii = scs_imp->ordptr[kever - 1]; ii <= scs_imp->ordptr[kever] - 1; ii++)
4993     {
4994         kf = &scs_imp->ordclk[ii - 1];
4995         C2F(curblk).kfun = *kf;
4996
4997         if (scs_imp->funtyp[*kf - 1] > -1)
4998         {
4999             /* continuous state */
5000             if (scs_imp->blocks[*kf - 1].nx > 0)
5001             {
5002                 scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
5003                 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5004             }
5005
5006             scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5007
5008             flag = 1;
5009             /*callf(told, xd, x, x,x,&flag);*/
5010             callf(told, &scs_imp->blocks[*kf - 1], &flag);
5011             if (flag < 0)
5012             {
5013                 *ierr = 5 - flag;
5014                 return;
5015             }
5016         }
5017
5018         /* Initialize tvec */
5019         if (scs_imp->blocks[*kf - 1].nevout > 0)
5020         {
5021             if (scs_imp->funtyp[*kf - 1] < 0)
5022             {
5023                 i = synchro_nev(scs_imp, *kf, ierr);
5024                 if (*ierr != 0)
5025                 {
5026                     return;
5027                 }
5028                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5029                 putevs(told, &i2, &ierr1);
5030                 if (ierr1 != 0)
5031                 {
5032                     /* event conflict */
5033                     *ierr = 3;
5034                     return;
5035                 }
5036                 edoit(told, kiwa);
5037                 if (*ierr != 0)
5038                 {
5039                     return;
5040                 }
5041             }
5042         }
5043     }
5044 } /* edoit_ */
5045 /*--------------------------------------------------------------------------*/
5046 /* Subroutine odoit */
5047 static void odoit(double *told, double *xt, double *xtd, double *residual)
5048 {
5049     /* update blocks derivative of continuous time block */
5050     /*     Copyright INRIA */
5051
5052     int i2 = 0;
5053     scicos_flag flag = 0;
5054     int keve = 0, kiwa = 0;
5055     int ierr1 = 0, i = 0;
5056     int ii = 0, jj = 0;
5057
5058     ScicosImport *scs_imp = NULL;
5059     int *kf = NULL;
5060
5061     scs_imp = getscicosimportptr();
5062
5063     /* Function Body */
5064     kiwa = 0;
5065     for (jj = 0; jj < * (scs_imp->noord); jj++)
5066     {
5067         kf = &scs_imp->oord[jj];
5068         C2F(curblk).kfun = *kf;
5069         /* continuous state */
5070         if (scs_imp->blocks[*kf - 1].nx > 0)
5071         {
5072             scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5073             scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5074             scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5075         }
5076
5077         scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[jj + * (scs_imp->noord)];
5078         if (scs_imp->funtyp[*kf - 1] > -1)
5079         {
5080             flag = 1;
5081             /*callf(told, xtd, xt, residual,x,&flag);*/
5082             callf(told, &scs_imp->blocks[*kf - 1], &flag);
5083             if (flag < 0)
5084             {
5085                 *ierr = 5 - flag;
5086                 return;
5087             }
5088         }
5089
5090         if (scs_imp->blocks[*kf - 1].nevout > 0)
5091         {
5092             if (scs_imp->funtyp[*kf - 1] < 0)
5093             {
5094                 if (scs_imp->blocks[*kf - 1].nmode > 0)
5095                 {
5096                     i2 = scs_imp->blocks[*kf - 1].mode[0] + scs_imp->clkptr[*kf - 1] - 1;
5097                 }
5098                 else
5099                 {
5100                     i = synchro_nev(scs_imp, *kf, ierr);
5101                     if (*ierr != 0)
5102                     {
5103                         return;
5104                     }
5105                     i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5106                 }
5107                 putevs(told, &i2, &ierr1);
5108                 if (ierr1 != 0)
5109                 {
5110                     /* event conflict */
5111                     *ierr = 3;
5112                     return;
5113                 }
5114                 ozdoit(told, xt, xtd, &kiwa);
5115                 if (*ierr != 0)
5116                 {
5117                     return;
5118                 }
5119             }
5120         }
5121     }
5122
5123     /*  update states derivatives */
5124     for (ii = 0; ii < * (scs_imp->noord); ii++)
5125     {
5126         kf = &scs_imp->oord[ii];
5127         C2F(curblk).kfun = *kf;
5128         if (scs_imp->blocks[*kf - 1].nx > 0 || \
5129                 *scs_imp->blocks[*kf - 1].work != NULL)
5130         {
5131             /* work tests if a hidden state exists, used for delay block */
5132             flag = 0;
5133             /* continuous state */
5134             if (scs_imp->blocks[*kf - 1].nx > 0)
5135             {
5136                 scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5137                 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5138                 scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5139             }
5140             scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[ii + * (scs_imp->noord)];
5141             /*callf(told, xtd, xt, residual,xt,&flag);*/
5142             callf(told, &scs_imp->blocks[*kf - 1], &flag);
5143
5144             if (flag < 0)
5145             {
5146                 *ierr = 5 - flag;
5147                 return;
5148             }
5149         }
5150     }
5151
5152     for (i = 0; i < kiwa; i++)
5153     {
5154         keve = iwa[i];
5155         for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1; ii++)
5156         {
5157             kf = &scs_imp->ordclk[ii - 1];
5158             C2F(curblk).kfun = *kf;
5159             if (scs_imp->blocks[*kf - 1].nx > 0 || \
5160                     *scs_imp->blocks[*kf - 1].work != NULL)
5161             {
5162                 /* work tests if a hidden state exists */
5163                 flag = 0;
5164                 /* continuous state */
5165                 if (scs_imp->blocks[*kf - 1].nx > 0)
5166                 {
5167                     scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5168                     scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5169                     scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5170                 }
5171                 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5172                 /*callf(told, xtd, xt, residual,xt,&flag);*/
5173                 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5174
5175                 if (flag < 0)
5176                 {
5177                     *ierr = 5 - flag;
5178                     return;
5179                 }
5180             }
5181         }
5182     }
5183 } /* odoit_ */
5184 /*--------------------------------------------------------------------------*/
5185 /* Subroutine reinitdoit */
5186 static void reinitdoit(double *told)
5187 {
5188     /* update blocks xproperties of continuous time block */
5189     /*     Copyright INRIA */
5190
5191     int i2 = 0;
5192     scicos_flag flag = 0;
5193     int keve = 0, kiwa = 0;
5194     int ierr1 = 0, i = 0;
5195     int ii = 0, jj = 0;
5196
5197     ScicosImport *scs_imp = NULL;
5198     int *kf = NULL;
5199
5200     scs_imp = getscicosimportptr();
5201
5202     /* Function Body */
5203     kiwa = 0;
5204     for (jj = 0; jj < * (scs_imp->noord); jj++)
5205     {
5206         kf = &scs_imp->oord[jj];
5207         C2F(curblk).kfun = *kf;
5208         /* continuous state */
5209         if (scs_imp->blocks[*kf - 1].nx > 0)
5210         {
5211             scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
5212             scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5213         }
5214         scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[jj + * (scs_imp->noord)];
5215         if (scs_imp->funtyp[*kf - 1] > -1)
5216         {
5217             flag = 1;
5218             /*callf(told, xd, x, x,x,&flag);*/
5219             callf(told, &scs_imp->blocks[*kf - 1], &flag);
5220             if (flag < 0)
5221             {
5222                 *ierr = 5 - flag;
5223                 return;
5224             }
5225         }
5226
5227         if (scs_imp->blocks[*kf - 1].nevout > 0 && scs_imp->funtyp[*kf - 1] < 0)
5228         {
5229             i = synchro_nev(scs_imp, *kf, ierr);
5230             if (*ierr != 0)
5231             {
5232                 return;
5233             }
5234             if (scs_imp->blocks[*kf - 1].nmode > 0)
5235             {
5236                 scs_imp->blocks[*kf - 1].mode[0] = i;
5237             }
5238             i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5239             putevs(told, &i2, &ierr1);
5240             if (ierr1 != 0)
5241             {
5242                 /* event conflict */
5243                 *ierr = 3;
5244                 return;
5245             }
5246             doit(told);
5247             if (*ierr != 0)
5248             {
5249                 return;
5250             }
5251         }
5252     }
5253
5254     /* reinitialize */
5255     for (ii = 0; ii < * (scs_imp->noord); ii++)
5256     {
5257         kf = &scs_imp->oord[ii];
5258         C2F(curblk).kfun = *kf;
5259         if (scs_imp->blocks[*kf - 1].nx > 0)
5260         {
5261             flag = 7;
5262             scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
5263             scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5264             scs_imp->blocks[*kf - 1].res = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5265             scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[ii + * (scs_imp->noord)];
5266             scs_imp->blocks[*kf - 1].xprop = &scs_imp->xprop[-1 + scs_imp->xptr[*kf - 1]];
5267             /*callf(told, xd, x, xd,x,&flag);*/
5268             callf(told, &scs_imp->blocks[*kf - 1], &flag);
5269             if (flag < 0)
5270             {
5271                 *ierr = 5 - flag;
5272                 return;
5273             }
5274         }
5275     }
5276
5277     for (i = 0; i < kiwa; i++)
5278     {
5279         keve = iwa[i];
5280         for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1; ii++)
5281         {
5282             kf = &scs_imp->ordclk[ii - 1];
5283             C2F(curblk).kfun = *kf;
5284             if (scs_imp->blocks[*kf - 1].nx > 0)
5285             {
5286                 flag = 7;
5287                 scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
5288                 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5289                 scs_imp->blocks[*kf - 1].res = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5290                 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5291                 scs_imp->blocks[*kf - 1].xprop = &scs_imp->xprop[-1 + scs_imp->xptr[*kf - 1]];
5292                 /*callf(told, xd, x, xd,x,&flag);*/
5293                 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5294                 if (flag < 0)
5295                 {
5296                     *ierr = 5 - flag;
5297                     return;
5298                 }
5299             }
5300         }
5301     }
5302 } /* reinitdoit_ */
5303 /*--------------------------------------------------------------------------*/
5304 /* Subroutine ozdoit */
5305 static void ozdoit(double *told, double *xt, double *xtd, int *kiwa)
5306 {
5307     /* update blocks output of continuous time block on discrete activations */
5308     /*     Copyright INRIA */
5309
5310     int i2 = 0;
5311     scicos_flag flag = 0;
5312     int nord = 0;
5313     int ierr1 = 0, i = 0;
5314     int ii = 0, kever = 0;
5315
5316     ScicosImport *scs_imp = NULL;
5317     int *kf = NULL;
5318
5319     scs_imp = getscicosimportptr();
5320
5321     /* Function Body */
5322     kever = *pointi;
5323     *pointi = evtspt[kever];
5324     evtspt[kever] = -1;
5325
5326     nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever - 1];
5327     if (nord == 0)
5328     {
5329         return;
5330     }
5331     iwa[*kiwa] = kever;
5332     ++(*kiwa);
5333
5334     for (ii = scs_imp->ordptr[kever - 1]; ii <= scs_imp->ordptr[kever] - 1; ii++)
5335     {
5336         kf = &scs_imp->ordclk[ii - 1];
5337         C2F(curblk).kfun = *kf;
5338         if (scs_imp->funtyp[*kf - 1] > -1)
5339         {
5340             /* continuous state */
5341             if (scs_imp->blocks[*kf - 1].nx > 0)
5342             {
5343                 scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5344                 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5345             }
5346             scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5347             flag = 1;
5348             /*callf(told, xtd, xt, xt,x,&flag);*/
5349             callf(told, &scs_imp->blocks[*kf - 1], &flag);
5350             if (flag < 0)
5351             {
5352                 *ierr = 5 - flag;
5353                 return;
5354             }
5355         }
5356
5357         /* Initialize tvec */
5358         if (scs_imp->blocks[*kf - 1].nevout > 0)
5359         {
5360             if (scs_imp->funtyp[*kf - 1] < 0)
5361             {
5362                 if (phase == 1 || scs_imp->blocks[*kf - 1].nmode == 0)
5363                 {
5364                     i = synchro_nev(scs_imp, *kf, ierr);
5365                     if (*ierr != 0)
5366                     {
5367                         return;
5368                     }
5369                 }
5370                 else
5371                 {
5372                     i = scs_imp->blocks[*kf - 1].mode[0];
5373                 }
5374                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5375                 putevs(told, &i2, &ierr1);
5376                 if (ierr1 != 0)
5377                 {
5378                     /* event conflict */
5379                     *ierr = 3;
5380                     return;
5381                 }
5382                 ozdoit(told, xt, xtd, kiwa);
5383             }
5384         }
5385     }
5386 } /* ozdoit_ */
5387 /*--------------------------------------------------------------------------*/
5388 /* Subroutine zdoit */
5389 static void zdoit(double *told, double *xt, double *xtd, double *g)
5390 {
5391     /* update blocks zcross of continuous time block  */
5392     /*     Copyright INRIA */
5393     int i2 = 0;
5394     scicos_flag flag = 0;
5395     int keve = 0, kiwa = 0;
5396     int ierr1 = 0, i = 0, j = 0;
5397     int ii = 0, jj = 0;
5398
5399     ScicosImport *scs_imp = NULL;
5400     int *kf = NULL;
5401
5402     scs_imp = getscicosimportptr();
5403
5404     /* Function Body */
5405     for (i = 0; i < * (scs_imp->ng); i++)
5406     {
5407         g[i] = 0.;
5408     }
5409
5410     kiwa = 0;
5411     for (jj = 0; jj < * (scs_imp->nzord); jj++)
5412     {
5413         kf = &scs_imp->zord[jj];
5414         C2F(curblk).kfun = *kf;
5415         /* continuous state */
5416         if (scs_imp->blocks[*kf - 1].nx > 0)
5417         {
5418             scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5419             scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5420         }
5421         scs_imp->blocks[*kf - 1].nevprt = scs_imp->zord[jj + * (scs_imp->nzord)];
5422
5423         if (scs_imp->funtyp[*kf - 1] > -1)
5424         {
5425             flag = 1;
5426             /*callf(told, xtd, xt, xt,xt,&flag);*/
5427             callf(told, &scs_imp->blocks[*kf - 1], &flag);
5428             if (flag < 0)
5429             {
5430                 *ierr = 5 - flag;
5431                 return;
5432             }
5433         }
5434
5435         /* Initialize tvec */
5436         if (scs_imp->blocks[*kf - 1].nevout > 0)
5437         {
5438             if (scs_imp->funtyp[*kf - 1] < 0)
5439             {
5440                 if (phase == 1 || scs_imp->blocks[*kf - 1].nmode == 0)
5441                 {
5442                     i = synchro_nev(scs_imp, *kf, ierr);
5443                     if (*ierr != 0)
5444                     {
5445                         return;
5446                     }
5447                 }
5448                 else
5449                 {
5450                     i = scs_imp->blocks[*kf - 1].mode[0];
5451                 }
5452                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5453                 putevs(told, &i2, &ierr1);
5454                 if (ierr1 != 0)
5455                 {
5456                     /* event conflict */
5457                     *ierr = 3;
5458                     return;
5459                 }
5460                 ozdoit(told, xt, xtd, &kiwa);
5461                 if (*ierr != 0)
5462                 {
5463                     return;
5464                 }
5465             }
5466         }
5467     }
5468
5469     /* update zero crossing surfaces */
5470     for (ii = 0; ii < * (scs_imp->nzord); ii++)
5471     {
5472         kf = &scs_imp->zord[ii];
5473         C2F(curblk).kfun = *kf;
5474         if (scs_imp->blocks[*kf - 1].ng > 0)
5475         {
5476             /* update g array ptr */
5477             scs_imp->blocks[*kf - 1].g = &g[scs_imp->zcptr[*kf - 1] - 1];
5478             if (scs_imp->funtyp[*kf - 1] > 0)
5479             {
5480                 flag = 9;
5481                 /* continuous state */
5482                 if (scs_imp->blocks[*kf - 1].nx > 0)
5483                 {
5484                     scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5485                     scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5486                 }
5487                 scs_imp->blocks[*kf - 1].nevprt = scs_imp->zord[ii + * (scs_imp->nzord)];
5488                 /*callf(told, xtd, xt, xtd,g,&flag);*/
5489                 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5490                 if (flag < 0)
5491                 {
5492                     *ierr = 5 - flag;
5493                     return;
5494                 }
5495             }
5496             else
5497             {
5498                 j = synchro_g_nev(scs_imp, g, *kf, ierr);
5499                 if (*ierr != 0)
5500                 {
5501                     return;
5502                 }
5503                 if ( (phase == 1) && (scs_imp->blocks[*kf - 1].nmode > 0) )
5504                 {
5505                     scs_imp->blocks[*kf - 1].mode[0] = j;
5506                 }
5507             }
5508
5509             // scs_imp->blocks[*kf-1].g = &scs_imp->g[scs_imp->zcptr[*kf]-1];
5510
5511         }
5512     }
5513
5514     for (i = 0; i < kiwa; i++)
5515     {
5516         keve = iwa[i];
5517         for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1; ii++)
5518         {
5519             kf = &scs_imp->ordclk[ii - 1];
5520             C2F(curblk).kfun = *kf;
5521             if (scs_imp->blocks[*kf - 1].ng > 0)
5522             {
5523                 /* update g array ptr */
5524                 scs_imp->blocks[*kf - 1].g = &g[scs_imp->zcptr[*kf - 1] - 1];
5525                 if (scs_imp->funtyp[*kf - 1] > 0)
5526                 {
5527                     flag = 9;
5528                     /* continuous state */
5529                     if (scs_imp->blocks[*kf - 1].nx > 0)
5530                     {
5531                         scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5532                         scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5533                     }
5534                     scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5535                     /*callf(told, xtd, xt, xtd,g,&flag);*/
5536                     callf(told, &scs_imp->blocks[*kf - 1], &flag);
5537                     if (flag < 0)
5538                     {
5539                         *ierr = 5 - flag;
5540                         return;
5541                     }
5542                 }
5543                 else
5544                 {
5545                     j = synchro_g_nev(scs_imp, g, *kf, ierr);
5546                     if (*ierr != 0)
5547                     {
5548                         return;
5549                     }
5550                     if ((phase == 1) && (scs_imp->blocks[*kf - 1].nmode > 0))
5551                     {
5552                         scs_imp->blocks[*kf - 1].mode[0] = j;
5553                     }
5554                 }
5555
5556                 //scs_imp->blocks[*kf-1].g = &scs_imp->g[scs_imp->zcptr[*kf]-1];
5557             }
5558         }
5559     }
5560 } /* zdoit_ */
5561 /*--------------------------------------------------------------------------*/
5562 /* Subroutine Jdoit */
5563 void Jdoit(double *told, double *xt, double *xtd, double *residual, int *job)
5564 {
5565     /* update blocks jacobian of continuous time block  */
5566     /*     Copyright INRIA */
5567
5568     int i2 = 0;
5569     scicos_flag flag = 0;
5570     int keve = 0, kiwa = 0;
5571     int ierr1 = 0, i = 0;
5572     int ii = 0, jj = 0;
5573
5574     ScicosImport *scs_imp = NULL;
5575     int *kf = NULL;
5576
5577     scs_imp = getscicosimportptr();
5578
5579     /* Function Body */
5580     kiwa = 0;
5581     for (jj = 0; jj < * (scs_imp->noord); jj++)
5582     {
5583         kf = &scs_imp->oord[jj];
5584         C2F(curblk).kfun = *kf;
5585         scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[jj + * (scs_imp->noord)];
5586         if (scs_imp->funtyp[*kf - 1] > -1)
5587         {
5588             flag = 1;
5589             /* applying desired output */
5590             if ((*job == 2) && (scs_imp->oord[jj] == AJacobian_block))
5591             {
5592             }
5593             else
5594                 /* continuous state */
5595                 if (scs_imp->blocks[*kf - 1].nx > 0)
5596                 {
5597                     scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5598                     scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5599                     scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5600                 }
5601
5602             /*callf(told, xtd, xt, residual,x,&flag);*/
5603             callf(told, &scs_imp->blocks[*kf - 1], &flag);
5604             if (flag < 0)
5605             {
5606                 *ierr = 5 - flag;
5607                 return;
5608             }
5609         }
5610
5611         if (scs_imp->blocks[*kf - 1].nevout > 0)
5612         {
5613             if (scs_imp->funtyp[*kf - 1] < 0)
5614             {
5615                 if (scs_imp->blocks[*kf - 1].nmode > 0)
5616                 {
5617                     i2 = scs_imp->blocks[*kf - 1].mode[0] + scs_imp->clkptr[*kf - 1] - 1;
5618                 }
5619                 else
5620                 {
5621                     i = synchro_nev(scs_imp, *kf, ierr);
5622                     if (*ierr != 0)
5623                     {
5624                         return;
5625                     }
5626                     i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5627                 }
5628                 putevs(told, &i2, &ierr1);
5629                 if (ierr1 != 0)
5630                 {
5631                     /* event conflict */
5632                     *ierr = 3;
5633                     return;
5634                 }
5635                 ozdoit(told, xt, xtd, &kiwa);
5636             }
5637         }
5638     }
5639
5640     /* update states derivatives */
5641     for (ii = 0; ii < * (scs_imp->noord); ii++)
5642     {
5643         kf = &scs_imp->oord[ii];
5644         C2F(curblk).kfun = *kf;
5645         if (scs_imp->blocks[*kf - 1].nx > 0 || \
5646                 *scs_imp->blocks[*kf - 1].work != NULL)
5647         {
5648             /* work tests if a hidden state exists, used for delay block */
5649             flag = 0;
5650             if (((*job == 1) && (scs_imp->oord[ii] == AJacobian_block)) || (*job != 1))
5651             {
5652                 if (*job == 1)
5653                 {
5654                     flag = 10;
5655                 }
5656                 /* continuous state */
5657                 if (scs_imp->blocks[*kf - 1].nx > 0)
5658                 {
5659                     scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5660                     scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5661                     scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5662                 }
5663                 scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[ii + * (scs_imp->noord)];
5664                 /*callf(told, xtd, xt, residual,xt,&flag);*/
5665                 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5666             }
5667             if (flag < 0)
5668             {
5669                 *ierr = 5 - flag;
5670                 return;
5671             }
5672         }
5673     }
5674
5675     for (i = 0; i < kiwa; i++)
5676     {
5677         keve = iwa[i];
5678         for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1;  ii++)
5679         {
5680             kf = &scs_imp->ordclk[ii - 1];
5681             C2F(curblk).kfun = *kf;
5682             if (scs_imp->blocks[*kf - 1].nx > 0 || \
5683                     *scs_imp->blocks[*kf - 1].work != NULL)
5684             {
5685                 /* work tests if a hidden state exists */
5686                 flag = 0;
5687                 if (((*job == 1) && (scs_imp->oord[ii - 1] == AJacobian_block)) || (*job != 1))
5688                 {
5689                     if (*job == 1)
5690                     {
5691                         flag = 10;
5692                     }
5693                     /* continuous state */
5694                     if (scs_imp->blocks[*kf - 1].nx > 0)
5695                     {
5696                         scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5697                         scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5698                         scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5699                     }
5700                     scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5701                     /*callf(told, xtd, xt, residual,xt,&flag);*/
5702                     callf(told, &scs_imp->blocks[*kf - 1], &flag);
5703                 }
5704                 if (flag < 0)
5705                 {
5706                     *ierr = 5 - flag;
5707                     return;
5708                 }
5709             }
5710         }
5711     }
5712 } /* Jdoit_ */
5713 /*--------------------------------------------------------------------------*/
5714 /* Subroutine synchro_nev */
5715 static int synchro_nev(ScicosImport *scs_imp, int kf, int *ierr)
5716 {
5717     /* synchro blocks computation  */
5718     /*     Copyright INRIA */
5719     SCSREAL_COP *outtbdptr = NULL;     /*to store double of outtb*/
5720     SCSINT8_COP *outtbcptr = NULL;     /*to store int8 of outtb*/
5721     SCSINT16_COP *outtbsptr = NULL;    /*to store int16 of outtb*/
5722     SCSINT32_COP *outtblptr = NULL;    /*to store int32 of outtb*/
5723     SCSUINT8_COP *outtbucptr = NULL;   /*to store unsigned int8 of outtb */
5724     SCSUINT16_COP *outtbusptr = NULL;  /*to store unsigned int16 of outtb */
5725     SCSUINT32_COP *outtbulptr = NULL;  /*to store unsigned int32 of outtb */
5726
5727     int cond = 0;
5728     int i = 0; /* return 0 by default */
5729
5730     /* variable for param */
5731     int *outtbtyp = 0;
5732     void **outtbptr = NULL;
5733     int *funtyp = 0;
5734     int *inplnk = 0;
5735     int *inpptr = 0;
5736
5737     /* get param ptr */
5738     outtbtyp = scs_imp->outtbtyp;
5739     outtbptr = scs_imp->outtbptr;
5740     funtyp = scs_imp->funtyp;
5741     inplnk = scs_imp->inplnk;
5742     inpptr = scs_imp->inpptr;
5743
5744     /* if-then-else blk */
5745     if (funtyp[kf - 1] == -1)
5746     {
5747         switch (outtbtyp[-1 + inplnk[inpptr[kf - 1] - 1]])
5748         {
5749             case SCSREAL_N    :
5750                 outtbdptr = (SCSREAL_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5751                 cond = (*outtbdptr <= 0.);
5752                 break;
5753
5754             case SCSCOMPLEX_N :
5755                 outtbdptr = (SCSCOMPLEX_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5756                 cond = (*outtbdptr <= 0.);
5757                 break;
5758
5759             case SCSINT8_N    :
5760                 outtbcptr = (SCSINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5761                 cond = (*outtbcptr <= 0);
5762                 break;
5763
5764             case SCSINT16_N   :
5765                 outtbsptr = (SCSINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5766                 cond = (*outtbsptr <= 0);
5767                 break;
5768
5769             case SCSINT32_N   :
5770                 outtblptr = (SCSINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5771                 cond = (*outtblptr <= 0);
5772                 break;
5773
5774             case SCSUINT8_N   :
5775                 outtbucptr = (SCSUINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5776                 cond = (*outtbucptr <= 0);
5777                 break;
5778
5779             case SCSUINT16_N  :
5780                 outtbusptr = (SCSUINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5781                 cond = (*outtbusptr <= 0);
5782                 break;
5783
5784             case SCSUINT32_N  :
5785                 outtbulptr = (SCSUINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5786                 cond = (*outtbulptr <= 0);
5787                 break;
5788
5789             default  : /* Add a message here */
5790                 *ierr = 25;
5791                 return 0;
5792         }
5793         if (cond)
5794         {
5795             i = 2;
5796         }
5797         else
5798         {
5799             i = 1;
5800         }
5801     }
5802     /* eselect blk */
5803     else if (funtyp[kf - 1] == -2)
5804     {
5805         switch (outtbtyp[-1 + inplnk[inpptr[kf - 1] - 1]])
5806         {
5807             case SCSREAL_N    :
5808                 outtbdptr = (SCSREAL_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5809                 i = Max(Min((int) * outtbdptr, scs_imp->blocks[kf - 1].nevout), 1);
5810                 break;
5811
5812             case SCSCOMPLEX_N :
5813                 outtbdptr = (SCSCOMPLEX_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5814                 i = Max(Min((int) * outtbdptr, scs_imp->blocks[kf - 1].nevout), 1);
5815                 break;
5816
5817             case SCSINT8_N    :
5818                 outtbcptr = (SCSINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5819                 i = Max(Min((int) * outtbcptr,
5820                             scs_imp->blocks[kf - 1].nevout), 1);