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