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