Scicos: only call CVDense() when necessary
[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         /* Call CVDense to specify the CVDENSE dense linear solver, only for solvers needing CVode's Newton method */
1622         if (solver == CVode_BDF_Newton || solver == CVode_Adams_Newton)
1623         {
1624             flag = CVDense(ode_mem, *neq);
1625         }
1626         if (check_flag(&flag, "CVDense", 1))
1627         {
1628             *ierr = 300 + (-flag);
1629             freeall
1630             return;
1631         }
1632
1633         if (hmax > 0)
1634         {
1635             flag = ODESetMaxStep(ode_mem, (realtype) hmax);
1636             if (check_flag(&flag, "CVodeSetMaxStep", 1))
1637             {
1638                 *ierr = 300 + (-flag);
1639                 freeall;
1640                 return;
1641             }
1642         }
1643         /* Set the Jacobian routine to Jac (user-supplied)
1644         flag = CVDlsSetDenseJacFn(ode_mem, Jac);
1645         if (check_flag(&flag, "CVDlsSetDenseJacFn", 1)) return(1);  */
1646
1647     }/* testing if neq>0 */
1648
1649     /* Function Body */
1650     C2F(coshlt).halt = 0;
1651     *ierr = 0;
1652
1653     C2F(xscion)(&inxsci);
1654     /*     initialization */
1655     C2F(realtimeinit)(told, &C2F(rtfactor).scale);
1656
1657     phase = 1;
1658     hot = 0;
1659
1660     jj = 0;
1661     for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
1662     {
1663         if (Blocks[C2F(curblk).kfun - 1].ng >= 1)
1664         {
1665             zcros[jj] = C2F(curblk).kfun;
1666             ++jj;
1667         }
1668     }
1669     /*     . ng >= jj required */
1670     if (jj != ng)
1671     {
1672         zcros[jj] = -1;
1673     }
1674     /*     initialization (propagation of constant blocks outputs) */
1675     idoit(told);
1676     if (*ierr != 0)
1677     {
1678         freeall;
1679         return;
1680     }
1681     /*--discrete zero crossings----dzero--------------------*/
1682     if (ng > 0) /* storing ZC signs just after a solver call*/
1683     {
1684         /*zdoit(told, g, x, x);*/
1685         zdoit(told, x, x, g);
1686         if (*ierr != 0)
1687         {
1688             freeall;
1689             return;
1690         }
1691         for (jj = 0; jj < ng; ++jj)
1692             if (g[jj] >= 0)
1693             {
1694                 jroot[jj] = 5;
1695             }
1696             else
1697             {
1698                 jroot[jj] = -5;
1699             }
1700     }
1701     /*--discrete zero crossings----dzero--------------------*/
1702
1703     /*     main loop on time */
1704     while (*told < *tf)
1705     {
1706         while (ismenu()) //** if the user has done something, do the actions
1707         {
1708             int ierr2 = 0;
1709             SeqSync = GetCommand(CommandToUnstack); //** get to the action
1710             CommandLength = (int)strlen(CommandToUnstack);
1711             syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
1712         }
1713         if (C2F(coshlt).halt != 0)
1714         {
1715             if (C2F(coshlt).halt == 2)
1716             {
1717                 *told = *tf;    /* end simulation */
1718             }
1719             C2F(coshlt).halt = 0;
1720             freeall;
1721             return;
1722         }
1723         if (*pointi == 0)
1724         {
1725             t = *tf;
1726         }
1727         else
1728         {
1729             t = tevts[*pointi];
1730         }
1731         if (fabs(t - *told) < ttol)
1732         {
1733             t = *told;
1734             /*     update output part */
1735         }
1736         if (*told > t)
1737         {
1738             /*     !  scheduling problem */
1739             *ierr = 1;
1740             freeall;
1741             return;
1742         }
1743         if (*told != t)
1744         {
1745             if (xptr[nblk + 1] == 1)
1746             {
1747                 /*     .     no continuous state */
1748                 if (*told + deltat + ttol > t)
1749                 {
1750                     *told = t;
1751                 }
1752                 else
1753                 {
1754                     *told += deltat;
1755                 }
1756                 /*     .     update outputs of 'c' type blocks with no continuous state */
1757                 if (*told >= *tf)
1758                 {
1759                     /*     .     we are at the end, update continuous part before leaving */
1760                     if (ncord > 0)
1761                     {
1762                         cdoit(told);
1763                         freeall;
1764                         return;
1765                     }
1766                 }
1767             }
1768             else
1769             {
1770                 /*     integrate */
1771                 rhotmp = *tf + ttol;
1772                 if (*pointi != 0)
1773                 {
1774                     kpo = *pointi;
1775 L20:
1776                     if (critev[kpo] == 1)
1777                     {
1778                         rhotmp = tevts[kpo];
1779                         goto L30;
1780                     }
1781                     kpo = evtspt[kpo];
1782                     if (kpo != 0)
1783                     {
1784                         goto L20;
1785                     }
1786 L30:
1787                     if (rhotmp < tstop)
1788                     {
1789                         hot = 0;
1790                     }
1791                 }
1792                 tstop = rhotmp;
1793                 t = Min(*told + deltat, Min(t, *tf + ttol));
1794
1795                 if (ng > 0 &&  hot == 0 && nmod > 0)
1796                 {
1797                     zdoit(told, x, x, g);
1798                     if (*ierr != 0)
1799                     {
1800                         freeall;
1801                         return;
1802                     }
1803                 }
1804
1805                 if (hot == 0) /* hot==0 : cold restart*/
1806                 {
1807                     flag = ODESetStopTime(ode_mem, (realtype)tstop);  /* Setting the stop time*/
1808                     if (check_flag(&flag, "CVodeSetStopTime", 1))
1809                     {
1810                         *ierr = 300 + (-flag);
1811                         freeall;
1812                         return;
1813                     }
1814
1815                     flag = ODEReInit(ode_mem, (realtype)(*told), y);
1816                     if (check_flag(&flag, "CVodeReInit", 1))
1817                     {
1818                         *ierr = 300 + (-flag);
1819                         freeall;
1820                         return;
1821                     }
1822                 }
1823
1824                 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1825                 {
1826                     sciprint(_("****SUNDIALS.Cvode from: %f to %f hot= %d  \n"), *told, t, hot);
1827                 }
1828
1829                 /*--discrete zero crossings----dzero--------------------*/
1830                 /*--check for Dzeros after Mode settings or ddoit()----*/
1831                 Discrete_Jump = 0;
1832
1833                 if (ng > 0 && hot == 0)
1834                 {
1835                     zdoit(told, x, x, g);
1836                     if (*ierr != 0)
1837                     {
1838                         freeall;
1839                         return;
1840                     }
1841                     for (jj = 0; jj < ng; ++jj)
1842                     {
1843                         if ((g[jj] >= 0.0) && (jroot[jj] == -5))
1844                         {
1845                             Discrete_Jump = 1;
1846                             jroot[jj] = 1;
1847                         }
1848                         else if ((g[jj] < 0.0) && (jroot[jj] == 5))
1849                         {
1850                             Discrete_Jump = 1;
1851                             jroot[jj] = -1;
1852                         }
1853                         else
1854                         {
1855                             jroot[jj] = 0;
1856                         }
1857                     }
1858                 }
1859                 /*--discrete zero crossings----dzero--------------------*/
1860
1861                 if (Discrete_Jump == 0) /* if there was a dzero, its event should be activated*/
1862                 {
1863                     phase = 2;
1864                     flag = ODE(ode_mem, t, y, told, ODE_NORMAL);
1865                     if (*ierr != 0)
1866                     {
1867                         freeall;
1868                         return;
1869                     }
1870                     phase = 1;
1871                 }
1872                 else
1873                 {
1874                     flag = CV_ROOT_RETURN; /* in order to handle discrete jumps */
1875                 }
1876
1877                 /*     .     update outputs of 'c' type  blocks if we are at the end*/
1878                 if (*told >= *tf)
1879                 {
1880                     if (ncord > 0)
1881                     {
1882                         cdoit(told);
1883                         freeall;
1884                         return;
1885                     }
1886                 }
1887
1888                 if (flag >= 0)
1889                 {
1890                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1891                     {
1892                         sciprint(_("****SUNDIALS.Cvode reached: %f\n"), *told);
1893                     }
1894                     hot = 1;
1895                     cnt = 0;
1896                 }
1897                 else if ( flag == CV_TOO_MUCH_WORK ||  flag == CV_CONV_FAILURE || flag == CV_ERR_FAILURE)
1898                 {
1899                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1900                     {
1901                         sciprint(_("****SUNDIALS.Cvode: too much work at time=%g (stiff region, change RTOL and ATOL)\n"), *told);
1902                     }
1903                     hot = 0;
1904                     cnt++;
1905                     if (cnt > 5)
1906                     {
1907                         *ierr = 300 + (-flag);
1908                         freeall;
1909                         return;
1910                     }
1911                 }
1912                 else
1913                 {
1914                     if (flag < 0)
1915                     {
1916                         *ierr = 300 + (-flag);    /* raising errors due to internal errors, otherwise error due to flagr*/
1917                     }
1918                     freeall;
1919                     return;
1920                 }
1921
1922                 if (flag == CV_ZERO_DETACH_RETURN)
1923                 {
1924                     hot = 0;
1925                 };  /* new feature of sundials, detects zero-detaching */
1926
1927                 if (flag == CV_ROOT_RETURN)
1928                 {
1929                     /*     .        at a least one root has been found */
1930                     hot = 0;
1931                     if (Discrete_Jump == 0)
1932                     {
1933                         flagr = ODEGetRootInfo(ode_mem, jroot);
1934                         if (check_flag(&flagr, "CVodeGetRootInfo", 1))
1935                         {
1936                             *ierr = 300 + (-flagr);
1937                             freeall;
1938                             return;
1939                         }
1940                     }
1941                     /*     .        at a least one root has been found */
1942                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1943                     {
1944                         sciprint(_("root found at t=: %f\n"), *told);
1945                     }
1946                     /*     .        update outputs affecting ztyp blocks ONLY FOR OLD BLOCKS */
1947                     zdoit(told, x, xd, g);
1948                     if (*ierr != 0)
1949                     {
1950                         freeall;
1951                         return;
1952                     }
1953                     for (jj = 0; jj < ng; ++jj)
1954                     {
1955                         C2F(curblk).kfun = zcros[ jj];
1956                         if (C2F(curblk).kfun == -1)
1957                         {
1958                             break;
1959                         }
1960                         kev = 0;
1961
1962                         for (j = zcptr[C2F(curblk).kfun] - 1 ;
1963                                 j < zcptr[C2F(curblk).kfun + 1] - 1 ; ++j)
1964                         {
1965                             if (jroot[j] != 0)
1966                             {
1967                                 kev = 1;
1968                                 break;
1969                             }
1970                         }
1971                         /*   */
1972                         if (kev != 0)
1973                         {
1974                             Blocks[C2F(curblk).kfun - 1].jroot = &jroot[zcptr[C2F(curblk).kfun] - 1];
1975                             if (funtyp[C2F(curblk).kfun] > 0)
1976                             {
1977
1978                                 if (Blocks[C2F(curblk).kfun - 1].nevout > 0)
1979                                 {
1980                                     flag__ = 3;
1981                                     if (Blocks[C2F(curblk).kfun - 1].nx > 0)
1982                                     {
1983                                         Blocks[C2F(curblk).kfun - 1].x  = &x[xptr[C2F(curblk).kfun] - 1];
1984                                         Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
1985                                     }
1986                                     /* call corresponding block to determine output event (kev) */
1987                                     Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
1988                                     /*callf(told, xd, x, x,g,&flag__);*/
1989                                     callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1990                                     if (flag__ < 0)
1991                                     {
1992                                         *ierr = 5 - flag__;
1993                                         freeall;
1994                                         return;
1995                                     }
1996                                     /*     .              update event agenda */
1997                                     for (k = 0; k < Blocks[C2F(curblk).kfun - 1].nevout; ++k)
1998                                     {
1999                                         if (Blocks[C2F(curblk).kfun - 1].evout[k] >= 0.)
2000                                         {
2001                                             i3 = k + clkptr[C2F(curblk).kfun] ;
2002                                             addevs(Blocks[C2F(curblk).kfun - 1].evout[k] + (*told), &i3, &ierr1);
2003                                             if (ierr1 != 0)
2004                                             {
2005                                                 /*     .                       nevts too small */
2006                                                 *ierr = 3;
2007                                                 freeall;
2008                                                 return;
2009                                             }
2010                                         }
2011                                     }
2012                                 }
2013                                 /*     .              update state */
2014                                 if (Blocks[C2F(curblk).kfun - 1].nx > 0)
2015                                 {
2016                                     /*     .              call corresponding block to update state */
2017                                     flag__ = 2;
2018                                     Blocks[C2F(curblk).kfun - 1].x  = &x[xptr[C2F(curblk).kfun] - 1];
2019                                     Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
2020                                     Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
2021                                     /*callf(told, xd, x, x,g,&flag__);*/
2022                                     callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
2023
2024                                     if (flag__ < 0)
2025                                     {
2026                                         *ierr = 5 - flag__;
2027                                         freeall;
2028                                         return;
2029                                     }
2030                                 }
2031                             }
2032                         }
2033                     }
2034                 }
2035             }
2036             /*--discrete zero crossings----dzero--------------------*/
2037             if (ng > 0) /* storing ZC signs just after a sundials call*/
2038             {
2039                 zdoit(told, x, x, g);
2040                 if (*ierr != 0)
2041                 {
2042                     freeall;
2043                     return;
2044                 }
2045                 for (jj = 0; jj < ng; ++jj)
2046                 {
2047                     if (g[jj] >= 0)
2048                     {
2049                         jroot[jj] = 5;
2050                     }
2051                     else
2052                     {
2053                         jroot[jj] = -5;
2054                     }
2055                 }
2056             }
2057             /*--discrete zero crossings----dzero--------------------*/
2058
2059             C2F(realtime)(told);
2060         }
2061         else
2062         {
2063             /*     .  t==told */
2064             if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2065             {
2066                 sciprint(_("Event: %d activated at t=%f\n"), *pointi, *told);
2067                 for (kev = 0; kev < nblk; kev++)
2068                 {
2069                     if (Blocks[kev].nmode > 0)
2070                     {
2071                         sciprint(_("mode of block %d=%d, "), kev, Blocks[kev].mode[0]);
2072                     }
2073                 }
2074                 sciprint(_("**mod**\n"));
2075             }
2076
2077             ddoit(told);
2078             if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2079             {
2080                 sciprint(_("End of activation\n"));
2081             }
2082             if (*ierr != 0)
2083             {
2084                 freeall;
2085                 return;
2086             }
2087
2088         }
2089         /*     end of main loop on time */
2090     }
2091     freeall;
2092 } /* cossim_ */
2093
2094 /*--------------------------------------------------------------------------*/
2095 static void cossimdaskr(double *told)
2096 {
2097     /* System generated locals */
2098     int i3;
2099     //** used for the [stop] button
2100     static char CommandToUnstack[1024];
2101     static int CommandLength = 0;
2102     static int SeqSync = 0;
2103     static int one = 1;
2104
2105     /* Local variables */
2106     static scicos_flag flag__ = 0;
2107     static int ierr1 = 0;
2108     static int j = 0, k = 0;
2109     static double t = 0.;
2110     static int jj = 0;
2111     static double rhotmp = 0., tstop = 0.;
2112     static int inxsci = 0;
2113     static int kpo = 0, kev = 0;
2114
2115     int *jroot = NULL, *zcros = NULL;
2116     int *Mode_save = NULL;
2117     int Mode_change = 0;
2118
2119     int flag = 0, flagr = 0;
2120     N_Vector   yy = NULL, yp = NULL;
2121     realtype reltol = 0., abstol = 0.;
2122     int Discrete_Jump = 0;
2123     N_Vector IDx = NULL;
2124     realtype *scicos_xproperty = NULL;
2125     DlsMat TJacque = NULL;
2126
2127     void *dae_mem = NULL;
2128     UserData data = NULL;
2129     IDAMem copy_IDA_mem = NULL;
2130     int maxnj = 0, maxnit = 0, maxnh = 0;
2131     /*-------------------- Analytical Jacobian memory allocation ----------*/
2132     int  Jn = 0, Jnx = 0, Jno = 0, Jni = 0, Jactaille = 0;
2133     double uround = 0.;
2134     int cnt = 0, N_iters = 0;
2135     /* Saving solver number */
2136     int solver = C2F(cmsolver).solver;
2137     /* Flags for initial values calculation */
2138     int DAE_YA_YDP_INIT = 1;
2139     int DAE_Y_INIT      = 2;
2140     /* Defining function pointers, for more readability*/
2141     void(* DAEFree) (void**);
2142     int (* DAESolve) (void*, realtype, realtype*, N_Vector, N_Vector, int);
2143     int (* DAEReInit) (void*, realtype, N_Vector, N_Vector);
2144     int (* DAESetId) (void*, N_Vector);
2145     int (* DAECalcIC) (void*, int, realtype);
2146     int (* DAESetMaxStep) (void*, realtype);
2147     int (* DAESetUserData) (void*, void*);
2148     int (* DAESetStopTime) (void*, realtype);
2149     int (* DAEGetRootInfo) (void*, int*);
2150     int (* DAESStolerances) (void*, realtype, realtype);
2151     int (* DAEGetConsistentIC) (void*, N_Vector, N_Vector);
2152     int (* DAESetMaxNumSteps) (void*, long int);
2153     int (* DAESetMaxNumJacsIC) (void*, int);
2154     int (* DAESetMaxNumItersIC) (void*, int);
2155     int (* DAESetMaxNumStepsIC) (void*, int);
2156     int (* DAESetLineSearchOffIC) (void*, int);
2157     /* For DAEs, the generic flags for stop mode depend on the used solver */
2158     int DAE_NORMAL = 0, DAE_ONE_STEP = 0;
2159     DAE_NORMAL   = (solver == IDA_BDF_Newton) ? 1 : 0;  /* IDA_NORMAL   = 1, DDAS_NORMAL   = 0 */
2160     DAE_ONE_STEP = (solver == IDA_BDF_Newton) ? 2 : 1;  /* IDA_ONE_STEP = 2, DDAS_ONE_STEP = 1 */
2161     switch (solver)
2162     {
2163         case IDA_BDF_Newton:
2164             DAEFree = &IDAFree;
2165             DAESolve = &IDASolve;
2166             DAESetId = &IDASetId;
2167             DAEReInit = &IDAReInit;
2168             DAECalcIC = &IDACalcIC;
2169             DAESetMaxStep = &IDASetMaxStep;
2170             DAESetUserData = &IDASetUserData;
2171             DAESetStopTime = &IDASetStopTime;
2172             DAEGetRootInfo = &IDAGetRootInfo;
2173             DAESStolerances = &IDASStolerances;
2174             DAESetMaxNumSteps = &IDASetMaxNumSteps;
2175             DAEGetConsistentIC = &IDAGetConsistentIC;
2176             DAESetMaxNumJacsIC = &IDASetMaxNumJacsIC;
2177             DAESetMaxNumItersIC = &IDASetMaxNumItersIC;
2178             DAESetMaxNumStepsIC = &IDASetMaxNumStepsIC;
2179             DAESetLineSearchOffIC = &IDASetLineSearchOffIC;
2180             break;
2181         case DDaskr_BDF_Newton:
2182         case DDaskr_BDF_GMRes:
2183             DAEFree = &DDaskrFree;
2184             DAESolve = &DDaskrSolve;
2185             DAESetId = &DDaskrSetId;
2186             DAEReInit = &DDaskrReInit;
2187             DAECalcIC = &DDaskrCalcIC;
2188             DAESetMaxStep = &DDaskrSetMaxStep;
2189             DAESetUserData = &DDaskrSetUserData;
2190             DAESetStopTime = &DDaskrSetStopTime;
2191             DAEGetRootInfo = &DDaskrGetRootInfo;
2192             DAESStolerances = &DDaskrSStolerances;
2193             DAESetMaxNumSteps = &DDaskrSetMaxNumSteps;
2194             DAEGetConsistentIC = &DDaskrGetConsistentIC;
2195             DAESetMaxNumJacsIC = &DDaskrSetMaxNumJacsIC;
2196             DAESetMaxNumItersIC = &DDaskrSetMaxNumItersIC;
2197             DAESetMaxNumStepsIC = &DDaskrSetMaxNumStepsIC;
2198             DAESetLineSearchOffIC = &DDaskrSetLineSearchOffIC;
2199             break;
2200         default: // Unknown solver number
2201             *ierr = 1000;
2202             return;
2203     }
2204
2205     /* Set extension of Sundials for scicos */
2206     set_sundials_with_extension(TRUE);
2207
2208     // CI=1.0;   /* for function Get_Jacobian_ci */
2209     jroot = NULL;
2210     if (ng > 0)
2211     {
2212         if ((jroot = MALLOC(sizeof(int) * ng)) == NULL )
2213         {
2214             *ierr = 10000;
2215             return;
2216         }
2217     }
2218     for ( jj = 0 ; jj < ng ; jj++ )
2219     {
2220         jroot[jj] = 0 ;
2221     }
2222
2223     zcros = NULL;
2224     if (ng > 0)
2225     {
2226         if ((zcros = MALLOC(sizeof(int) * ng)) == NULL )
2227         {
2228             *ierr = 10000;
2229             if (ng > 0)
2230             {
2231                 FREE(jroot);
2232             }
2233             return;
2234         }
2235     }
2236
2237     Mode_save = NULL;
2238     if (nmod > 0)
2239     {
2240         if ((Mode_save = MALLOC(sizeof(int) * nmod)) == NULL )
2241         {
2242             *ierr = 10000;
2243             if (ng > 0)
2244             {
2245                 FREE(jroot);
2246             }
2247             if (ng > 0)
2248             {
2249                 FREE(zcros);
2250             }
2251             return;
2252         }
2253     }
2254
2255     reltol = (realtype) rtol;
2256     abstol = (realtype) Atol;  /*  Ith(abstol,1) = (realtype) Atol;*/
2257
2258     if (*neq > 0)
2259     {
2260         yy = NULL;
2261         yy = N_VNewEmpty_Serial(*neq);
2262         if (check_flag((void *)yy, "N_VNew_Serial", 0))
2263         {
2264             if (ng > 0)
2265             {
2266                 FREE(jroot);
2267             }
2268             if (ng > 0)
2269             {
2270                 FREE(zcros);
2271             }
2272             if (nmod > 0)
2273             {
2274                 FREE(Mode_save);
2275             }
2276         }
2277         NV_DATA_S(yy) = x;
2278
2279         yp = NULL;
2280         yp = N_VNewEmpty_Serial(*neq);
2281         if (check_flag((void *)yp, "N_VNew_Serial", 0))
2282         {
2283             if (*neq > 0)
2284             {
2285                 N_VDestroy_Serial(yy);
2286             }
2287             if (ng > 0)
2288             {
2289                 FREE(jroot);
2290             }
2291             if (ng > 0)
2292             {
2293                 FREE(zcros);
2294             }
2295             if (nmod > 0)
2296             {
2297                 FREE(Mode_save);
2298             }
2299             return;
2300         }
2301         NV_DATA_S(yp) = xd;
2302
2303         IDx = NULL;
2304         IDx = N_VNew_Serial(*neq);
2305         if (check_flag((void *)IDx, "N_VNew_Serial", 0))
2306         {
2307             *ierr = 10000;
2308             if (*neq > 0)
2309             {
2310                 N_VDestroy_Serial(yp);
2311             }
2312             if (*neq > 0)
2313             {
2314                 N_VDestroy_Serial(yy);
2315             }
2316             if (ng > 0)
2317             {
2318                 FREE(jroot);
2319             }
2320             if (ng > 0)
2321             {
2322                 FREE(zcros);
2323             }
2324             if (nmod > 0)
2325             {
2326                 FREE(Mode_save);
2327             }
2328             return;
2329         }
2330
2331         /* Call the Create and Init functions to initialize DAE memory */
2332         dae_mem = NULL;
2333         if (solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes)
2334         {
2335             dae_mem = DDaskrCreate(neq, ng, solver);
2336         }
2337         else
2338         {
2339             dae_mem = IDACreate();
2340         }
2341         if (check_flag((void *)dae_mem, "IDACreate", 0))
2342         {
2343             if (*neq > 0)
2344             {
2345                 N_VDestroy_Serial(IDx);
2346             }
2347             if (*neq > 0)
2348             {
2349                 N_VDestroy_Serial(yp);
2350             }
2351             if (*neq > 0)
2352             {
2353                 N_VDestroy_Serial(yy);
2354             }
2355             if (ng > 0)
2356             {
2357                 FREE(jroot);
2358             }
2359             if (ng > 0)
2360             {
2361                 FREE(zcros);
2362             }
2363             if (nmod > 0)
2364             {
2365                 FREE(Mode_save);
2366             }
2367             return;
2368         }
2369         if (C2F(cmsolver).solver == 100)
2370         {
2371             copy_IDA_mem = (IDAMem) dae_mem;
2372         }
2373
2374         if (solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes)
2375         {
2376             flag = DDaskrSetErrHandlerFn(dae_mem, SundialsErrHandler, NULL);
2377         }
2378         else
2379         {
2380             flag = IDASetErrHandlerFn(dae_mem, SundialsErrHandler, NULL);
2381         }
2382         if (check_flag(&flag, "IDASetErrHandlerFn", 1))
2383         {
2384             *ierr = 200 + (-flag);
2385             if (*neq > 0)
2386             {
2387                 DAEFree(&dae_mem);
2388             }
2389             if (*neq > 0)
2390             {
2391                 N_VDestroy_Serial(IDx);
2392             }
2393             if (*neq > 0)
2394             {
2395                 N_VDestroy_Serial(yp);
2396             }
2397             if (*neq > 0)
2398             {
2399                 N_VDestroy_Serial(yy);
2400             }
2401             if (ng > 0)
2402             {
2403                 FREE(jroot);
2404             }
2405             if (ng > 0)
2406             {
2407                 FREE(zcros);
2408             }
2409             if (nmod > 0)
2410             {
2411                 FREE(Mode_save);
2412             }
2413             return;
2414         }
2415
2416         if (solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes)
2417         {
2418             flag = DDaskrInit(dae_mem, simblkddaskr, T0, yy, yp, jacpsol, psol);
2419         }
2420         else
2421         {
2422             flag = IDAInit(dae_mem, simblkdaskr, T0, yy, yp);
2423         }
2424         if (check_flag(&flag, "IDAInit", 1))
2425         {
2426             *ierr = 200 + (-flag);
2427             if (*neq > 0)
2428             {
2429                 DAEFree(&dae_mem);
2430             }
2431             if (*neq > 0)
2432             {
2433                 N_VDestroy_Serial(IDx);
2434             }
2435             if (*neq > 0)
2436             {
2437                 N_VDestroy_Serial(yp);
2438             }
2439             if (*neq > 0)
2440             {
2441                 N_VDestroy_Serial(yy);
2442             }
2443             if (ng > 0)
2444             {
2445                 FREE(jroot);
2446             }
2447             if (ng > 0)
2448             {
2449                 FREE(zcros);
2450             }
2451             if (nmod > 0)
2452             {
2453                 FREE(Mode_save);
2454             }
2455             return;
2456         }
2457
2458         flag = DAESStolerances(dae_mem, reltol, abstol);
2459         if (check_flag(&flag, "IDASStolerances", 1))
2460         {
2461             *ierr = 200 + (-flag);
2462             if (*neq > 0)
2463             {
2464                 DAEFree(&dae_mem);
2465             }
2466             if (*neq > 0)
2467             {
2468                 N_VDestroy_Serial(IDx);
2469             }
2470             if (*neq > 0)
2471             {
2472                 N_VDestroy_Serial(yp);
2473             }
2474             if (*neq > 0)
2475             {
2476                 N_VDestroy_Serial(yy);
2477             }
2478             if (ng > 0)
2479             {
2480                 FREE(jroot);
2481             }
2482             if (ng > 0)
2483             {
2484                 FREE(zcros);
2485             }
2486             if (nmod > 0)
2487             {
2488                 FREE(Mode_save);
2489             }
2490             return;
2491         }
2492
2493         if (solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes)
2494         {
2495             flag = DDaskrRootInit(dae_mem, ng, grblkddaskr);
2496         }
2497         else
2498         {
2499             flag = IDARootInit(dae_mem, ng, grblkdaskr);
2500         }
2501         if (check_flag(&flag, "IDARootInit", 1))
2502         {
2503             *ierr = 200 + (-flag);
2504             if (*neq > 0)
2505             {
2506                 DAEFree(&dae_mem);
2507             }
2508             if (*neq > 0)
2509             {
2510                 N_VDestroy_Serial(IDx);
2511             }
2512             if (*neq > 0)
2513             {
2514                 N_VDestroy_Serial(yp);
2515             }
2516             if (*neq > 0)
2517             {
2518                 N_VDestroy_Serial(yy);
2519             }
2520             if (ng > 0)
2521             {
2522                 FREE(jroot);
2523             }
2524             if (ng > 0)
2525             {
2526                 FREE(zcros);
2527             }
2528             if (nmod > 0)
2529             {
2530                 FREE(Mode_save);
2531             }
2532             return;
2533         }
2534
2535         if (solver == IDA_BDF_Newton)
2536         {
2537             flag = IDADense(dae_mem, *neq);
2538         }
2539         if (check_flag(&flag, "IDADense", 1))
2540         {
2541             *ierr = 200 + (-flag);
2542             if (*neq > 0)if (solver == IDA_BDF_Newton)
2543                 {
2544                     IDAFree(&dae_mem);
2545                 }
2546             if (*neq > 0)
2547             {
2548                 N_VDestroy_Serial(IDx);
2549             }
2550             if (*neq > 0)
2551             {
2552                 N_VDestroy_Serial(yp);
2553             }
2554             if (*neq > 0)
2555             {
2556                 N_VDestroy_Serial(yy);
2557             }
2558             if (ng > 0)
2559             {
2560                 FREE(jroot);
2561             }
2562             if (ng > 0)
2563             {
2564                 FREE(zcros);
2565             }
2566             if (nmod > 0)
2567             {
2568                 FREE(Mode_save);
2569             }
2570             return;
2571         }
2572
2573         data = NULL;
2574         if ((data = (UserData) MALLOC(sizeof(*data))) == NULL)
2575         {
2576             *ierr = 10000;
2577             if (*neq > 0)
2578             {
2579                 IDAFree(&dae_mem);
2580             }
2581             if (*neq > 0)
2582             {
2583                 N_VDestroy_Serial(IDx);
2584             }
2585             if (*neq > 0)
2586             {
2587                 N_VDestroy_Serial(yp);
2588             }
2589             if (*neq > 0)
2590             {
2591                 N_VDestroy_Serial(yy);
2592             }
2593             if (ng > 0)
2594             {
2595                 FREE(jroot);
2596             }
2597             if (ng > 0)
2598             {
2599                 FREE(zcros);
2600             }
2601             if (nmod > 0)
2602             {
2603                 FREE(Mode_save);
2604             }
2605             return;
2606         }
2607         data->dae_mem = dae_mem;
2608         data->ewt   = NULL;
2609         data->iwork = NULL;
2610         data->rwork = NULL;
2611         data->gwork = NULL;
2612
2613         data->ewt   = N_VNew_Serial(*neq);
2614         if (check_flag((void *)data->ewt, "N_VNew_Serial", 0))
2615         {
2616             *ierr = 200 + (-flag);
2617             if (*neq > 0)
2618             {
2619                 FREE(data);
2620             }
2621             if (*neq > 0)
2622             {
2623                 DAEFree(&dae_mem);
2624             }
2625             if (*neq > 0)
2626             {
2627                 N_VDestroy_Serial(IDx);
2628             }
2629             if (*neq > 0)
2630             {
2631                 N_VDestroy_Serial(yp);
2632             }
2633             if (*neq > 0)
2634             {
2635                 N_VDestroy_Serial(yy);
2636             }
2637             if (ng > 0)
2638             {
2639                 FREE(jroot);
2640             }
2641             if (ng > 0)
2642             {
2643                 FREE(zcros);
2644             }
2645             if (nmod > 0)
2646             {
2647                 FREE(Mode_save);
2648             }
2649             return;
2650         }
2651         if ( ng > 0 )
2652         {
2653             if ((data->gwork = (double *) MALLOC(ng * sizeof(double))) == NULL)
2654             {
2655                 if (*neq > 0)
2656                 {
2657                     N_VDestroy_Serial(data->ewt);
2658                 }
2659                 if (*neq > 0)
2660                 {
2661                     FREE(data);
2662                 }
2663                 if (*neq > 0)
2664                 {
2665                     DAEFree(&dae_mem);
2666                 }
2667                 if (*neq > 0)
2668                 {
2669                     N_VDestroy_Serial(IDx);
2670                 }
2671                 if (*neq > 0)
2672                 {
2673                     N_VDestroy_Serial(yp);
2674                 }
2675                 if (*neq > 0)
2676                 {
2677                     N_VDestroy_Serial(yy);
2678                 }
2679                 if (ng > 0)
2680                 {
2681                     FREE(jroot);
2682                 }
2683                 if (ng > 0)
2684                 {
2685                     FREE(zcros);
2686                 }
2687                 if (nmod > 0)
2688                 {
2689                     FREE(Mode_save);
2690                 }
2691                 return;
2692             }
2693         }
2694         /*Jacobian_Flag=0; */
2695         if (AJacobian_block > 0) /* set by the block with A-Jac in flag-4 using Set_Jacobian_flag(1); */
2696         {
2697             Jn = *neq;
2698             Jnx = Blocks[AJacobian_block - 1].nx;
2699             Jno = Blocks[AJacobian_block - 1].nout;
2700             Jni = Blocks[AJacobian_block - 1].nin;
2701         }
2702         else
2703         {
2704             Jn = *neq;
2705             Jnx = 0;
2706             Jno = 0;
2707             Jni = 0;
2708         }
2709         Jactaille = 3 * Jn + (Jn + Jni) * (Jn + Jno) + Jnx * (Jni + 2 * Jn + Jno) + (Jn - Jnx) * (2 * (Jn - Jnx) + Jno + Jni) + 2 * Jni * Jno;
2710
2711         if ((data->rwork = (double *) MALLOC(Jactaille * sizeof(double))) == NULL)
2712         {
2713             if ( ng > 0 )
2714             {
2715                 FREE(data->gwork);
2716             }
2717             if (*neq > 0)
2718             {
2719                 N_VDestroy_Serial(data->ewt);
2720             }
2721             if (*neq > 0)
2722             {
2723                 FREE(data);
2724             }
2725             if (*neq > 0)
2726             {
2727                 DAEFree(&dae_mem);
2728             }
2729             if (*neq > 0)
2730             {
2731                 N_VDestroy_Serial(IDx);
2732             }
2733             if (*neq > 0)
2734             {
2735                 N_VDestroy_Serial(yp);
2736             }
2737             if (*neq > 0)
2738             {
2739                 N_VDestroy_Serial(yy);
2740             }
2741             if (ng > 0)
2742             {
2743                 FREE(jroot);
2744             }
2745             if (ng > 0)
2746             {
2747                 FREE(zcros);
2748             }
2749             if (nmod > 0)
2750             {
2751                 FREE(Mode_save);
2752             }
2753             *ierr = 10000;
2754             return;
2755         }
2756
2757         if (solver == IDA_BDF_Newton)
2758         {
2759             flag = IDADlsSetDenseJacFn(dae_mem, Jacobians);
2760         }
2761         if (check_flag(&flag, "IDADlsSetDenseJacFn", 1))
2762         {
2763             *ierr = 200 + (-flag);
2764             freeallx
2765             return;
2766         }
2767
2768         TJacque = (DlsMat) NewDenseMat(*neq, *neq);
2769
2770         flag = DAESetUserData(dae_mem, data);
2771         if (check_flag(&flag, "IDASetUserData", 1))
2772         {
2773             *ierr = 200 + (-flag);
2774             freeallx
2775             return;
2776         }
2777
2778         if (hmax > 0)
2779         {
2780             flag = DAESetMaxStep(dae_mem, (realtype) hmax);
2781             if (check_flag(&flag, "IDASetMaxStep", 1))
2782             {
2783                 *ierr = 200 + (-flag);
2784                 freeallx
2785                 return;
2786             }
2787         }
2788
2789         maxnj = 100; /* setting the maximum number of Jacobian evaluations during a Newton step */
2790         flag = DAESetMaxNumJacsIC(dae_mem, maxnj);
2791         if (check_flag(&flag, "IDASetMaxNumJacsIC", 1))
2792         {
2793             *ierr = 200 + (-flag);
2794             freeallx
2795             return;
2796         }
2797
2798         maxnit = 10; /* setting the maximum number of Newton iterations in any attempt to solve CIC */
2799         if (C2F(cmsolver).solver == 102)
2800         {
2801             maxnit = 15;    /* By default, the Krylov max iterations should be 15 */
2802         }
2803         flag = DAESetMaxNumItersIC(dae_mem, maxnit);
2804         if (check_flag(&flag, "IDASetMaxNumItersIC", 1))
2805         {
2806             *ierr = 200 + (-flag);
2807             freeallx
2808             return;
2809         }
2810
2811         /* setting the maximum number of steps in an integration interval */
2812         maxnh = 2000;
2813         flag = DAESetMaxNumSteps(dae_mem, maxnh);
2814         if (check_flag(&flag, "IDASetMaxNumSteps", 1))
2815         {
2816             *ierr = 200 + (-flag);
2817             freeallx
2818             return;
2819         }
2820
2821     } /* testing if neq>0 */
2822
2823     uround = 1.0;
2824     do
2825     {
2826         uround = uround * 0.5;
2827     }
2828     while ( 1.0 + uround != 1.0);
2829     uround = uround * 2.0;
2830     SQuround = sqrt(uround);
2831     /* Function Body */
2832
2833     C2F(coshlt).halt = 0;
2834     *ierr = 0;
2835     /*     hot = .false. */
2836     phase = 1;
2837     hot = 0;
2838
2839     /*      stuck=.false. */
2840     C2F(xscion)(&inxsci);
2841     /*     initialization */
2842     C2F(realtimeinit)(told, &C2F(rtfactor).scale);
2843     /*     ATOL and RTOL are scalars */
2844
2845     jj = 0;
2846     for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
2847     {
2848         if (Blocks[C2F(curblk).kfun - 1].ng >= 1)
2849         {
2850             zcros[jj] = C2F(curblk).kfun;
2851             ++jj;
2852         }
2853     }
2854     /*     . ng >= jj required */
2855     if (jj != ng)
2856     {
2857         zcros[jj] = -1;
2858     }
2859     /*     initialization (propagation of constant blocks outputs) */
2860     idoit(told);
2861     if (*ierr != 0)
2862     {
2863         freeallx;
2864         return;
2865     }
2866
2867     /*--discrete zero crossings----dzero--------------------*/
2868     if (ng > 0) /* storing ZC signs just after a solver call*/
2869     {
2870         zdoit(told, x, x, g);
2871         if (*ierr != 0)
2872         {
2873             freeallx;
2874             return;
2875         }
2876         for (jj = 0; jj < ng; ++jj)
2877             if (g[jj] >= 0)
2878             {
2879                 jroot[jj] = 5;
2880             }
2881             else
2882             {
2883                 jroot[jj] = -5;
2884             }
2885     }
2886     /*     main loop on time */
2887     while (*told < *tf)
2888     {
2889         while (ismenu()) //** if the user has done something, do the actions
2890         {
2891             int ierr2 = 0;
2892             SeqSync = GetCommand(CommandToUnstack); //** get to the action
2893             CommandLength = (int)strlen(CommandToUnstack);
2894             syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
2895         }
2896         if (C2F(coshlt).halt != 0)
2897         {
2898             if (C2F(coshlt).halt == 2)
2899             {
2900                 *told = *tf;    /* end simulation */
2901             }
2902             C2F(coshlt).halt = 0;
2903             freeallx;
2904             return;
2905         }
2906         if (*pointi == 0)
2907         {
2908             t = *tf;
2909         }
2910         else
2911         {
2912             t = tevts[*pointi];
2913         }
2914         if (fabs(t - *told) < ttol)
2915         {
2916             t = *told;
2917             /*     update output part */
2918         }
2919         if (*told > t)
2920         {
2921             /*     !  scheduling problem */
2922             *ierr = 1;
2923             freeallx;
2924             return;
2925         }
2926         if (*told != t)
2927         {
2928             if (xptr[nblk + 1] == 1)
2929             {
2930                 /*     .     no continuous state */
2931                 if (*told + deltat + ttol > t)
2932                 {
2933                     *told = t;
2934                 }
2935                 else
2936                 {
2937                     *told += deltat;
2938                 }
2939                 /*     .     update outputs of 'c' type blocks with no continuous state */
2940                 if (*told >= *tf)
2941                 {
2942                     /*     .     we are at the end, update continuous part before leaving */
2943                     cdoit(told);
2944                     freeallx;
2945                     return;
2946                 }
2947             }
2948             else
2949             {
2950                 rhotmp = *tf + ttol;
2951                 if (*pointi != 0)
2952                 {
2953                     kpo = *pointi;
2954 L20:
2955                     if (critev[kpo] == 1)
2956                     {
2957                         rhotmp = tevts[kpo];
2958                         goto L30;
2959                     }
2960                     kpo = evtspt[kpo];
2961                     if (kpo != 0)
2962                     {
2963                         goto L20;
2964                     }
2965 L30:
2966                     if (rhotmp < tstop)
2967                     {
2968                         hot = 0;/* Cold-restart the solver if the new TSTOP isn't beyong the previous one*/
2969                     }
2970                 }
2971                 tstop = rhotmp;
2972                 t = Min(*told + deltat, Min(t, *tf + ttol));
2973
2974                 if (hot == 0)  /* CIC calculation when hot==0 */
2975                 {
2976
2977                     /* Setting the stop time*/
2978                     flag = DAESetStopTime(dae_mem, (realtype)tstop);
2979                     if (check_flag(&flag, "IDASetStopTime", 1))
2980                     {
2981                         *ierr = 200 + (-flag);
2982                         freeallx;
2983                         return;
2984                     }
2985
2986                     if (ng > 0 && nmod > 0)
2987                     {
2988                         phase = 1;
2989                         zdoit(told, x, xd, g);
2990                         if (*ierr != 0)
2991                         {
2992                             freeallx;
2993                             return;
2994                         }
2995                     }
2996
2997                     /*----------ID setting/checking------------*/
2998                     N_VConst(ONE, IDx); /* Initialize id to 1's. */
2999                     scicos_xproperty = NV_DATA_S(IDx);
3000                     reinitdoit(told);
3001                     if (*ierr > 0)
3002                     {
3003                         freeallx;
3004                         return;
3005                     }
3006                     for (jj = 0; jj < *neq; jj++)
3007                     {
3008                         if (xprop[jj] ==  1)
3009                         {
3010                             scicos_xproperty[jj] = ONE;
3011                         }
3012                         if (xprop[jj] == -1)
3013                         {
3014                             scicos_xproperty[jj] = ZERO;
3015                         }
3016                     }
3017                     /* CI=0.0;CJ=100.0; // for functions Get_Jacobian_ci and Get_Jacobian_cj
3018                     Jacobians(*neq, (realtype) (*told), yy, yp, bidon, (realtype) CJ, data, TJacque, tempv1, tempv2, tempv3);
3019                     for (jj=0;jj<*neq;jj++){
3020                     Jacque_col=DENSE_COL(TJacque,jj);
3021                     CI=ZERO;
3022                     for (kk=0;kk<*neq;kk++){
3023                     if ((Jacque_col[kk]-Jacque_col[kk]!=0)) {
3024                     CI=-ONE;
3025                     break;
3026                     }else{
3027                     if (Jacque_col[kk]!=0){
3028                     CI=ONE;
3029                     break;
3030                     }
3031                     }
3032                     }
3033                     if (CI>=ZERO){  scicos_xproperty[jj]=CI;}else{fprintf(stderr,"\nWarinng! Xproperties are not match for i=%d!",jj);}
3034                     } */
3035                     /* printf("\n"); for(jj=0;jj<*neq;jj++) { printf("x%d=%g ",jj,scicos_xproperty[jj]); }*/
3036
3037                     flag = DAESetId(dae_mem, IDx);
3038                     if (check_flag(&flag, "IDASetId", 1))
3039                     {
3040                         *ierr = 200 + (-flag);
3041                         freeallx
3042                         return;
3043                     }
3044                     // CI=1.0;  // for function Get_Jacobian_ci
3045                     /*--------------------------------------------*/
3046                     // maxnj=100; /* setting the maximum number of Jacobian evaluation during a Newton step */
3047                     // flag=DAESetMaxNumJacsIC(dae_mem, maxnj);
3048                     // if (check_flag(&flag, "IDASetMaxNumJacsIC", 1)) {
3049                     //   *ierr=200+(-flag);
3050                     //   freeallx;
3051                     //   return;
3052                     // };
3053                     // flag=DAESetLineSearchOffIC(dae_mem,FALSE);  /* (def=false)  */
3054                     // if (check_flag(&flag, "IDASetLineSearchOffIC", 1)) {
3055                     //   *ierr=200+(-flag);
3056                     //   freeallx;
3057                     //   return;
3058                     // };
3059                     // flag=DAESetMaxNumItersIC(dae_mem, 10);/* (def=10) setting the maximum number of Newton iterations in any attempt to solve CIC */
3060                     // if (check_flag(&flag, "IDASetMaxNumItersIC", 1)) {
3061                     //   *ierr=200+(-flag);
3062                     //   freeallx;
3063                     //   return;
3064                     // };
3065
3066                     N_iters = 4 + nmod * 4;
3067                     for (j = 0; j <= N_iters; j++)
3068                     {
3069                         /* counter to reevaluate the
3070                                                                         modes in  mode->CIC->mode->CIC-> loop
3071                                                                         do it once in the absence of mode (nmod=0) */
3072                         /* updating the modes through Flag==9, Phase==1 */
3073
3074                         /* Serge Steer 29/06/2009 */
3075                         while (ismenu()) //** if the user has done something, do the actions
3076                         {
3077                             int ierr2 = 0;
3078                             SeqSync = GetCommand(CommandToUnstack); //** get to the action
3079                             CommandLength = (int)strlen(CommandToUnstack);
3080                             syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
3081                         }
3082                         if (C2F(coshlt).halt != 0)
3083                         {
3084                             if (C2F(coshlt).halt == 2)
3085                             {
3086                                 *told = *tf;    /* end simulation */
3087                             }
3088                             C2F(coshlt).halt = 0;
3089                             freeallx;
3090                             return;
3091                         }
3092
3093                         /* yy->PH */
3094                         flag = DAEReInit(dae_mem, (realtype)(*told), yy, yp);
3095                         if (check_flag(&flag, "IDAReInit", 1))
3096                         {
3097                             *ierr = 200 + (-flag);
3098                             freeallx;
3099                             return;
3100                         }
3101
3102                         phase = 2; /* IDACalcIC: PHI-> yy0: if (ok) yy0_cic-> PHI*/
3103                         if (C2F(cmsolver).solver == 100)
3104                         {
3105                             copy_IDA_mem->ida_kk = 1;
3106                         }
3107                         // the initial conditons y0 and yp0 do not satisfy the DAE
3108                         flagr = DAECalcIC(dae_mem, DAE_YA_YDP_INIT, (realtype)(t));
3109                         phase = 1;
3110                         flag = DAEGetConsistentIC(dae_mem, yy, yp); /* PHI->YY */
3111                         if (*ierr > 5)    /* *ierr>5 => singularity in block */
3112                         {
3113                             freeallx;
3114                             return;
3115                         }
3116
3117                         if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3118                         {
3119                             if (flagr >= 0)
3120                             {
3121                                 sciprint(_("**** SUNDIALS.IDA successfully initialized *****\n") );
3122                             }
3123                             else
3124                             {
3125                                 sciprint(_("**** SUNDIALS.IDA failed to initialize ->try again *****\n") );
3126                             }
3127                         }
3128                         /*-------------------------------------*/
3129                         /* saving the previous modes*/
3130                         for (jj = 0; jj < nmod; ++jj)
3131                         {
3132                             Mode_save[jj] = mod[jj];
3133                         }
3134                         if (ng > 0 && nmod > 0)
3135                         {
3136                             phase = 1;
3137                             zdoit(told, x, xd, g);
3138                             if (*ierr != 0)
3139                             {
3140                                 freeallx;
3141                                 return;
3142                             }
3143                         }
3144                         /*------------------------------------*/
3145                         Mode_change = 0;
3146                         for (jj = 0; jj < nmod; ++jj)
3147                         {
3148                             if (Mode_save[jj] != mod[jj])
3149                             {
3150                                 Mode_change = 1;
3151                                 break;
3152                             }
3153                         }
3154                         if (Mode_change == 0)
3155                         {
3156                             if (flagr >= 0 )
3157                             {
3158                                 break;  /*   if (flagr>=0) break;  else{ *ierr=200+(-flagr); freeallx;  return; }*/
3159                             }
3160                             else if (j >= (int)( N_iters / 2))
3161                             {
3162                                 /* DAESetMaxNumStepsIC(mem,10); */     /* maxnh (def=5) */
3163                                 DAESetMaxNumJacsIC(dae_mem, 10);      /* maxnj 100 (def=4)*/
3164                                 /* DAESetMaxNumItersIC(mem,100000); */ /* maxnit in IDANewtonIC (def=10) */
3165                                 DAESetLineSearchOffIC(dae_mem, TRUE); /* (def=false)  */
3166                                 /* DAESetNonlinConvCoefIC(mem,1.01);*/ /* (def=0.01-0.33*/
3167                                 flag = DAESetMaxNumItersIC(dae_mem, 1000);
3168                                 if (check_flag(&flag, "IDASetMaxNumItersIC", 1))
3169                                 {
3170                                     *ierr = 200 + (-flag);
3171                                     freeallx;
3172                                     return;
3173                                 };
3174                             }
3175                         }
3176                     }/* mode-CIC  counter*/
3177                     if (Mode_change == 1)
3178                     {
3179                         /* In this case, we try again by relaxing all modes and calling DAECalcIC again
3180                         /Masoud */
3181                         phase = 1;
3182                         if (C2F(cmsolver).solver == 100)
3183                         {
3184                             copy_IDA_mem->ida_kk = 1;
3185                         }
3186                         flagr = DAECalcIC(dae_mem, DAE_YA_YDP_INIT, (realtype)(t));
3187                         phase = 1;
3188                         flag = DAEGetConsistentIC(dae_mem, yy, yp); /* PHI->YY */
3189                         if ((flagr < 0) || (*ierr > 5)) /* *ierr>5 => singularity in block */
3190                         {
3191                             *ierr = 23;
3192                             freeallx;
3193                             return;
3194                         }
3195                     }
3196                     /*-----If flagr<0 the initialization solver has not converged-----*/
3197                     if (flagr < 0)
3198                     {
3199                         *ierr = 237;
3200                         freeallx;
3201                         return;
3202                     }
3203
3204                 } /* CIC calculation when hot==0 */
3205
3206                 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3207                 {
3208                     sciprint(_("****daskr from: %f to %f hot= %d  \n"), *told, t, hot);
3209                 }
3210
3211                 /*--discrete zero crossings----dzero--------------------*/
3212                 /*--check for Dzeros after Mode settings or ddoit()----*/
3213                 Discrete_Jump = 0;
3214                 if (ng > 0 && hot == 0)
3215                 {
3216                     zdoit(told, x, xd, g);
3217                     if (*ierr != 0)
3218                     {
3219                         freeallx;
3220                         return;
3221                     }
3222                     for (jj = 0; jj < ng; ++jj)
3223                     {
3224                         if ((g[jj] >= 0.0) && ( jroot[jj] == -5))
3225                         {
3226                             Discrete_Jump = 1;
3227                             jroot[jj] = 1;
3228                         }
3229                         else if ((g[jj] < 0.0) && ( jroot[jj] == 5))
3230                         {
3231                             Discrete_Jump = 1;
3232                             jroot[jj] = -1;
3233                         }
3234                         else
3235                         {
3236                             jroot[jj] = 0;
3237                         }
3238                     }
3239                 }
3240
3241                 /*--discrete zero crossings----dzero--------------------*/
3242                 if (Discrete_Jump == 0) /* if there was a dzero, its event should be activated*/
3243                 {
3244                     phase = 2;
3245                     flagr = DAESolve(dae_mem, t, told, yy, yp, DAE_NORMAL);
3246                     phase = 1;
3247                     if (*ierr != 0)
3248                     {
3249                         freeallx;
3250                         return;
3251                     }
3252                 }
3253                 else
3254                 {
3255                     flagr = IDA_ROOT_RETURN; /* in order to handle discrete jumps */
3256                 }
3257                 if (flagr >= 0)
3258                 {
3259                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3260                     {
3261                         sciprint(_("****SUNDIALS.Ida reached: %f\n"), *told);
3262                     }
3263                     hot = 1;
3264                     cnt = 0;
3265                 }
3266                 else if ( flagr == IDA_TOO_MUCH_WORK ||  flagr == IDA_CONV_FAIL || flagr == IDA_ERR_FAIL)
3267                 {
3268                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3269                     {
3270                         sciprint(_("**** SUNDIALS.Ida: too much work at time=%g (stiff region, change RTOL and ATOL)\n"), *told);
3271                     }
3272                     hot = 0;
3273                     cnt++;
3274                     if (cnt > 5)
3275                     {
3276                         *ierr = 200 + (-flagr);
3277                         freeallx;
3278                         return;
3279                     }
3280                 }
3281                 else
3282                 {
3283                     if (flagr < 0)
3284                     {
3285                         *ierr = 200 + (-flagr);    /* raising errors due to internal errors, otherwise error due to flagr*/
3286                     }
3287                     freeallx;
3288                     return;
3289                 }
3290
3291                 /*     update outputs of 'c' type  blocks if we are at the end*/
3292                 if (*told >= *tf)
3293                 {
3294                     cdoit(told);
3295                     freeallx;
3296                     return;
3297                 }
3298
3299                 if (flagr == IDA_ZERO_DETACH_RETURN)
3300                 {
3301                     hot = 0;
3302                 }; /* new feature of sundials, detects unmasking */
3303                 if (flagr == IDA_ROOT_RETURN)
3304                 {
3305                     /*     .        at least one root has been found */
3306                     hot = 0;
3307                     if (Discrete_Jump == 0)
3308                     {
3309                         flagr = DAEGetRootInfo(dae_mem, jroot);
3310                         if (check_flag(&flagr, "IDAGetRootInfo", 1))
3311                         {
3312                             *ierr = 200 + (-flagr);
3313                             freeallx;
3314                             return;
3315                         }
3316                     }
3317
3318                     if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3319                     {
3320                         sciprint(_("root found at t=: %f\n"), *told);
3321                     }
3322                     /*     .        update outputs affecting ztyp blocks  ONLY FOR OLD BLOCKS*/
3323                     zdoit(told, x, xd, g);
3324                     if (*ierr != 0)
3325                     {
3326                         freeallx;
3327                         return;
3328                     }
3329                     for (jj = 0; jj < ng; ++jj)
3330                     {
3331                         C2F(curblk).kfun = zcros[jj];
3332                         if (C2F(curblk).kfun == -1)
3333                         {
3334                             break;
3335                         }
3336                         kev = 0;
3337                         for (j = zcptr[C2F(curblk).kfun] - 1 ;
3338                                 j < zcptr[C2F(curblk).kfun + 1] - 1 ; ++j)
3339                         {
3340                             if (jroot[j] != 0)
3341                             {
3342                                 kev = 1;
3343                                 break;
3344                             }
3345                         }
3346                         if (kev != 0)
3347                         {
3348                             Blocks[C2F(curblk).kfun - 1].jroot = &jroot[zcptr[C2F(curblk).kfun] - 1];
3349                             if (funtyp[C2F(curblk).kfun] > 0)
3350                             {
3351                                 if (Blocks[C2F(curblk).kfun - 1].nevout > 0)
3352                                 {
3353                                     flag__ = 3;
3354                                     if (Blocks[C2F(curblk).kfun - 1].nx > 0)
3355                                     {
3356                                         Blocks[C2F(curblk).kfun - 1].x  = &x[xptr[C2F(curblk).kfun] - 1];
3357                                         Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
3358                                     }
3359                                     /*     call corresponding block to determine output event (kev) */
3360                                     Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
3361                                     /*callf(told, xd, x, x,g,&flag__);*/
3362                                     callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
3363                                     if (flag__ < 0)
3364                                     {
3365                                         *ierr = 5 - flag__;
3366                                         freeallx;
3367                                         return;
3368                                     }
3369                                     /*     update event agenda */
3370                                     for (k = 0; k < Blocks[C2F(curblk).kfun - 1].nevout; ++k)
3371                                     {
3372                                         if (Blocks[C2F(curblk).kfun - 1].evout[k] >= 0)
3373                                         {
3374                                             i3 = k + clkptr[C2F(curblk).kfun] ;
3375                                             addevs(Blocks[C2F(curblk).kfun - 1].evout[k] + (*told), &i3, &ierr1);
3376                                             if (ierr1 != 0)
3377                                             {
3378                                                 /*     .                       nevts too small */
3379                                                 *ierr = 3;
3380                                                 freeallx;
3381                                                 return;
3382                                             }
3383                                         }
3384                                     }
3385                                 }
3386                                 /* update state */
3387                                 if ((Blocks[C2F(curblk).kfun - 1].nx > 0) || (*Blocks[C2F(curblk).kfun - 1].work != NULL) )
3388                                 {
3389                                     /* call corresponding block to update state */
3390                                     flag__ = 2;
3391                                     if (Blocks[C2F(curblk).kfun - 1].nx > 0)
3392                                     {
3393                                         Blocks[C2F(curblk).kfun - 1].x  = &x[xptr[C2F(curblk).kfun] - 1];
3394                                         Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
3395                                     }
3396                                     Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
3397
3398                                     Blocks[C2F(curblk).kfun - 1].xprop = &xprop[-1 + xptr[C2F(curblk).kfun]];
3399                                     /*callf(told, xd, x, x,g,&flag__);*/
3400                                     callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
3401
3402                                     if (flag__ < 0)
3403                                     {
3404                                         *ierr = 5 - flag__;
3405                                         freeallx;
3406                                         return;
3407                                     }
3408                                     for (j = 0; j < *neq; j++) /* Adjust xprop for IDx */
3409                                     {
3410                                         if (xprop[j] ==  1)
3411                                         {
3412                                             scicos_xproperty[j] = ONE;
3413                                         }
3414                                         if (xprop[j] == -1)
3415                                         {
3416                                             scicos_xproperty[j] = ZERO;
3417                                         }
3418                                     }
3419                                 }
3420                             }
3421                         }
3422                     }
3423                 }
3424                 /* Serge Steer 29/06/2009 */
3425                 while (ismenu()) //** if the user has done something, do the actions
3426                 {
3427                     int ierr2 = 0;
3428                     SeqSync = GetCommand(CommandToUnstack); //** get to the action
3429                     CommandLength = (int)strlen(CommandToUnstack);
3430                     syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
3431                 }
3432
3433                 if (C2F(coshlt).halt != 0)
3434                 {
3435                     if (C2F(coshlt).halt == 2)
3436                     {
3437                         *told = *tf;    /* end simulation */
3438                     }
3439                     C2F(coshlt).halt = 0;
3440                     freeallx;
3441                     return;
3442                 }
3443                 /* if(*pointi!=0){
3444                 t=tevts[*pointi];
3445                 if(*told<t-ttol){
3446                 cdoit(told);
3447                 goto L15;
3448                 }
3449                 }else{
3450                 if(*told<*tf){
3451                 cdoit(told);
3452                 goto L15;
3453                 }
3454                 }*/
3455
3456                 /*--discrete zero crossings----dzero--------------------*/
3457                 if (ng > 0) /* storing ZC signs just after a ddaskr call*/
3458                 {
3459                     zdoit(told, x, xd, g);
3460                     if (*ierr != 0)
3461                     {
3462                         freeallx;
3463                         return;
3464                     }
3465                     for (jj = 0; jj < ng; ++jj)
3466                     {
3467                         if (g[jj] >= 0)
3468                         {
3469                             jroot[jj] = 5;
3470                         }
3471                         else
3472                         {
3473                             jroot[jj] = -5;
3474                         }
3475                     }
3476                 }
3477                 /*--discrete zero crossings----dzero--------------------*/
3478             }
3479             C2F(realtime)(told);
3480         }
3481         else
3482         {
3483             /*     .  t==told */
3484             if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3485             {
3486                 sciprint(_("Event: %d activated at t=%f\n"), *pointi, *told);
3487             }
3488
3489             ddoit(told);
3490             if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3491             {
3492                 sciprint(_("End of activation\n"));
3493             }
3494             if (*ierr != 0)
3495             {
3496                 freeallx;
3497                 return;
3498             }
3499         }
3500         /*     end of main loop on time */
3501     }
3502     freeallx;
3503 } /* cossimdaskr_ */
3504 /*--------------------------------------------------------------------------*/
3505 /* Subroutine cosend */
3506 static void cosend(double *told)
3507 {
3508     /* Local variables */
3509     static scicos_flag flag__ = 0;
3510
3511     static int kfune = 0;
3512
3513     /* Function Body */
3514     *ierr = 0;
3515     /*     loop on blocks */
3516     for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
3517     {
3518         flag__ = 5;
3519         Blocks[C2F(curblk).kfun - 1].nevprt = 0;
3520         if (funtyp[C2F(curblk).kfun] >= 0)
3521         {
3522             if (Blocks[C2F(curblk).kfun - 1].nx > 0)
3523             {
3524                 Blocks[C2F(curblk).kfun - 1].x  = &x[xptr[C2F(curblk).kfun] - 1];
3525                 Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
3526             }
3527             /*callf(told, xd, x, x,x,&flag__);*/
3528             callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
3529             if (flag__ < 0 && *ierr == 0)
3530             {
3531                 *ierr = 5 - flag__;
3532                 kfune = C2F(curblk).kfun;
3533             }
3534         }
3535     }
3536     if (*ierr != 0)
3537     {
3538         C2F(curblk).kfun = kfune;
3539         return;
3540     }
3541 } /* cosend_ */
3542 /*--------------------------------------------------------------------------*/
3543 /* callf */
3544 void callf(double *t, scicos_block *block, scicos_flag *flag)
3545 {
3546     double* args[SZ_SIZE];
3547     int sz[SZ_SIZE];
3548     double intabl[TB_SIZE];
3549     double outabl[TB_SIZE];
3550
3551     int ii = 0, in = 0, out = 0, ki = 0, ko = 0, no = 0, ni = 0, k = 0, j = 0;
3552     int szi = 0, flagi = 0;
3553     double *ptr_d = NULL;
3554
3555     /* function pointers type def */
3556     voidf loc ;
3557     ScicosF0 loc0;
3558     ScicosF loc1;
3559     /*  ScicosFm1 loc3;*/
3560     ScicosF2 loc2;
3561     ScicosF2z loc2z;
3562     ScicosFi loci1;
3563     ScicosFi2 loci2;
3564     ScicosFi2z loci2z;
3565     ScicosF4 loc4;
3566
3567     int solver = C2F(cmsolver).solver;
3568     int cosd   = C2F(cosdebug).cosd;
3569     /*int kf     = C2F(curblk).kfun;*/
3570     scicos_time     = *t;
3571     block_error     = (int*) flag;
3572
3573     /* debug block is never called */
3574     /*if (kf==(debug_block+1)) return;*/
3575     if (block->type == 99)
3576     {
3577         return;
3578     }
3579
3580     /* flag 7 implicit initialization */
3581     flagi = (int) * flag;
3582     /* change flag to zero if flagi==7 for explicit block */
3583     if (flagi == 7 && block->type < 10000)
3584     {
3585         *flag = 0;
3586     }
3587
3588     /* display information for debugging mode */
3589     if (cosd > 1)
3590     {
3591         if (cosd != 3)
3592         {
3593             sciprint(_("block %d is called "), C2F(curblk).kfun);
3594             sciprint(_("with flag %d "), *flag);
3595             sciprint(_("at time %f \n"), *t);
3596         }
3597         if (debug_block > -1)
3598         {
3599             if (cosd != 3)
3600             {
3601                 sciprint(_("Entering the block \n"));
3602             }
3603             call_debug_scicos(block, flag, flagi, debug_block);
3604             if (*flag < 0)
3605             {
3606                 return;    /* error in debug block */
3607             }
3608         }
3609     }
3610
3611     C2F(scsptr).ptr = block->scsptr;
3612
3613     /* get pointer of the function */
3614     loc = block->funpt;
3615
3616     /* continuous state */
3617     if ((solver == IDA_BDF_Newton || solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes) && block->type < 10000 && *flag == 0)
3618     {
3619         ptr_d = block->xd;
3620         block->xd  = block->res;
3621     }
3622
3623     /* switch loop */
3624     //sciprint("callf type=%d flag=%d\n",block->type,flagi);
3625     switch (block->type)
3626     {
3627             /*******************/
3628             /* function type 0 */
3629             /*******************/
3630         case 0 :
3631         {
3632             /* This is for compatibility */
3633             /* jroot is returned in g for old type */
3634             if (block->nevprt < 0)
3635             {
3636                 for (j = 0; j < block->ng; ++j)
3637                 {
3638                     block->g[j] = (double)block->jroot[j];
3639                 }
3640             }
3641
3642             /* concatenated entries and concatened outputs */
3643             /* catenate inputs if necessary */
3644             ni = 0;
3645             if (block->nin > 1)
3646             {
3647                 ki = 0;
3648                 for (in = 0; in < block->nin; in++)
3649                 {
3650                     szi = block->insz[in] * block->insz[in + block->nin];
3651                     for (ii = 0; ii < szi; ii++)
3652                     {
3653                         intabl[ki++] = *((double *)(block->inptr[in]) + ii);
3654                     }
3655                     ni = ni + szi;
3656                 }
3657                 args[0] = &(intabl[0]);
3658             }
3659             else
3660             {
3661                 if (block->nin == 0)
3662                 {
3663                     args[0] = NULL;
3664                 }
3665                 else
3666                 {
3667                     args[0] = (double *)(block->inptr[0]);
3668                     ni = block->insz[0] * block->insz[1];
3669                 }
3670             }
3671
3672             /* catenate outputs if necessary */
3673             no = 0;
3674             if (block->nout > 1)
3675             {
3676                 ko = 0;
3677                 for (out = 0; out < block->nout; out++)
3678                 {
3679                     szi = block->outsz[out] * block->outsz[out + block->nout];
3680                     for (ii = 0; ii < szi; ii++)
3681                     {
3682                         outabl[ko++] = *((double *)(block->outptr[out]) + ii);
3683                     }
3684                     no = no + szi;
3685                 }
3686                 args[1] = &(outabl[0]);
3687             }
3688             else
3689             {
3690                 if (block->nout == 0)
3691                 {
3692                     args[1] = NULL;
3693                 }
3694                 else
3695                 {
3696                     args[1] = (double *)(block->outptr[0]);
3697                     no = block->outsz[0] * block->outsz[1];
3698                 }
3699             }
3700
3701             loc0 = (ScicosF0) loc;
3702
3703             (*loc0)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3704                     block->z, &block->nz,
3705                     block->evout, &block->nevout, block->rpar, &block->nrpar,
3706                     block->ipar, &block->nipar, (double *)args[0], &ni,
3707                     (double *)args[1], &no);
3708
3709             /* split output vector on each port if necessary */
3710             if (block->nout > 1)
3711             {
3712                 ko = 0;
3713                 for (out = 0; out < block->nout; out++)
3714                 {
3715                     szi = block->outsz[out] * block->outsz[out + block->nout];
3716                     for (ii = 0; ii < szi; ii++)
3717                     {
3718                         *((double *)(block->outptr[out]) + ii) = outabl[ko++];
3719                     }
3720                 }
3721             }
3722
3723             /* adjust values of output register */
3724             for (in = 0; in < block->nevout; ++in)
3725             {
3726                 block->evout[in] = block->evout[in] - *t;
3727             }
3728
3729             break;
3730         }
3731
3732         /*******************/
3733         /* function type 1 */
3734         /*******************/
3735         case 1 :
3736         {
3737             /* This is for compatibility */
3738             /* jroot is returned in g for old type */
3739             if (block->nevprt < 0)
3740             {
3741                 for (j = 0; j < block->ng; ++j)
3742                 {
3743                     block->g[j] = (double)block->jroot[j];
3744                 }
3745             }
3746
3747             /* one entry for each input or output */
3748             for (in = 0 ; in < block->nin ; in++)
3749             {
3750                 args[in] = block->inptr[in];
3751                 sz[in] = block->insz[in];
3752             }
3753             for (out = 0; out < block->nout; out++)
3754             {
3755                 args[in + out] = block->outptr[out];
3756                 sz[in + out] = block->outsz[out];
3757             }
3758             /* with zero crossing */
3759             if (block->ztyp > 0)
3760             {
3761                 args[block->nin + block->nout] = block->g;
3762                 sz[block->nin + block->nout] = block->ng;
3763             }
3764
3765             loc1 = (ScicosF) loc;
3766
3767             (*loc1)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3768                     block->z, &block->nz,
3769                     block->evout, &block->nevout, block->rpar, &block->nrpar,
3770                     block->ipar, &block->nipar,
3771                     (double *)args[0], &sz[0],
3772                     (double *)args[1], &sz[1], (double *)args[2], &sz[2],
3773                     (double *)args[3], &sz[3], (double *)args[4], &sz[4],
3774                     (double *)args[5], &sz[5], (double *)args[6], &sz[6],
3775                     (double *)args[7], &sz[7], (double *)args[8], &sz[8],
3776                     (double *)args[9], &sz[9], (double *)args[10], &sz[10],
3777                     (double *)args[11], &sz[11], (double *)args[12], &sz[12],
3778                     (double *)args[13], &sz[13], (double *)args[14], &sz[14],
3779                     (double *)args[15], &sz[15], (double *)args[16], &sz[16],
3780                     (double *)args[17], &sz[17]);
3781
3782             /* adjust values of output register */
3783             for (in = 0; in < block->nevout; ++in)
3784             {
3785                 block->evout[in] = block->evout[in] - *t;
3786             }
3787
3788             break;
3789         }
3790
3791         /*******************/
3792         /* function type 2 */
3793         /*******************/
3794         case 2 :
3795         {
3796             /* This is for compatibility */
3797             /* jroot is returned in g for old type */
3798             if (block->nevprt < 0)
3799             {
3800                 for (j = 0; j < block->ng; ++j)
3801                 {
3802                     block->g[j] = (double)block->jroot[j];
3803                 }
3804             }
3805
3806             /* no zero crossing */
3807             if (block->ztyp == 0)
3808             {
3809                 loc2 = (ScicosF2) loc;
3810                 (*loc2)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3811                         block->z, &block->nz,
3812                         block->evout, &block->nevout, block->rpar, &block->nrpar,
3813                         block->ipar, &block->nipar, (double **)block->inptr,
3814                         block->insz, &block->nin,
3815                         (double **)block->outptr, block->outsz, &block->nout);
3816             }
3817             /* with zero crossing */
3818             else
3819             {
3820                 loc2z = (ScicosF2z) loc;
3821                 (*loc2z)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3822                          block->z, &block->nz,
3823                          block->evout, &block->nevout, block->rpar, &block->nrpar,
3824                          block->ipar, &block->nipar, (double **)block->inptr,
3825                          block->insz, &block->nin,
3826                          (double **)block->outptr, block->outsz, &block->nout,
3827                          block->g, &block->ng);
3828             }
3829
3830             /* adjust values of output register */
3831             for (in = 0; in < block->nevout; ++in)
3832             {
3833                 block->evout[in] = block->evout[in] - *t;
3834             }
3835
3836             break;
3837         }
3838
3839         /*******************/
3840         /* function type 4 */
3841         /*******************/
3842         case 4 :
3843         {
3844             /* get pointer of the function type 4*/
3845             loc4 = (ScicosF4) loc;
3846
3847             (*loc4)(block, *flag);
3848
3849             break;
3850         }
3851
3852         /***********************/
3853         /* function type 10001 */
3854         /***********************/
3855         case 10001 :
3856         {
3857             /* This is for compatibility */
3858             /* jroot is returned in g for old type */
3859             if (block->nevprt < 0)
3860             {
3861                 for (j = 0; j < block->ng; ++j)
3862                 {
3863                     block->g[j] = (double)block->jroot[j];
3864                 }
3865             }
3866
3867             /* implicit block one entry for each input or output */
3868             for (in = 0 ; in < block->nin ; in++)
3869             {
3870                 args[in] = block->inptr[in];
3871                 sz[in] = block->insz[in];
3872             }
3873             for (out = 0; out < block->nout; out++)
3874             {
3875                 args[in + out] = block->outptr[out];
3876                 sz[in + out] = block->outsz[out];
3877             }
3878             /* with zero crossing */
3879             if (block->ztyp > 0)
3880             {
3881                 args[block->nin + block->nout] = block->g;
3882                 sz[block->nin + block->nout] = block->ng;
3883             }
3884
3885             loci1 = (ScicosFi) loc;
3886             (*loci1)(flag, &block->nevprt, t, block->res, block->xd, block->x,
3887                      &block->nx, block->z, &block->nz,
3888                      block->evout, &block->nevout, block->rpar, &block->nrpar,
3889                      block->ipar, &block->nipar,
3890                      (double *)args[0], &sz[0],
3891                      (double *)args[1], &sz[1], (double *)args[2], &sz[2],
3892                      (double *)args[3], &sz[3], (double *)args[4], &sz[4],
3893                      (double *)args[5], &sz[5], (double *)args[6], &sz[6],
3894                      (double *)args[7], &sz[7], (double *)args[8], &sz[8],
3895                      (double *)args[9], &sz[9], (double *)args[10], &sz[10],
3896                      (double *)args[11], &sz[11], (double *)args[12], &sz[12],
3897                      (double *)args[13], &sz[13], (double *)args[14], &sz[14],
3898                      (double *)args[15], &sz[15], (double *)args[16], &sz[16],
3899                      (double *)args[17], &sz[17]);
3900
3901             /* adjust values of output register */
3902             for (in = 0; in < block->nevout; ++in)
3903             {
3904                 block->evout[in] = block->evout[in] - *t;
3905             }
3906
3907             break;
3908         }
3909
3910         /***********************/
3911         /* function type 10002 */
3912         /***********************/
3913         case 10002 :
3914         {
3915             /* This is for compatibility */
3916             /* jroot is returned in g for old type */
3917             if (block->nevprt < 0)
3918             {
3919                 for (j = 0; j < block->ng; ++j)
3920                 {
3921                     block->g[j] = (double)block->jroot[j];
3922                 }
3923             }
3924
3925             /* implicit block, inputs and outputs given by a table of pointers */
3926             /* no zero crossing */
3927             if (block->ztyp == 0)
3928             {
3929                 loci2 = (ScicosFi2) loc;
3930                 (*loci2)(flag, &block->nevprt, t, block->res,
3931                          block->xd, block->x, &block->nx,
3932                          block->z, &block->nz,
3933                          block->evout, &block->nevout, block->rpar, &block->nrpar,
3934                          block->ipar, &block->nipar, (double **)block->inptr,
3935                          block->insz, &block->nin,
3936                          (double **)block->outptr, block->outsz, &block->nout);
3937             }
3938             /* with zero crossing */
3939             else
3940             {
3941                 loci2z = (ScicosFi2z) loc;
3942                 (*loci2z)(flag, &block->nevprt, t, block->res,
3943                           block->xd, block->x, &block->nx,
3944                           block->z, &block->nz,
3945                           block->evout, &block->nevout, block->rpar, &block->nrpar,
3946                           block->ipar, &block->nipar,
3947                           (double **)block->inptr, block->insz, &block->nin,
3948                           (double **)block->outptr, block->outsz, &block->nout,
3949                           block->g, &block->ng);
3950             }
3951
3952             /* adjust values of output register */
3953             for (in = 0; in < block->nevout; ++in)
3954             {
3955                 block->evout[in] = block->evout[in] - *t;
3956             }
3957
3958             break;
3959         }
3960
3961         /***********************/
3962         /* function type 10004 */
3963         /***********************/
3964         case 10004 :
3965         {
3966             /* get pointer of the function type 4*/
3967             loc4 = (ScicosF4) loc;
3968
3969             (*loc4)(block, *flag);
3970
3971             break;
3972         }
3973
3974         /***********/
3975         /* default */
3976         /***********/
3977         default :
3978         {
3979             sciprint(_("Undefined Function type\n"));
3980             *flag = -1000;
3981             return; /* exit */
3982         }
3983     }
3984     // sciprint("callf end  flag=%d\n",*flag);
3985     /* Implicit Solver & explicit block & flag==0 */
3986     /* adjust continuous state vector after call */
3987     if ((solver == IDA_BDF_Newton || solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes) && block->type < 10000 && *flag == 0)
3988     {
3989         block->xd  = ptr_d;
3990         if (flagi != 7)
3991         {
3992             for (k = 0; k < block->nx; k++)
3993             {
3994                 block->res[k] = block->res[k] - block->xd[k];
3995             }
3996         }
3997         else
3998         {
3999             for (k = 0; k < block->nx; k++)
4000             {
4001                 block->xd[k] = block->res[k];
4002             }
4003         }
4004     }
4005
4006     /* debug block */
4007     if (cosd > 1)
4008     {
4009         if (debug_block > -1)
4010         {
4011             if (*flag < 0)
4012             {
4013                 return;    /* error in block */
4014             }
4015             if (cosd != 3)
4016             {
4017                 sciprint(_("Leaving block %d \n"), C2F(curblk).kfun);
4018             }
4019             call_debug_scicos(block, flag, flagi, debug_block);
4020             /*call_debug_scicos(flag,kf,flagi,debug_block);*/
4021         }
4022     }
4023 } /* callf */
4024 /*--------------------------------------------------------------------------*/
4025 /* call_debug_scicos */
4026 static void call_debug_scicos(scicos_block *block, scicos_flag *flag, int flagi, int deb_blk)
4027 {
4028     voidf loc ;
4029     int solver = C2F(cmsolver).solver, k = 0;
4030     ScicosF4 loc4;
4031     double *ptr_d = NULL;
4032
4033     C2F(cosdebugcounter).counter = C2F(cosdebugcounter).counter + 1;
4034     C2F(scsptr).ptr = Blocks[deb_blk].scsptr;
4035
4036     loc  = Blocks[deb_blk].funpt; /* GLOBAL */
4037     loc4 = (ScicosF4) loc;
4038
4039     /* continuous state */
4040     if ((solver == IDA_BDF_Newton || solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes) && block->type < 10000 && *flag == 0)
4041     {
4042         ptr_d = block->xd;
4043         block->xd  = block->res;
4044     }
4045
4046     (*loc4)(block, *flag);
4047
4048     /* Implicit Solver & explicit block & flag==0 */
4049     /* adjust continuous state vector after call */
4050     if ((solver == IDA_BDF_Newton || solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes) && block->type < 10000 && *flag == 0)
4051     {
4052         block->xd  = ptr_d;
4053         if (flagi != 7)
4054         {
4055             for (k = 0; k < block->nx; k++)
4056             {
4057                 block->res[k] = block->res[k] - block->xd[k];
4058             }
4059         }
4060         else
4061         {
4062             for (k = 0; k < block->nx; k++)
4063             {
4064                 block->xd[k] = block->res[k];
4065             }
4066         }
4067     }
4068
4069     if (*flag < 0)
4070     {
4071         sciprint(_("Error in the Debug block \n"));
4072     }
4073 } /* call_debug_scicos */
4074 /*--------------------------------------------------------------------------*/
4075 /* simblk */
4076 static int simblk(realtype t, N_Vector yy, N_Vector yp, void *f_data)
4077 {
4078     double tx = 0., *x = NULL, *xd = NULL;
4079     int i = 0, nantest = 0;
4080
4081     tx = (double) t;
4082     x =  (double *) NV_DATA_S(yy);
4083     xd = (double *) NV_DATA_S(yp);
4084
4085     for (i = 0; i < *neq; i++)
4086     {
4087         xd[i] = 0;    /* à la place de "C2F(dset)(neq, &c_b14,xcdot , &c__1);"*/
4088     }
4089     C2F(ierode).iero = 0;
4090     *ierr = 0;
4091     odoit(&tx, x, xd, xd);
4092     C2F(ierode).iero = *ierr;
4093
4094     if (*ierr == 0)
4095     {
4096         nantest = 0;
4097         for (i = 0; i < *neq; i++) /* NaN checking */
4098         {
4099             if ((xd[i] - xd[i] != 0))
4100             {
4101                 sciprint(_("\nWarning: The computing function #%d returns a NaN/Inf"), i);
4102                 nantest = 1;
4103                 break;
4104             }
4105         }
4106         if (nantest == 1)
4107         {
4108             return 349;    /* recoverable error; */
4109         }
4110     }
4111
4112     return (abs(*ierr)); /* ierr>0 recoverable error; ierr>0 unrecoverable error; ierr=0: ok*/
4113
4114 } /* simblk */
4115 /*--------------------------------------------------------------------------*/
4116 /* grblk */
4117 static int grblk(realtype t, N_Vector yy, realtype *gout, void *g_data)
4118 {
4119     double tx = 0., *x = NULL;
4120     int jj = 0, nantest = 0;
4121
4122     tx = (double) t;
4123     x =  (double *) NV_DATA_S(yy);
4124
4125     C2F(ierode).iero = 0;
4126     *ierr = 0;
4127
4128     zdoit(&tx, x, x, (double*) gout);
4129
4130     if (*ierr == 0)
4131     {
4132         nantest = 0;
4133         for (jj = 0; jj < ng; jj++)
4134             if (gout[jj] - gout[jj] != 0)
4135             {
4136                 sciprint(_("\nWarning: The zero_crossing function #%d returns a NaN/Inf"), jj);
4137                 nantest = 1;
4138                 break;
4139             } /* NaN checking */
4140         if (nantest == 1)
4141         {
4142             return 350;    /* recoverable error; */
4143         }
4144     }
4145     C2F(ierode).iero = *ierr;
4146
4147     return 0;
4148 } /* grblk */
4149 /*--------------------------------------------------------------------------*/
4150 /* simblklsodar */
4151 static void simblklsodar(int * nequations, realtype * tOld, realtype * actual, realtype * res)
4152 {
4153     double tx = 0.;
4154     int i = 0;
4155
4156     tx = (double) * tOld;
4157
4158     for (i = 0; i < *nequations; ++i)
4159     {
4160         res[i] = 0;    /* à la place de "C2F(dset)(neq, &c_b14,xcdot , &c__1);"*/
4161     }
4162     C2F(ierode).iero = 0;
4163     *ierr = 0;
4164     odoit(&tx, actual, res, res);
4165     C2F(ierode).iero = *ierr;
4166
4167     if (*ierr == 0)
4168     {
4169         for (i = 0; i < *nequations; i++) /* NaN checking */
4170         {
4171             if ((res[i] - res[i] != 0))
4172             {
4173                 sciprint(_("\nWarning: The computing function #%d returns a NaN/Inf"), i);
4174             }
4175         }
4176     }
4177 } /* simblklsodar */
4178 /*--------------------------------------------------------------------------*/
4179 /* grblklsodar */
4180 static void grblklsodar(int * nequations, realtype * tOld, realtype * actual, int * ngc, realtype * res)
4181 {
4182     double tx = 0.;
4183     int jj = 0;
4184
4185     tx = (double) * tOld;
4186
4187     C2F(ierode).iero = 0;
4188     *ierr = 0;
4189
4190     zdoit(&tx, actual, actual, res);
4191
4192     if (*ierr == 0)
4193     {
4194         for (jj = 0; jj < *ngc; jj++)
4195         {
4196             if (res[jj] - res[jj] != 0)
4197             {
4198                 sciprint(_("\nWarning: The zero_crossing function #%d returns a NaN/Inf"), jj);
4199             } /* NaN checking */
4200         }
4201     }
4202 } /* grblklsodar */
4203 /*--------------------------------------------------------------------------*/
4204 /* simblkdaskr */
4205 static int simblkdaskr(realtype tres, N_Vector yy, N_Vector yp, N_Vector resval, void *rdata)
4206 {
4207     double tx = 0.;
4208     double *xc = NULL, *xcdot = NULL, *residual = NULL;
4209     realtype alpha = 0.;
4210
4211     UserData data;
4212
4213     realtype hh = 0.;
4214     int qlast = 0;
4215     int jj = 0, flag = 0, nantest = 0;
4216
4217     data = (UserData) rdata;
4218
4219     if (get_phase_simulation() == 1)
4220     {
4221         /* Just to update mode in a very special case, i.e., when initialization using modes fails.
4222         in this case, we relax all modes and try again one more time.
4223         */
4224         zdoit(&tx, NV_DATA_S(yy), NV_DATA_S(yp), NULL);
4225     }
4226
4227     hh = ZERO;
4228     flag = IDAGetCurrentStep(data->dae_mem, &hh);
4229     if (flag < 0)
4230     {
4231         *ierr = 200 + (-flag);
4232         return (*ierr);
4233     };
4234
4235     qlast = 0;
4236     flag = IDAGetCurrentOrder(data->dae_mem, &qlast);
4237     if (flag < 0)
4238     {
4239         *ierr = 200 + (-flag);
4240         return (*ierr);
4241     };
4242
4243     alpha = ZERO;
4244     for (jj = 0; jj < qlast; jj++)
4245     {
4246         alpha = alpha - ONE / (jj + 1);
4247     }
4248     if (hh != 0)
4249         // CJ=-alpha/hh;  // For function Get_Jacobian_cj
4250     {
4251         CJJ = -alpha / hh;
4252     }
4253     else
4254     {
4255         *ierr = 217;
4256         return (*ierr);
4257     }
4258     xc = (double *)  NV_DATA_S(yy);
4259     xcdot = (double *) NV_DATA_S(yp);
4260     residual = (double *) NV_DATA_S(resval);
4261     tx = (double) tres;
4262
4263     C2F(dcopy)(neq, xcdot, &c__1, residual, &c__1);
4264     *ierr = 0;
4265     C2F(ierode).iero = 0;
4266     odoit(&tx, xc, xcdot, residual);
4267
4268     C2F(ierode).iero = *ierr;
4269
4270     if (*ierr == 0)
4271     {
4272         nantest = 0;
4273         for (jj = 0; jj < *neq; jj++)
4274             if (residual[jj] - residual[jj] != 0) /* NaN checking */
4275             {
4276                 //sciprint(_("\nWarning: The residual function #%d returns a NaN"),jj);
4277                 nantest = 1;
4278                 break;
4279             }
4280         if (nantest == 1)
4281         {
4282             return 257;    /* recoverable error; */
4283         }
4284     }
4285
4286     return (abs(*ierr)); /* ierr>0 recoverable error; ierr>0 unrecoverable error; ierr=0: ok*/
4287 }/* simblkdaskr */
4288 /*--------------------------------------------------------------------------*/
4289 /* grblkdaskr */
4290 static int grblkdaskr(realtype t, N_Vector yy, N_Vector yp, realtype *gout, void *g_data)
4291 {
4292     double tx = 0.;
4293     int jj = 0, nantest = 0;
4294
4295     tx = (double) t;
4296
4297     *ierr = 0;
4298     C2F(ierode).iero = 0;
4299     zdoit(&tx, NV_DATA_S(yy), NV_DATA_S(yp), (double *) gout);
4300     if (*ierr == 0)
4301     {
4302         nantest = 0; /* NaN checking */
4303         for (jj = 0; jj < ng; jj++)
4304         {
4305             if (gout[jj] - gout[jj] != 0)
4306             {
4307                 sciprint(_("\nWarning: The zero-crossing function #%d returns a NaN"), jj);
4308                 nantest = 1;
4309                 break;
4310             }
4311         }
4312         if (nantest == 1)
4313         {
4314             return 258; /* recoverable error; */
4315         }
4316     }
4317     C2F(ierode).iero = *ierr;
4318     return (*ierr);
4319 }/* grblkdaskr */
4320 /*--------------------------------------------------------------------------*/
4321 /* simblkddaskr */
4322 static void simblkddaskr(realtype *tOld, realtype *actual, realtype *actualP, realtype *res, int *flag, double *dummy1, int *dummy2)
4323 {
4324     double tx = 0.;
4325
4326     int jj = 0;
4327
4328     if (get_phase_simulation() == 1)
4329     {
4330         /* Just to update mode in a very special case, i.e., when initialization using modes fails.
4331         in this case, we relax all modes and try again one more time.
4332         */
4333         zdoit(&tx, actual, actualP, NULL);
4334     }
4335
4336     CJJ = 6;
4337
4338     tx = (double) * tOld;
4339     *flag = 0;
4340
4341     C2F(dcopy)(neq, actualP, &c__1, res, &c__1);
4342     *ierr = 0;
4343     C2F(ierode).iero = 0;
4344     odoit(&tx, actual, actualP, res);
4345     C2F(ierode).iero = *ierr;
4346
4347     if (*ierr == 0)
4348     {
4349         for (jj = 0; jj < *neq; jj++)
4350             if (res[jj] - res[jj] != 0) /* NaN checking */
4351             {
4352                 sciprint(_("\nWarning: The residual function #%d returns a NaN"), jj);
4353                 *flag = -1; /* recoverable error; */
4354                 return;
4355             }
4356     }
4357     else
4358     {
4359         *flag = -2;
4360         return;
4361     }
4362     /* *flag=-1 recoverable error; *flag=-2 unrecoverable error; *flag=0: ok*/
4363
4364 }/* simblkddaskr */
4365 /*--------------------------------------------------------------------------*/
4366 /* grblkddaskr */
4367 static void grblkddaskr(int *nequations, realtype *tOld, realtype *actual, int *ngc, realtype *res, double *dummy1, int *dummy2)
4368 {
4369     double tx = 0.;
4370     int jj = 0;
4371
4372     tx = (double) * tOld;
4373
4374     *ierr = 0;
4375     C2F(ierode).iero = 0;
4376     zdoit(&tx, actual, actual, res);
4377     C2F(ierode).iero = *ierr;
4378
4379     if (*ierr == 0)
4380     {
4381         /* NaN checking */
4382         for (jj = 0; jj < *ngc; jj++)
4383         {
4384             if (res[jj] - res[jj] != 0)
4385             {
4386                 sciprint(_("\nWarning: The zero-crossing function #%d returns a NaN"), jj);
4387                 return;
4388             }
4389         }
4390     }
4391     else
4392     {
4393         sciprint(_("\nError: Problem in the evaluation of a root function"));
4394         return;
4395     }
4396
4397 }/* grblkddaskr */
4398 /*--------------------------------------------------------------------------*/
4399 /* jacpsol */
4400 static void jacpsol(realtype *res, int *ires, int *neq, realtype *tOld, realtype *actual, realtype *actualP,
4401                     realtype *rewt, realtype *savr, realtype *wk, realtype *h, realtype *cj, realtype *wp, int *iwp,
4402                     int *ier, double *dummy1, int *dummy2)
4403 {
4404     /* Here, we compute the system preconditioner matrix P, which is actually the jacobian matrix,
4405        so P(i,j) = dres(i)/dactual(j) + cj*dres(i)/dactualP(j), and we LU-decompose it. */
4406     int i = 0, j = 0, nrow = 0, info = 0;
4407     realtype tx = 0, del = 0, delinv = 0, ysave = 0, ypsave = 0;
4408     realtype * e = NULL;
4409
4410     tx = *tOld;
4411
4412     /* Work array used to evaluate res(*tOld, actual + small_increment, actualP + small_increment).
4413        savr already contains res(*tOld, actual, actualP). */
4414     e = (realtype *) calloc(*neq, sizeof(realtype));
4415
4416     for (i = 0; i < *neq; ++i)
4417     {
4418         del =  max (SQuround * max(fabs(actual[i]), fabs(*h * actualP[i])), 1. / rewt[i]);
4419         del *= (*h * actualP[i] >= 0) ? 1 : -1;
4420         del =  (actual[i] + del) - actual[i];
4421         ysave  = actual[i];
4422         ypsave = actualP[i];
4423         actual[i]  += del;
4424         actualP[i] += *cj * del;
4425         *ierr = 0;
4426         C2F(ierode).iero = 0;
4427         simblkddaskr(tOld, actual, actualP, e, ires, dummy1, dummy2);
4428         C2F(ierode).iero = *ierr;
4429         if (*ires == 0)
4430         {
4431             delinv = 1. / del;
4432             for (j = 0; j < *neq; ++j)
4433             {
4434                 wp[nrow + j] = (e[j] - savr[j]) * delinv;
4435                 /* NaN test */
4436                 if (wp[nrow + j] - wp[nrow + j] != 0)
4437                 {
4438                     sciprint(_("\nWarning: The preconditioner evaluation function returns a NaN at index #%d."), nrow + j);
4439                     *ier = -1;
4440                 }
4441             }
4442             nrow       += *neq;
4443             actual[i]  =  ysave;
4444             actualP[i] =  ypsave;
4445         }
4446         else
4447         {
4448             sciprint(_("\nError: The preconditioner evaluation function failed."));
4449             *ier = -1;
4450             free(e);
4451             return;
4452         }
4453     }
4454
4455     /* Proceed to LU factorization of P. */
4456     C2F(dgefa) (wp, neq, neq, iwp, &info);
4457     if (info != 0)
4458     {
4459         sciprint(_("\nError: Failed to factor the preconditioner."));
4460         *ier = -1;
4461     }
4462
4463     free(e);
4464 }/* jacpsol */
4465 /*--------------------------------------------------------------------------*/
4466 /* psol */
4467 static void psol(int *neq, realtype *tOld, realtype *actual, realtype *actualP,
4468                  realtype *savr, realtype *wk, realtype *cj, realtype *wght, realtype *wp,
4469                  int *iwp, realtype *b, realtype *eplin, int *ier, double *dummy1, int *dummy2)
4470 {
4471     /* This function "applies" the inverse of the preconditioner to 'b' (computes P^-1*b).
4472        It is done by solving P*x = b using the linpack routine 'dgesl'. */
4473     int i = 0, job = 0;
4474
4475     C2F(dgesl) (wp, neq, neq, iwp, b, &job);
4476
4477     /* NaN test */
4478     for (i = 0; i < *neq; ++i)
4479     {
4480         if (b[i] - b[i] != 0)
4481         {
4482             sciprint(_("\nWarning: The preconditioner application function returns a NaN at index #%d."), i);
4483             /* Indicate a recoverable error, meaning that the step will be retried with the same step size
4484                but with a call to 'jacpsol' to update necessary data, unless the Jacobian data is current,
4485                in which case the step will be retried with a smaller step size. */
4486             *ier = 1;
4487         }
4488     }
4489 }/* psol */
4490 /*--------------------------------------------------------------------------*/
4491 /* Subroutine addevs */
4492 static void addevs(double t, int *evtnb, int *ierr1)
4493 {
4494     static int i = 0, j = 0;
4495
4496     /* Function Body */
4497     *ierr1 = 0;
4498     if (evtspt[*evtnb] != -1)
4499     {
4500         if ((evtspt[*evtnb] == 0) && (*pointi == *evtnb))
4501         {
4502             tevts[*evtnb] = t;
4503             return;
4504         }
4505         else
4506         {
4507             if (*pointi == *evtnb)
4508             {
4509                 *pointi = evtspt[*evtnb]; /* remove from chain */
4510             }
4511             else
4512             {
4513                 i = *pointi;
4514                 while (*evtnb != evtspt[i])
4515                 {
4516                     i = evtspt[i];
4517                 }
4518                 evtspt[i] = evtspt[*evtnb]; /* remove old evtnb from chain */
4519                 if (TCritWarning == 0)
4520                 {
4521                     sciprint(_("\n Warning: an event is reprogrammed at t=%g by removing another"), t );
4522                     sciprint(_("\n         (already programmed) event. There may be an error in"));
4523                     sciprint(_("\n         your model. Please check your model\n"));
4524                     TCritWarning = 1;
4525                 }
4526                 do_cold_restart(); /* the erased event could be a critical
4527                                                                    event, so do_cold_restart is added to
4528                                                                    refresh the critical event table */
4529             }
4530             evtspt[*evtnb] = 0;
4531             tevts[*evtnb] = t;
4532         }
4533     }
4534     else
4535     {
4536         evtspt[*evtnb] = 0;
4537         tevts[*evtnb] = t;
4538     }
4539     if (*pointi == 0)
4540     {
4541         *pointi = *evtnb;
4542         return;
4543     }
4544     if (t < tevts[*pointi])
4545     {
4546         evtspt[*evtnb] = *pointi;
4547         *pointi = *evtnb;
4548         return;
4549     }
4550     i = *pointi;
4551
4552 L100:
4553     if (evtspt[i] == 0)
4554     {
4555         evtspt[i] = *evtnb;
4556         return;
4557     }
4558     if (t >= tevts[evtspt[i]])
4559     {
4560         j = evtspt[i];
4561         if (evtspt[j] == 0)
4562         {
4563             evtspt[j] = *evtnb;
4564             return;
4565         }
4566         i = j;
4567         goto L100;
4568     }
4569     else
4570     {
4571         evtspt[*evtnb] = evtspt[i];
4572         evtspt[i] = *evtnb;
4573     }
4574 } /* addevs */
4575 /*--------------------------------------------------------------------------*/
4576 /* Subroutine putevs */
4577 void putevs(double *t, int *evtnb, int *ierr1)
4578 {
4579     /* Function Body */
4580     *ierr1 = 0;
4581     if (evtspt[*evtnb] != -1)
4582     {
4583         *ierr1 = 1;
4584         return;
4585     }
4586     else
4587     {
4588         evtspt[*evtnb] = 0;
4589         tevts[*evtnb] = *t;
4590     }
4591     if (*pointi == 0)
4592     {
4593         *pointi = *evtnb;
4594         return;
4595     }
4596     evtspt[*evtnb] = *pointi;
4597     *pointi = *evtnb;
4598 } /* putevs */
4599 /*--------------------------------------------------------------------------*/
4600 /* Subroutine idoit */
4601 static void idoit(double *told)
4602 {
4603     /* initialisation (propagation of constant blocks outputs) */
4604     /*     Copyright INRIA */
4605
4606     int i2 = 0;
4607     scicos_flag flag = 0;
4608     int i = 0, j = 0;
4609     int ierr1 = 0;
4610
4611     ScicosImport *scs_imp = NULL;
4612     int *kf = NULL;
4613
4614     scs_imp = getscicosimportptr();
4615
4616     flag = 1;
4617     for (j = 0; j < * (scs_imp->niord); j++)
4618     {
4619         kf = &scs_imp->iord[j];
4620         C2F(curblk).kfun = *kf; /* */
4621         if (scs_imp->funtyp[*kf - 1] > -1)
4622         {
4623             /* continuous state */
4624             if (scs_imp->blocks[*kf - 1].nx > 0)
4625             {
4626                 scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4627                 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4628             }
4629             scs_imp->blocks[*kf - 1].nevprt = scs_imp->iord[j + * (scs_imp->niord)];
4630             /*callf(told, xd, x, x,x,&flag);*/
4631             callf(told, &scs_imp->blocks[*kf - 1], &flag);
4632             if (flag < 0)
4633             {
4634                 *ierr = 5 - flag;
4635                 return;
4636             }
4637         }
4638         if (scs_imp->blocks[*kf - 1].nevout > 0)
4639         {
4640             if (scs_imp->funtyp[*kf - 1] < 0)
4641             {
4642                 i = synchro_nev(scs_imp, *kf, ierr);
4643                 if (*ierr != 0)
4644                 {
4645                     return;
4646                 }
4647                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4648                 putevs(told, &i2, &ierr1);
4649                 if (ierr1 != 0)
4650                 {
4651                     /* event conflict */
4652                     *ierr = 3;
4653                     return;
4654                 }
4655                 doit(told);
4656                 if (*ierr != 0)
4657                 {
4658                     return;
4659                 }
4660             }
4661         }
4662     }
4663 } /* idoit_ */
4664 /*--------------------------------------------------------------------------*/
4665 /* Subroutine doit */
4666 static void doit(double *told)
4667 {
4668     /* propagation of blocks outputs on discrete activations */
4669     /*     Copyright INRIA */
4670
4671     int i = 0, i2 = 0;
4672     scicos_flag flag = 0;
4673     int nord = 0;
4674     int ierr1 = 0;
4675     int ii = 0, kever = 0;
4676
4677     ScicosImport *scs_imp = NULL;
4678     int *kf = NULL;
4679
4680     scs_imp = getscicosimportptr();
4681
4682     kever = *pointi;
4683     *pointi = evtspt[kever];
4684     evtspt[kever] = -1;
4685
4686     nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever - 1];
4687     if (nord == 0)
4688     {
4689         return;
4690     }
4691
4692     for (ii = scs_imp->ordptr[kever - 1]; ii <= scs_imp->ordptr[kever] - 1 ; ii++)
4693     {
4694         kf = &scs_imp->ordclk[ii - 1];
4695         C2F(curblk).kfun = *kf;
4696         if (scs_imp->funtyp[*kf - 1] > -1)
4697         {
4698             /* continuous state */
4699             if (scs_imp->blocks[*kf - 1].nx > 0)
4700             {
4701                 scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4702                 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4703             }
4704             scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
4705             flag = 1;
4706             /*callf(told, xd, x, x,x,&flag);*/
4707             callf(told, &scs_imp->blocks[*kf - 1], &flag);
4708             if (flag < 0)
4709             {
4710                 *ierr = 5 - flag;
4711                 return;
4712             }
4713         }
4714
4715         /* Initialize tvec */
4716         if (scs_imp->blocks[*kf - 1].nevout > 0)
4717         {
4718             if (scs_imp->funtyp[*kf - 1] < 0)
4719             {
4720                 i = synchro_nev(scs_imp, *kf, ierr);
4721                 if (*ierr != 0)
4722                 {
4723                     return;
4724                 }
4725                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4726                 putevs(told, &i2, &ierr1);
4727                 if (ierr1 != 0)
4728                 {
4729                     /* event conflict */
4730                     *ierr = 3;
4731                     return;
4732                 }
4733                 doit(told);
4734                 if (*ierr != 0)
4735                 {
4736                     return;
4737                 }
4738             }
4739         }
4740     }
4741 } /* doit_ */
4742 /*--------------------------------------------------------------------------*/
4743 /* Subroutine cdoit */
4744 static void cdoit(double *told)
4745 {
4746     /* propagation of continuous blocks outputs */
4747     /*     Copyright INRIA */
4748
4749     int i2 = 0;
4750     scicos_flag flag = 0;
4751     int ierr1 = 0;
4752     int i = 0, j = 0;
4753
4754     ScicosImport *scs_imp = NULL;
4755     int *kf = NULL;
4756
4757     scs_imp = getscicosimportptr();
4758
4759     /* Function Body */
4760     for (j = 0; j < * (scs_imp->ncord); j++)
4761     {
4762         kf = &scs_imp->cord[j];
4763         C2F(curblk).kfun = *kf;
4764         /* continuous state */
4765         if (scs_imp->blocks[*kf - 1].nx > 0)
4766         {
4767             scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4768             scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4769         }
4770         scs_imp->blocks[*kf - 1].nevprt = scs_imp->cord[j + * (scs_imp->ncord)];
4771         if (scs_imp->funtyp[*kf - 1] > -1)
4772         {
4773             flag = 1;
4774             /*callf(told, xd, x, x,x,&flag);*/
4775             callf(told, &scs_imp->blocks[*kf - 1], &flag);
4776             if (flag < 0)
4777             {
4778                 *ierr = 5 - flag;
4779                 return;
4780             }
4781         }
4782
4783         /* Initialize tvec */
4784         if (scs_imp->blocks[*kf - 1].nevout > 0)
4785         {
4786             if (scs_imp->funtyp[*kf - 1] < 0)
4787             {
4788                 i = synchro_nev(scs_imp, *kf, ierr);
4789                 if (*ierr != 0)
4790                 {
4791                     return;
4792                 }
4793                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4794                 putevs(told, &i2, &ierr1);
4795                 if (ierr1 != 0)
4796                 {
4797                     /* event conflict */
4798                     *ierr = 3;
4799                     return;
4800                 }
4801                 doit(told);
4802                 if (*ierr != 0)
4803                 {
4804                     return;
4805                 }
4806             }
4807         }
4808     }
4809 } /* cdoit_ */
4810 /*--------------------------------------------------------------------------*/
4811 /* Subroutine ddoit */
4812 static void ddoit(double *told)
4813 {
4814     /* update states & event out on discrete activations */
4815     /*     Copyright INRIA */
4816
4817     int i2 = 0, j = 0;
4818     scicos_flag flag = 0;
4819     int kiwa = 0;
4820     int i = 0, i3 = 0, ierr1 = 0;
4821     int ii = 0, keve = 0;
4822
4823     ScicosImport *scs_imp = NULL;
4824     int *kf = NULL;
4825
4826     scs_imp = getscicosimportptr();
4827
4828     /* Function Body */
4829     kiwa = 0;
4830     edoit(told, &kiwa);
4831     if (*ierr != 0)
4832     {
4833         return;
4834     }
4835
4836     /* update continuous and discrete states on event */
4837     if (kiwa == 0)
4838     {
4839         return;
4840     }
4841     for (i = 0; i < kiwa; i++)
4842     {
4843         keve = iwa[i];
4844         if (critev[keve] != 0)
4845         {
4846             hot = 0;
4847         }
4848         i2 = scs_imp->ordptr[keve] - 1;
4849         for (ii = scs_imp->ordptr[keve - 1]; ii <= i2; ii++)
4850         {
4851             kf = &scs_imp->ordclk[ii - 1];
4852             C2F(curblk).kfun = *kf;
4853             /* continuous state */
4854             if (scs_imp->blocks[*kf - 1].nx > 0)
4855             {
4856                 scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4857                 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4858             }
4859
4860             scs_imp->blocks[*kf - 1].nevprt = scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1];
4861
4862             if (scs_imp->blocks[*kf - 1].nevout > 0)
4863             {
4864                 if (scs_imp->funtyp[*kf - 1] >= 0)
4865                 {
4866                     /* initialize evout */
4867                     for (j = 0; j < scs_imp->blocks[*kf - 1].nevout; j++)
4868                     {
4869                         scs_imp->blocks[*kf - 1].evout[j] = -1;
4870                     }
4871                     flag = 3;
4872
4873                     if (scs_imp->blocks[*kf - 1].nevprt > 0) /* if event has continuous origin don't call*/
4874                     {
4875                         /*callf(told, xd, x, x ,x,&flag);*/
4876                         callf(told, &scs_imp->blocks[*kf - 1], &flag);
4877                         if (flag < 0)
4878                         {
4879                             *ierr = 5 - flag;
4880                             return;
4881                         }
4882                     }
4883
4884                     for (j = 0; j < scs_imp->blocks[*kf - 1].nevout; j++)
4885                     {
4886                         if (scs_imp->blocks[*kf - 1].evout[j] >= 0.)
4887                         {
4888                             i3 = j + scs_imp->clkptr[*kf - 1] ;
4889                             addevs(scs_imp->blocks[*kf - 1].evout[j] + (*told), &i3, &ierr1);
4890                             if (ierr1 != 0)
4891                             {
4892                                 /* event conflict */
4893                                 *ierr = 3;
4894                                 return;
4895                             }
4896                         }
4897                     }
4898                 }
4899             }
4900
4901             if (scs_imp->blocks[*kf - 1].nevprt > 0)
4902             {
4903                 if (scs_imp->blocks[*kf - 1].nx + scs_imp->blocks[*kf - 1].nz + scs_imp->blocks[*kf - 1].noz > 0 || \
4904                         *scs_imp->blocks[*kf - 1].work != NULL)
4905                 {
4906                     /*  if a hidden state exists, must also call (for new scope eg)  */
4907                     /*  to avoid calling non-real activations */
4908                     flag = 2;
4909                     /*callf(told, xd, x, x,x,&flag);*/
4910                     callf(told, &scs_imp->blocks[*kf - 1], &flag);
4911                     if (flag < 0)
4912                     {
4913                         *ierr = 5 - flag;
4914                         return;
4915                     }
4916                 }
4917             }
4918             else
4919             {
4920                 if (*scs_imp->blocks[*kf - 1].work != NULL)
4921                 {
4922                     flag = 2;
4923                     scs_imp->blocks[*kf - 1].nevprt = 0; /* in case some hidden continuous blocks need updating */
4924                     /*callf(told, xd, x, x,x,&flag);*/
4925                     callf(told, &scs_imp->blocks[*kf - 1], &flag);
4926                     if (flag < 0)
4927                     {
4928                         *ierr = 5 - flag;
4929                         return;
4930                     }
4931                 }
4932             }
4933         }
4934     }
4935 } /* ddoit_ */
4936 /*--------------------------------------------------------------------------*/
4937 /* Subroutine edoit */
4938 static void edoit(double *told, int *kiwa)
4939 {
4940     /* update blocks output on discrete activations */
4941     /*     Copyright INRIA */
4942
4943     int i2 = 0;
4944     scicos_flag flag = 0;
4945     int ierr1 = 0, i = 0;
4946     int kever = 0, ii = 0;
4947
4948     ScicosImport *scs_imp = NULL;
4949     int *kf = NULL;
4950     int nord = 0;
4951
4952     scs_imp = getscicosimportptr();
4953
4954     /* Function Body */
4955     kever = *pointi;
4956
4957     *pointi = evtspt[kever];
4958     evtspt[kever] = -1;
4959
4960     nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever - 1];
4961     if (nord == 0)
4962     {
4963         return;
4964     }
4965     iwa[*kiwa] = kever;
4966     ++(*kiwa);
4967     for (ii = scs_imp->ordptr[kever - 1]; ii <= scs_imp->ordptr[kever] - 1; ii++)
4968     {
4969         kf = &scs_imp->ordclk[ii - 1];
4970         C2F(curblk).kfun = *kf;
4971
4972         if (scs_imp->funtyp[*kf - 1] > -1)
4973         {
4974             /* continuous state */
4975             if (scs_imp->blocks[*kf - 1].nx > 0)
4976             {
4977                 scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4978                 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4979             }
4980
4981             scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
4982
4983             flag = 1;
4984             /*callf(told, xd, x, x,x,&flag);*/
4985             callf(told, &scs_imp->blocks[*kf - 1], &flag);
4986             if (flag < 0)
4987             {
4988                 *ierr = 5 - flag;
4989                 return;
4990             }
4991         }
4992
4993         /* Initialize tvec */
4994         if (scs_imp->blocks[*kf - 1].nevout > 0)
4995         {
4996             if (scs_imp->funtyp[*kf - 1] < 0)
4997             {
4998                 i = synchro_nev(scs_imp, *kf, ierr);
4999                 if (*ierr != 0)
5000                 {
5001                     return;
5002                 }
5003                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5004                 putevs(told, &i2, &ierr1);
5005                 if (ierr1 != 0)
5006                 {
5007                     /* event conflict */
5008                     *ierr = 3;
5009                     return;
5010                 }
5011                 edoit(told, kiwa);
5012                 if (*ierr != 0)
5013                 {
5014                     return;
5015                 }
5016             }
5017         }
5018     }
5019 } /* edoit_ */
5020 /*--------------------------------------------------------------------------*/
5021 /* Subroutine odoit */
5022 static void odoit(double *told, double *xt, double *xtd, double *residual)
5023 {
5024     /* update blocks derivative of continuous time block */
5025     /*     Copyright INRIA */
5026
5027     int i2 = 0;
5028     scicos_flag flag = 0;
5029     int keve = 0, kiwa = 0;
5030     int ierr1 = 0, i = 0;
5031     int ii = 0, jj = 0;
5032
5033     ScicosImport *scs_imp = NULL;
5034     int *kf = NULL;
5035
5036     scs_imp = getscicosimportptr();
5037
5038     /* Function Body */
5039     kiwa = 0;
5040     for (jj = 0; jj < * (scs_imp->noord); jj++)
5041     {
5042         kf = &scs_imp->oord[jj];
5043         C2F(curblk).kfun = *kf;
5044         /* continuous state */
5045         if (scs_imp->blocks[*kf - 1].nx > 0)
5046         {
5047             scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5048             scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5049             scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5050         }
5051
5052         scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[jj + * (scs_imp->noord)];
5053         if (scs_imp->funtyp[*kf - 1] > -1)
5054         {
5055             flag = 1;
5056             /*callf(told, xtd, xt, residual,x,&flag);*/
5057             callf(told, &scs_imp->blocks[*kf - 1], &flag);
5058             if (flag < 0)
5059             {
5060                 *ierr = 5 - flag;
5061                 return;
5062             }
5063         }
5064
5065         if (scs_imp->blocks[*kf - 1].nevout > 0)
5066         {
5067             if (scs_imp->funtyp[*kf - 1] < 0)
5068             {
5069                 if (scs_imp->blocks[*kf - 1].nmode > 0)
5070                 {
5071                     i2 = scs_imp->blocks[*kf - 1].mode[0] + scs_imp->clkptr[*kf - 1] - 1;
5072                 }
5073                 else
5074                 {
5075                     i = synchro_nev(scs_imp, *kf, ierr);
5076                     if (*ierr != 0)
5077                     {
5078                         return;
5079                     }
5080                     i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5081                 }
5082                 putevs(told, &i2, &ierr1);
5083                 if (ierr1 != 0)
5084                 {
5085                     /* event conflict */
5086                     *ierr = 3;
5087                     return;
5088                 }
5089                 ozdoit(told, xt, xtd, &kiwa);
5090                 if (*ierr != 0)
5091                 {
5092                     return;
5093                 }
5094             }
5095         }
5096     }
5097
5098     /*  update states derivatives */
5099     for (ii = 0; ii < * (scs_imp->noord); ii++)
5100     {
5101         kf = &scs_imp->oord[ii];
5102         C2F(curblk).kfun = *kf;
5103         if (scs_imp->blocks[*kf - 1].nx > 0 || \
5104                 *scs_imp->blocks[*kf - 1].work != NULL)
5105         {
5106             /* work tests if a hidden state exists, used for delay block */
5107             flag = 0;
5108             /* continuous state */
5109             if (scs_imp->blocks[*kf - 1].nx > 0)
5110             {
5111                 scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5112                 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5113                 scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5114             }
5115             scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[ii + * (scs_imp->noord)];
5116             /*callf(told, xtd, xt, residual,xt,&flag);*/
5117             callf(told, &scs_imp->blocks[*kf - 1], &flag);
5118
5119             if (flag < 0)
5120             {
5121                 *ierr = 5 - flag;
5122                 return;
5123             }
5124         }
5125     }
5126
5127     for (i = 0; i < kiwa; i++)
5128     {
5129         keve = iwa[i];
5130         for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1; ii++)
5131         {
5132             kf = &scs_imp->ordclk[ii - 1];
5133             C2F(curblk).kfun = *kf;
5134             if (scs_imp->blocks[*kf - 1].nx > 0 || \
5135                     *scs_imp->blocks[*kf - 1].work != NULL)
5136             {
5137                 /* work tests if a hidden state exists */
5138                 flag = 0;
5139                 /* continuous state */
5140                 if (scs_imp->blocks[*kf - 1].nx > 0)
5141                 {
5142                     scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5143                     scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5144                     scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5145                 }
5146                 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5147                 /*callf(told, xtd, xt, residual,xt,&flag);*/
5148                 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5149
5150                 if (flag < 0)
5151                 {
5152                     *ierr = 5 - flag;
5153                     return;
5154                 }
5155             }
5156         }
5157     }
5158 } /* odoit_ */
5159 /*--------------------------------------------------------------------------*/
5160 /* Subroutine reinitdoit */
5161 static void reinitdoit(double *told)
5162 {
5163     /* update blocks xproperties of continuous time block */
5164     /*     Copyright INRIA */
5165
5166     int i2 = 0;
5167     scicos_flag flag = 0;
5168     int keve = 0, kiwa = 0;
5169     int ierr1 = 0, i = 0;
5170     int ii = 0, jj = 0;
5171
5172     ScicosImport *scs_imp = NULL;
5173     int *kf = NULL;
5174
5175     scs_imp = getscicosimportptr();
5176
5177     /* Function Body */
5178     kiwa = 0;
5179     for (jj = 0; jj < * (scs_imp->noord); jj++)
5180     {
5181         kf = &scs_imp->oord[jj];
5182         C2F(curblk).kfun = *kf;
5183         /* continuous state */
5184         if (scs_imp->blocks[*kf - 1].nx > 0)
5185         {
5186             scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
5187             scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5188         }
5189         scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[jj + * (scs_imp->noord)];
5190         if (scs_imp->funtyp[*kf - 1] > -1)
5191         {
5192             flag = 1;
5193             /*callf(told, xd, x, x,x,&flag);*/
5194             callf(told, &scs_imp->blocks[*kf - 1], &flag);
5195             if (flag < 0)
5196             {
5197                 *ierr = 5 - flag;
5198                 return;
5199             }
5200         }
5201
5202         if (scs_imp->blocks[*kf - 1].nevout > 0 && scs_imp->funtyp[*kf - 1] < 0)
5203         {
5204             i = synchro_nev(scs_imp, *kf, ierr);
5205             if (*ierr != 0)
5206             {
5207                 return;
5208             }
5209             if (scs_imp->blocks[*kf - 1].nmode > 0)
5210             {
5211                 scs_imp->blocks[*kf - 1].mode[0] = i;
5212             }
5213             i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5214             putevs(told, &i2, &ierr1);
5215             if (ierr1 != 0)
5216             {
5217                 /* event conflict */
5218                 *ierr = 3;
5219                 return;
5220             }
5221             doit(told);
5222             if (*ierr != 0)
5223             {
5224                 return;
5225             }
5226         }
5227     }
5228
5229     /* reinitialize */
5230     for (ii = 0; ii < * (scs_imp->noord); ii++)
5231     {
5232         kf = &scs_imp->oord[ii];
5233         C2F(curblk).kfun = *kf;
5234         if (scs_imp->blocks[*kf - 1].nx > 0)
5235         {
5236             flag = 7;
5237             scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
5238             scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5239             scs_imp->blocks[*kf - 1].res = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5240             scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[ii + * (scs_imp->noord)];
5241             scs_imp->blocks[*kf - 1].xprop = &scs_imp->xprop[-1 + scs_imp->xptr[*kf - 1]];
5242             /*callf(told, xd, x, xd,x,&flag);*/
5243             callf(told, &scs_imp->blocks[*kf - 1], &flag);
5244             if (flag < 0)
5245             {
5246                 *ierr = 5 - flag;
5247                 return;
5248             }
5249         }
5250     }
5251
5252     for (i = 0; i < kiwa; i++)
5253     {
5254         keve = iwa[i];
5255         for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1; ii++)
5256         {
5257             kf = &scs_imp->ordclk[ii - 1];
5258             C2F(curblk).kfun = *kf;
5259             if (scs_imp->blocks[*kf - 1].nx > 0)
5260             {
5261                 flag = 7;
5262                 scs_imp->blocks[*kf - 1].x  = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
5263                 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5264                 scs_imp->blocks[*kf - 1].res = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5265                 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5266                 scs_imp->blocks[*kf - 1].xprop = &scs_imp->xprop[-1 + scs_imp->xptr[*kf - 1]];
5267                 /*callf(told, xd, x, xd,x,&flag);*/
5268                 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5269                 if (flag < 0)
5270                 {
5271                     *ierr = 5 - flag;
5272                     return;
5273                 }
5274             }
5275         }
5276     }
5277 } /* reinitdoit_ */
5278 /*--------------------------------------------------------------------------*/
5279 /* Subroutine ozdoit */
5280 static void ozdoit(double *told, double *xt, double *xtd, int *kiwa)
5281 {
5282     /* update blocks output of continuous time block on discrete activations */
5283     /*     Copyright INRIA */
5284
5285     int i2 = 0;
5286     scicos_flag flag = 0;
5287     int nord = 0;
5288     int ierr1 = 0, i = 0;
5289     int ii = 0, kever = 0;
5290
5291     ScicosImport *scs_imp = NULL;
5292     int *kf = NULL;
5293
5294     scs_imp = getscicosimportptr();
5295
5296     /* Function Body */
5297     kever = *pointi;
5298     *pointi = evtspt[kever];
5299     evtspt[kever] = -1;
5300
5301     nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever - 1];
5302     if (nord == 0)
5303     {
5304         return;
5305     }
5306     iwa[*kiwa] = kever;
5307     ++(*kiwa);
5308
5309     for (ii = scs_imp->ordptr[kever - 1]; ii <= scs_imp->ordptr[kever] - 1; ii++)
5310     {
5311         kf = &scs_imp->ordclk[ii - 1];
5312         C2F(curblk).kfun = *kf;
5313         if (scs_imp->funtyp[*kf - 1] > -1)
5314         {
5315             /* continuous state */
5316             if (scs_imp->blocks[*kf - 1].nx > 0)
5317             {
5318                 scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5319                 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5320             }
5321             scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5322             flag = 1;
5323             /*callf(told, xtd, xt, xt,x,&flag);*/
5324             callf(told, &scs_imp->blocks[*kf - 1], &flag);
5325             if (flag < 0)
5326             {
5327                 *ierr = 5 - flag;
5328                 return;
5329             }
5330         }
5331
5332         /* Initialize tvec */
5333         if (scs_imp->blocks[*kf - 1].nevout > 0)
5334         {
5335             if (scs_imp->funtyp[*kf - 1] < 0)
5336             {
5337                 if (phase == 1 || scs_imp->blocks[*kf - 1].nmode == 0)
5338                 {
5339                     i = synchro_nev(scs_imp, *kf, ierr);
5340                     if (*ierr != 0)
5341                     {
5342                         return;
5343                     }
5344                 }
5345                 else
5346                 {
5347                     i = scs_imp->blocks[*kf - 1].mode[0];
5348                 }
5349                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5350                 putevs(told, &i2, &ierr1);
5351                 if (ierr1 != 0)
5352                 {
5353                     /* event conflict */
5354                     *ierr = 3;
5355                     return;
5356                 }
5357                 ozdoit(told, xt, xtd, kiwa);
5358             }
5359         }
5360     }
5361 } /* ozdoit_ */
5362 /*--------------------------------------------------------------------------*/
5363 /* Subroutine zdoit */
5364 static void zdoit(double *told, double *xt, double *xtd, double *g)
5365 {
5366     /* update blocks zcross of continuous time block  */
5367     /*     Copyright INRIA */
5368     int i2 = 0;
5369     scicos_flag flag = 0;
5370     int keve = 0, kiwa = 0;
5371     int ierr1 = 0, i = 0, j = 0;
5372     int ii = 0, jj = 0;
5373
5374     ScicosImport *scs_imp = NULL;
5375     int *kf = NULL;
5376
5377     scs_imp = getscicosimportptr();
5378
5379     /* Function Body */
5380     for (i = 0; i < * (scs_imp->ng); i++)
5381     {
5382         g[i] = 0.;
5383     }
5384
5385     kiwa = 0;
5386     for (jj = 0; jj < * (scs_imp->nzord); jj++)
5387     {
5388         kf = &scs_imp->zord[jj];
5389         C2F(curblk).kfun = *kf;
5390         /* continuous state */
5391         if (scs_imp->blocks[*kf - 1].nx > 0)
5392         {
5393             scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5394             scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5395         }
5396         scs_imp->blocks[*kf - 1].nevprt = scs_imp->zord[jj + * (scs_imp->nzord)];
5397
5398         if (scs_imp->funtyp[*kf - 1] > -1)
5399         {
5400             flag = 1;
5401             /*callf(told, xtd, xt, xt,xt,&flag);*/
5402             callf(told, &scs_imp->blocks[*kf - 1], &flag);
5403             if (flag < 0)
5404             {
5405                 *ierr = 5 - flag;
5406                 return;
5407             }
5408         }
5409
5410         /* Initialize tvec */
5411         if (scs_imp->blocks[*kf - 1].nevout > 0)
5412         {
5413             if (scs_imp->funtyp[*kf - 1] < 0)
5414             {
5415                 if (phase == 1 || scs_imp->blocks[*kf - 1].nmode == 0)
5416                 {
5417                     i = synchro_nev(scs_imp, *kf, ierr);
5418                     if (*ierr != 0)
5419                     {
5420                         return;
5421                     }
5422                 }
5423                 else
5424                 {
5425                     i = scs_imp->blocks[*kf - 1].mode[0];
5426                 }
5427                 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5428                 putevs(told, &i2, &ierr1);
5429                 if (ierr1 != 0)
5430                 {
5431                     /* event conflict */
5432                     *ierr = 3;
5433                     return;
5434                 }
5435                 ozdoit(told, xt, xtd, &kiwa);
5436                 if (*ierr != 0)
5437                 {
5438                     return;
5439                 }
5440             }
5441         }
5442     }
5443
5444     /* update zero crossing surfaces */
5445     for (ii = 0; ii < * (scs_imp->nzord); ii++)
5446     {
5447         kf = &scs_imp->zord[ii];
5448         C2F(curblk).kfun = *kf;
5449         if (scs_imp->blocks[*kf - 1].ng > 0)
5450         {
5451             /* update g array ptr */
5452             scs_imp->blocks[*kf - 1].g = &g[scs_imp->zcptr[*kf - 1] - 1];
5453             if (scs_imp->funtyp[*kf - 1] > 0)
5454             {
5455                 flag = 9;
5456                 /* continuous state */
5457                 if (scs_imp->blocks[*kf - 1].nx > 0)
5458                 {
5459                     scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5460                     scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5461                 }
5462                 scs_imp->blocks[*kf - 1].nevprt = scs_imp->zord[ii + * (scs_imp->nzord)];
5463                 /*callf(told, xtd, xt, xtd,g,&flag);*/
5464                 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5465                 if (flag < 0)
5466                 {
5467                     *ierr = 5 - flag;
5468                     return;
5469                 }
5470             }
5471             else
5472             {
5473                 j = synchro_g_nev(scs_imp, g, *kf, ierr);
5474                 if (*ierr != 0)
5475                 {
5476                     return;
5477                 }
5478                 if ( (phase == 1) && (scs_imp->blocks[*kf - 1].nmode > 0) )
5479                 {
5480                     scs_imp->blocks[*kf - 1].mode[0] = j;
5481                 }
5482             }
5483
5484             // scs_imp->blocks[*kf-1].g = &scs_imp->g[scs_imp->zcptr[*kf]-1];
5485
5486         }
5487     }
5488
5489     for (i = 0; i < kiwa; i++)
5490     {
5491         keve = iwa[i];
5492         for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1; ii++)
5493         {
5494             kf = &scs_imp->ordclk[ii - 1];
5495             C2F(curblk).kfun = *kf;
5496             if (scs_imp->blocks[*kf - 1].ng > 0)
5497             {
5498                 /* update g array ptr */
5499                 scs_imp->blocks[*kf - 1].g = &g[scs_imp->zcptr[*kf - 1] - 1];
5500                 if (scs_imp->funtyp[*kf - 1] > 0)
5501                 {
5502                     flag = 9;
5503                     /* continuous state */
5504                     if (scs_imp->blocks[*kf - 1].nx > 0)
5505                     {
5506                         scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5507                         scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5508                     }
5509                     scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5510                     /*callf(told, xtd, xt, xtd,g,&flag);*/
5511                     callf(told, &scs_imp->blocks[*kf - 1], &flag);
5512                     if (flag < 0)
5513                     {
5514                         *ierr = 5 - flag;
5515                         return;
5516                     }
5517                 }
5518                 else
5519                 {
5520                     j = synchro_g_nev(scs_imp, g, *kf, ierr);
5521                     if (*ierr != 0)
5522                     {
5523                         return;
5524                     }
5525                     if ((phase == 1) && (scs_imp->blocks[*kf - 1].nmode > 0))
5526                     {
5527                         scs_imp->blocks[*kf - 1].mode[0] = j;
5528                     }
5529                 }
5530
5531                 //scs_imp->blocks[*kf-1].g = &scs_imp->g[scs_imp->zcptr[*kf]-1];
5532             }
5533         }
5534     }
5535 } /* zdoit_ */
5536 /*--------------------------------------------------------------------------*/
5537 /* Subroutine Jdoit */
5538 void Jdoit(double *told, double *xt, double *xtd, double *residual, int *job)
5539 {
5540     /* update blocks jacobian of continuous time block  */
5541     /*     Copyright INRIA */
5542
5543     int i2 = 0;
5544     scicos_flag flag = 0;
5545     int keve = 0, kiwa = 0;
5546     int ierr1 = 0, i = 0;
5547     int ii = 0, jj = 0;
5548
5549     ScicosImport *scs_imp = NULL;
5550     int *kf = NULL;
5551
5552     scs_imp = getscicosimportptr();
5553
5554     /* Function Body */
5555     kiwa = 0;
5556     for (jj = 0; jj < * (scs_imp->noord); jj++)
5557     {
5558         kf = &scs_imp->oord[jj];
5559         C2F(curblk).kfun = *kf;
5560         scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[jj + * (scs_imp->noord)];
5561         if (scs_imp->funtyp[*kf - 1] > -1)
5562         {
5563             flag = 1;
5564             /* applying desired output */
5565             if ((*job == 2) && (scs_imp->oord[jj] == AJacobian_block))
5566             {
5567             }
5568             else
5569                 /* continuous state */
5570                 if (scs_imp->blocks[*kf - 1].nx > 0)
5571                 {
5572                     scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5573                     scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5574                     scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5575                 }
5576
5577             /*callf(told, xtd, xt, residual,x,&flag);*/
5578             callf(told, &scs_imp->blocks[*kf - 1], &flag);
5579             if (flag < 0)
5580             {
5581                 *ierr = 5 - flag;
5582                 return;
5583             }
5584         }
5585
5586         if (scs_imp->blocks[*kf - 1].nevout > 0)
5587         {
5588             if (scs_imp->funtyp[*kf - 1] < 0)
5589             {
5590                 if (scs_imp->blocks[*kf - 1].nmode > 0)
5591                 {
5592                     i2 = scs_imp->blocks[*kf - 1].mode[0] + scs_imp->clkptr[*kf - 1] - 1;
5593                 }
5594                 else
5595                 {
5596                     i = synchro_nev(scs_imp, *kf, ierr);
5597                     if (*ierr != 0)
5598                     {
5599                         return;
5600                     }
5601                     i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5602                 }
5603                 putevs(told, &i2, &ierr1);
5604                 if (ierr1 != 0)
5605                 {
5606                     /* event conflict */
5607                     *ierr = 3;
5608                     return;
5609                 }
5610                 ozdoit(told, xt, xtd, &kiwa);
5611             }
5612         }
5613     }
5614
5615     /* update states derivatives */
5616     for (ii = 0; ii < * (scs_imp->noord); ii++)
5617     {
5618         kf = &scs_imp->oord[ii];
5619         C2F(curblk).kfun = *kf;
5620         if (scs_imp->blocks[*kf - 1].nx > 0 || \
5621                 *scs_imp->blocks[*kf - 1].work != NULL)
5622         {
5623             /* work tests if a hidden state exists, used for delay block */
5624             flag = 0;
5625             if (((*job == 1) && (scs_imp->oord[ii] == AJacobian_block)) || (*job != 1))
5626             {
5627                 if (*job == 1)
5628                 {
5629                     flag = 10;
5630                 }
5631                 /* continuous state */
5632                 if (scs_imp->blocks[*kf - 1].nx > 0)
5633                 {
5634                     scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5635                     scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5636                     scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5637                 }
5638                 scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[ii + * (scs_imp->noord)];
5639                 /*callf(told, xtd, xt, residual,xt,&flag);*/
5640                 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5641             }
5642             if (flag < 0)
5643             {
5644                 *ierr = 5 - flag;
5645                 return;
5646             }
5647         }
5648     }
5649
5650     for (i = 0; i < kiwa; i++)
5651     {
5652         keve = iwa[i];
5653         for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1;  ii++)
5654         {
5655             kf = &scs_imp->ordclk[ii - 1];
5656             C2F(curblk).kfun = *kf;
5657             if (scs_imp->blocks[*kf - 1].nx > 0 || \
5658                     *scs_imp->blocks[*kf - 1].work != NULL)
5659             {
5660                 /* work tests if a hidden state exists */
5661                 flag = 0;
5662                 if (((*job == 1) && (scs_imp->oord[ii - 1] == AJacobian_block)) || (*job != 1))
5663                 {
5664                     if (*job == 1)
5665                     {
5666                         flag = 10;
5667                     }
5668                     /* continuous state */
5669                     if (scs_imp->blocks[*kf - 1].nx > 0)
5670                     {
5671                         scs_imp->blocks[*kf - 1].x  = &xt[scs_imp->xptr[*kf - 1] - 1];
5672                         scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5673                         scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5674                     }
5675                     scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5676                     /*callf(told, xtd, xt, residual,xt,&flag);*/
5677                     callf(told, &scs_imp->blocks[*kf - 1], &flag);
5678                 }
5679                 if (flag < 0)
5680                 {
5681                     *ierr = 5 - flag;
5682                     return;
5683                 }
5684             }
5685         }
5686     }
5687 } /* Jdoit_ */
5688 /*--------------------------------------------------------------------------*/
5689 /* Subroutine synchro_nev */
5690 static int synchro_nev(ScicosImport *scs_imp, int kf, int *ierr)
5691 {
5692     /* synchro blocks computation  */
5693     /*     Copyright INRIA */
5694     SCSREAL_COP *outtbdptr = NULL;     /*to store double of outtb*/
5695     SCSINT8_COP *outtbcptr = NULL;     /*to store int8 of outtb*/
5696     SCSINT16_COP *outtbsptr = NULL;    /*to store int16 of outtb*/
5697     SCSINT32_COP *outtblptr = NULL;    /*to store int32 of outtb*/
5698     SCSUINT8_COP *outtbucptr = NULL;   /*to store unsigned int8 of outtb */
5699     SCSUINT16_COP *outtbusptr = NULL;  /*to store unsigned int16 of outtb */
5700     SCSUINT32_COP *outtbulptr = NULL;  /*to store unsigned int32 of outtb */
5701
5702     int cond = 0;
5703     int i = 0; /* return 0 by default */
5704
5705     /* variable for param */
5706     int *outtbtyp = 0;
5707     void **outtbptr = NULL;
5708     int *funtyp = 0;
5709     int *inplnk = 0;
5710     int *inpptr = 0;
5711
5712     /* get param ptr */
5713     outtbtyp = scs_imp->outtbtyp;
5714     outtbptr = scs_imp->outtbptr;
5715     funtyp = scs_imp->funtyp;
5716     inplnk = scs_imp->inplnk;
5717     inpptr = scs_imp->inpptr;
5718
5719     /* if-then-else blk */
5720     if (funtyp[kf - 1] == -1)
5721     {
5722         switch (outtbtyp[-1 + inplnk[inpptr[kf - 1] - 1]])
5723         {
5724             case SCSREAL_N    :
5725                 outtbdptr = (SCSREAL_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5726                 cond = (*outtbdptr <= 0.);
5727                 break;
5728
5729             case SCSCOMPLEX_N :
5730                 outtbdptr = (SCSCOMPLEX_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5731                 cond = (*outtbdptr <= 0.);
5732                 break;
5733
5734             case SCSINT8_N    :
5735                 outtbcptr = (SCSINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5736                 cond = (*outtbcptr <= 0);
5737                 break;
5738
5739             case SCSINT16_N   :
5740                 outtbsptr = (SCSINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5741                 cond = (*outtbsptr <= 0);
5742                 break;
5743
5744             case SCSINT32_N   :
5745                 outtblptr = (SCSINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5746                 cond = (*outtblptr <= 0);
5747                 break;
5748
5749             case SCSUINT8_N   :
5750                 outtbucptr = (SCSUINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5751                 cond = (*outtbucptr <= 0);
5752                 break;
5753
5754             case SCSUINT16_N  :
5755                 outtbusptr = (SCSUINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5756                 cond = (*outtbusptr <= 0);
5757                 break;
5758
5759             case SCSUINT32_N  :
5760                 outtbulptr = (SCSUINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5761                 cond = (*outtbulptr <= 0);
5762                 break;
5763
5764             default  : /* Add a message here */
5765                 *ierr = 25;
5766                 return 0;
5767         }
5768         if (cond)
5769         {
5770             i = 2;
5771         }
5772         else
5773         {
5774             i = 1;
5775         }
5776     }
5777     /* eselect blk */
5778     else if (funtyp[kf - 1] == -2)
5779     {
5780         switch (outtbtyp[-1 + inplnk[inpptr[kf - 1] - 1]])
5781         {
5782             case SCSREAL_N    :
5783                 outtbdptr = (SCSREAL_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5784                 i = Max(Min((int) * outtbdptr, scs_imp->blocks[kf - 1].nevout), 1);
5785                 break;
5786
5787             case SCSCOMPLEX_N :
5788                 outtbdptr = (SCSCOMPLEX_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5789                 i = Max(Min((int) * outtbdptr, scs_imp->blocks[kf - 1].nevout), 1);
5790                 break;
5791
5792             case SCSINT8_N    :
5793                 outtbcptr = (SCSINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5794                 i = Max(Min((int) * outtbcptr,
5795                             scs_imp->blocks[kf - 1].nevout), 1);
5796                 break;
5797
5798             case SCSINT16_N   :
5799                 outtbsptr = (SCSINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5800                 i = Max(Min((int) * outtbsptr, scs_imp->blocks[kf - 1].nevout), 1);
5801                 break;
5802
5803             case SCSINT32_N   :
5804                 outtblptr = (SCSINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5805                 i = Max(Min((int) * outtblptr, scs_imp->blocks[kf - 1].nevout), 1);
5806                 break;
5807