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