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