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