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