3 * Copyright (C) INRIA - METALAU Project <scicos@inria.fr>
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.
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.
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.
19 * See the file ./license.txt
27 * input/output regular ports
30 * added object paramaters/states
33 * CVODE (Sundials) replaced LSODAR
34 * IDA (Sundials) replaced DASKR
36 /*--------------------------------------------------------------------------*/
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 */
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 */
60 #include "machine.h" /* C2F */
61 #include "scicos-def.h"
62 #include "stack-def.h"
66 #include "scicos_internal.h"
68 #include "core_math.h"
69 #include "storeCommand.h"
72 #include "sci_malloc.h"
78 #include "dynlib_scicos.h"
80 #include "configvariable_interface.h" /* getEntryPointPosition() and getEntryPointFromPosition() */
82 #include "lsodar.h" /* prototypes for lsodar fcts. and consts. */
83 #include "ddaskr.h" /* prototypes for ddaskr fcts. and consts. */
85 #if defined(linux) && defined(__i386__)
86 #include "setPrecisionFPU.h"
89 #include "localization.h"
90 #include "charEncoding.h"
91 /*--------------------------------------------------------------------------*/
98 double *gwork; /* just added for a very special use: a
99 space passing to grblkdakr for zero crossing surfaces
100 when updating mode variables during initialization */
103 SCICOS_IMPEXP SCSPTR_struct C2F(scsptr);
109 CVode_BDF_Functional,
111 CVode_Adams_Functional,
114 Implicit_Runge_Kutta,
115 IDA_BDF_Newton = 100,
116 DDaskr_BDF_Newton = 101,
117 DDaskr_BDF_GMRes = 102
119 /*--------------------------------------------------------------------------*/
122 if (*neq>0) ODEFree(&ode_mem); \
123 if (*neq>0) N_VDestroy_Serial(y); \
124 if ( ng>0 ) FREE(jroot); \
125 if ( ng>0 ) FREE(zcros);
128 /* TJacque allocated by sundials */
130 if (*neq>0) free(TJacque); \
131 if (*neq>0) FREE(data->rwork); \
132 if (( ng>0 )&& (*neq>0)) FREE(data->gwork); \
133 if (*neq>0) N_VDestroy_Serial(data->ewt); \
134 if (*neq>0) FREE(data); \
135 if (*neq>0) DAEFree(&dae_mem); \
136 if (*neq>0) N_VDestroy_Serial(IDx); \
137 if (*neq>0) N_VDestroy_Serial(yp); \
138 if (*neq>0) N_VDestroy_Serial(yy); \
139 if ( ng>0 ) FREE(jroot); \
140 if ( ng>0 ) FREE(zcros); \
141 if (nmod>0) FREE(Mode_save);
143 #define freeouttbptr \
154 N_VDestroy_Serial(y); \
155 N_VDestroy_Serial(fscale); \
156 N_VDestroy_Serial(yscale); \
160 #define ONE RCONST(1.0)
161 #define ZERO RCONST(0.0)
162 #define T0 RCONST(0.0)
163 /*--------------------------------------------------------------------------*/
164 /* Table of constant values */
165 static int c__90 = 90;
166 static int c__91 = 91;
168 static double c_b14 = 0.;
171 int TCritWarning = 0;
173 static double *t0 = NULL, *tf = NULL;
174 static double *x = NULL, *tevts = NULL, *xd = NULL, *g = NULL;
175 static double Atol = 0., rtol = 0., ttol = 0., deltat = 0., hmax = 0.;
176 static int *ierr = NULL;
177 static int *pointi = NULL;
178 static int *xptr = NULL, *modptr = NULL, *evtspt = NULL;
179 static int *funtyp = NULL, *inpptr = NULL, *outptr = NULL, *inplnk = NULL, *outlnk = NULL;
180 static int *clkptr = NULL, *ordptr = NULL, *ordclk = NULL, *cord = NULL, *iord = NULL, *oord = NULL, *zord = NULL, *critev = NULL, *zcptr = NULL;
181 static int nblk = 0, nordptr = 0, nlnk = 0, nx = 0, ng = 0, ncord = 0, noord = 0, nzord = 0;
182 static int nordclk = 0, niord = 0, nmod = 0;
183 static int nelem = 0;
184 static int *mod = NULL;
185 static int *neq = NULL;
186 static int *xprop = NULL; /* xproperty */
187 static int debug_block = 0;
188 static int *iwa = NULL;
191 /* declaration of ptr for typed port */
192 static void **outtbptr = NULL; /*pointer array of object of outtb*/
193 static int *outtbsz = NULL; /*size of object of outtb*/
194 static int *outtbtyp = NULL; /*type of object of outtb*/
196 SCSREAL_COP *outtbdptr = NULL; /*to store double of outtb*/
197 SCSINT8_COP *outtbcptr = NULL; /*to store int8 of outtb*/
198 SCSINT16_COP *outtbsptr = NULL; /*to store int16 of outtb*/
199 SCSINT32_COP *outtblptr = NULL; /*to store int32 of outtb*/
200 SCSUINT8_COP *outtbucptr = NULL; /*to store unsigned int8 of outtb */
201 SCSUINT16_COP *outtbusptr = NULL; /*to store unsigned int16 of outtb */
202 SCSUINT32_COP *outtbulptr = NULL; /*to store unsigned int32 of outtb */
204 static outtb_el *outtb_elem = NULL;
206 static scicos_block *Blocks = NULL;
208 /* pass to external variable for code generation */
209 /* reserved variable name */
210 int *block_error = NULL;
211 double scicos_time = 0.;
213 int Jacobian_Flag = 0;
214 // double CI = 0., CJ = 0.; // doubles returned by Get_Jacobian_ci and Get_Jacobian_cj respectively
215 double CJJ = 0.; // returned by Get_Jacobian_parameter
216 double SQuround = 0.;
218 static int AJacobian_block = 0;
221 /* Variable declaration moved to scicos.c because it was in the scicos-def.h therefore
222 * multiple declaration of the variable and linkers were complaining about duplicate
225 SCICOS_IMPEXP COSDEBUGCOUNTER_struct C2F(cosdebugcounter);
226 SCICOS_IMPEXP RTFACTOR_struct C2F(rtfactor);
227 SCICOS_IMPEXP SOLVER_struct C2F(cmsolver);
228 SCICOS_IMPEXP CURBLK_struct C2F(curblk);
229 SCICOS_IMPEXP COSDEBUG_struct C2F(cosdebug);
230 SCICOS_IMPEXP COSHLT_struct C2F(coshlt);
231 SCICOS_IMPEXP DBCOS_struct C2F(dbcos);
232 SCICOS_IMPEXP COSTOL_struct C2F(costol);
233 SCICOS_IMPEXP COSERR_struct coserr;
234 /*--------------------------------------------------------------------------*/
235 static void FREE_blocks();
236 static void cosini(double *told);
237 static void cosend(double *told);
238 static void cossim(double *told);
239 static void cossimdaskr(double *told);
240 static void doit(double *told);
241 static void idoit(double *told);
242 static void zdoit(double *told, double *xt, double *xtd, double *g);
243 static void cdoit(double *told);
244 static void ozdoit(double *told, double *xt, double *xtd, int *kiwa);
245 static void odoit(double *told, double *xt, double *xtd, double *residual);
246 static void ddoit(double *told);
247 static void edoit(double *told, int *kiwa);
248 static void reinitdoit(double *told);
249 static int CallKinsol(double *told);
250 static int simblk(realtype t, N_Vector yy, N_Vector yp, void *f_data);
251 static int grblkdaskr(realtype t, N_Vector yy, N_Vector yp, realtype *gout, void *g_data);
252 static int grblk(realtype t, N_Vector yy, realtype *gout, void *g_data);
253 static void simblklsodar(int * nequations, realtype * tOld, realtype * actual, realtype * res);
254 static void grblklsodar(int * nequations, realtype * tOld, realtype * actual, int * ngc, realtype * res);
255 static void simblkddaskr(realtype *tOld, realtype *actual, realtype *actualP, realtype *res, int *flag, double *dummy1, int *dummy2);
256 static void grblkddaskr(int *nequations, realtype *tOld, realtype *actual, int *ngc, realtype *res, double *dummy1, int *dummy2);
257 static void jacpsol(realtype *res, int *ires, int *nequations, realtype *tOld, realtype *actual, realtype *actualP,
258 realtype *rewt, realtype *savr, realtype *wk, realtype *h, realtype *cj, realtype *wp,
259 int *iwp, int *ier, double *dummy1, int *dummy2);
260 static void psol(int *nequations, realtype *tOld, realtype *actual, realtype *actualP,
261 realtype *savr, realtype *wk, realtype *cj, realtype *wght, realtype *wp,
262 int *iwp, realtype *b, realtype *eplin, int *ier, double *dummy1, int *dummy2);
263 static void addevs(double t, int *evtnb, int *ierr1);
264 static int synchro_g_nev(ScicosImport *scs_imp, double *g, int kf, int *ierr);
265 static void Multp(double *A, double *B, double *R, int ra, int rb, int ca, int cb);
266 static int read_id(ezxml_t *elements, char *id, double *value);
267 static int simblkdaskr(realtype tres, N_Vector yy, N_Vector yp, N_Vector resval, void *rdata);
268 static void SundialsErrHandler(int error_code, const char *module, const char *function, char *msg, void *user_data);
269 static int Jacobians(long int Neq, realtype tt, realtype cj, N_Vector yy,
270 N_Vector yp, N_Vector resvec, DlsMat Jacque, void *jdata,
271 N_Vector tempv1, N_Vector tempv2, N_Vector tempv3);
272 static void call_debug_scicos(scicos_block *block, scicos_flag *flag, int flagi, int deb_blk);
273 static int synchro_nev(ScicosImport *scs_imp, int kf, int *ierr);
274 /*--------------------------------------------------------------------------*/
275 extern int C2F(dset)(int *n, double *dx, double *dy, int *incy);
276 extern int C2F(dcopy)(int *, double *, int *, double *, int *);
277 extern int C2F(dgefa)(double *A, int *lead_dim_A, int *n, int *ipivots, int *info);
278 extern int C2F(dgesl)(double *A, int *lead_dim_A, int *n, int *ipivots, double *B, int *job);
279 extern int C2F(msgs)();
280 extern void C2F(clearscicosimport)();
281 /*--------------------------------------------------------------------------*/
282 void putevs(double *t, int *evtnb, int *ierr1);
283 void Jdoit(double *told, double *xt, double *xtd, double *residual, int *job);
284 int simblkKinsol(N_Vector yy, N_Vector resval, void *rdata);
285 /*--------------------------------------------------------------------------*/
286 int C2F(scicos)(double *x_in, int *xptr_in, double *z__,
287 void **work, int *zptr, int *modptr_in,
288 void **oz, int *ozsz, int *oztyp, int *ozptr,
289 char **iz, int *izptr, char **uid, int *uidptr, double *t0_in,
290 double *tf_in, double *tevts_in, int *evtspt_in,
291 int *nevts, int *pointi_in, void **outtbptr_in,
292 int *outtbsz_in, int *outtbtyp_in,
293 outtb_el *outtb_elem_in, int *nelem1, int *nlnk1,
294 void** funptr, int *funtyp_in, int *inpptr_in,
295 int *outptr_in, int *inplnk_in, int *outlnk_in,
296 double *rpar, int *rpptr, int *ipar, int *ipptr,
297 void **opar, int *oparsz, int *opartyp, int *opptr,
298 int *clkptr_in, int *ordptr_in, int *nordptr1,
299 int *ordclk_in, int *cord_in, int *ncord1,
300 int *iord_in, int *niord1, int *oord_in,
301 int *noord1, int *zord_in, int *nzord1,
302 int *critev_in, int *nblk1, int *ztyp,
303 int *zcptr_in, int *subscr, int *nsubs,
304 double *simpar, int *flag__, int *ierr_out)
306 int i1, kf, lprt, in, out, job = 1;
309 static int mxtb = 0, ierr0 = 0, kfun0 = 0, i = 0, j = 0, k = 0, jj = 0;
310 static int ni = 0, no = 0;
311 static int nz = 0, noz = 0, nopar = 0;
314 // Set FPU Flag to Extended for scicos simulation
315 // in order to override Java setting it to Double.
316 #if defined(linux) && defined(__i386__)
319 /* Copyright INRIA */
320 /* iz,izptr are used to pass block labels */
327 /* Parameter adjustments */
331 modptr = modptr_in - 1;
335 evtspt = evtspt_in - 1;
336 tevts = tevts_in - 1;
337 outtbptr = outtbptr_in;
338 outtbsz = outtbsz_in;
339 outtbtyp = outtbtyp_in;
340 outtb_elem = outtb_elem_in;
341 funtyp = funtyp_in - 1;
342 inpptr = inpptr_in - 1;
343 outptr = outptr_in - 1;
344 inplnk = inplnk_in - 1;
345 outlnk = outlnk_in - 1;
349 clkptr = clkptr_in - 1;
350 ordptr = ordptr_in - 1;
351 ordclk = ordclk_in - 1;
357 critev = critev_in - 1;
359 zcptr = zcptr_in - 1;
367 C2F(rtfactor).scale = simpar[5];
368 C2F(cmsolver).solver = (int) simpar[6];
381 nordclk = ordptr[nordptr] - 1; /* number of rows in ordclk is ordptr(nclkp1)-1 */
382 ng = zcptr[nblk + 1] - 1; /* computes number of zero crossing surfaces */
383 nmod = modptr[nblk + 1] - 1; /* computes number of modes */
384 nz = zptr[nblk + 1] - 1; /* number of discrete real states */
385 noz = ozptr[nblk + 1] - 1; /* number of discrete object states */
386 nopar = opptr[nblk + 1] - 1; /* number of object parameters */
387 nx = xptr[nblk + 1] - 1; /* number of object parameters */
390 xd = &x[xptr[nblk + 1] - 1];
392 /* check for hard coded maxsize */
393 for (i = 1; i <= nblk; ++i)
395 if (funtyp[i] < 10000)
401 funtyp[i] = funtyp[i] % 1000 + 10000;
403 ni = inpptr[i + 1] - inpptr[i];
404 no = outptr[i + 1] - outptr[i];
409 /* hard coded maxsize in callf.c */
410 C2F(msgs)(&c__90, &c__0);
411 C2F(curblk).kfun = i;
416 else if (funtyp[i] == 2 || funtyp[i] == 3)
418 /* hard coded maxsize in scicos.h */
419 if (ni + no > SZ_SIZE)
421 C2F(msgs)(&c__90, &c__0);
422 C2F(curblk).kfun = i;
432 for (j = 1; j <= ni; ++j)
434 k = inplnk[inpptr[i] - 1 + j];
435 mxtb = mxtb + (outtbsz[k - 1] * outtbsz[(k - 1) + nlnk]);
440 for (j = 1; j <= no; ++j)
442 k = outlnk[outptr[i] - 1 + j];
443 mxtb = mxtb + (outtbsz[k - 1] * outtbsz[(k - 1) + nlnk]);
448 C2F(msgs)(&c__91, &c__0);
449 C2F(curblk).kfun = i;
456 if (nx > 0) /* xprop */
458 if ((xprop = MALLOC(sizeof(int) * nx)) == NULL )
464 for (i = 0; i < nx; i++) /* initialize */
468 if (nmod > 0) /* mod */
470 if ((mod = MALLOC(sizeof(int) * nmod)) == NULL )
480 if (ng > 0) /* g becomes global */
482 if ((g = MALLOC(sizeof(double) * ng)) == NULL )
497 debug_block = -1; /* no debug block for start */
498 C2F(cosdebugcounter).counter = 0;
500 /** Create Block's array **/
501 if ((Blocks = MALLOC(sizeof(scicos_block) * nblk)) == NULL )
519 /** Setting blocks properties for each entry in Block's array **/
521 /* 1 : type and pointer on simulation function */
522 for (kf = 0; kf < nblk; ++kf) /*for each block */
524 void* p = funptr[kf];
525 C2F(curblk).kfun = kf + 1;
526 Blocks[kf].type = funtyp[kf + 1];
527 if (Blocks[kf].type < 0)
530 funtyp[kf + 1] *= -1; // Restore a positive 'funtyp' for later use
531 switch (-Blocks[kf].type)
534 Blocks[kf].funpt = (voidg)F2C(sciblk);
537 sciprint(_("type 1 function not allowed for scilab blocks\n"));
538 *ierr = 1000 + kf + 1;
542 sciprint(_("type 2 function not allowed for scilab blocks\n"));
543 *ierr = 1000 + kf + 1;
547 Blocks[kf].funpt = (voidg)sciblk2;
551 Blocks[kf].funpt = (voidg)sciblk4;
554 case 99: /* debugging block */
555 Blocks[kf].funpt = (voidg)sciblk4;
556 /*Blocks[kf].type=4;*/
561 Blocks[kf].funpt = (voidg)sciblk4;
562 Blocks[kf].type = 10004;
565 sciprint(_("Undefined Function type\n"));
566 *ierr = 1000 + kf + 1;
570 Blocks[kf].scsptr = p; /* set scilab function adress for sciblk */
574 //linked functions ( internal or external
575 Blocks[kf].funpt = (voidf)p;
576 Blocks[kf].scsptr = NULL; /* this is done for being able to test if a block
577 is a scilab block in the debugging phase when
582 // switch (funtyp[kf + 1])
585 // Blocks[kf].funpt = (voidg) F2C(sciblk);
588 // sciprint(_("type 1 function not allowed for scilab blocks\n"));
589 // *ierr = 1000 + kf + 1;
593 // sciprint(_("type 2 function not allowed for scilab blocks\n"));
594 // *ierr = 1000 + kf + 1;
598 // Blocks[kf].funpt = (voidg) sciblk2;
599 // Blocks[kf].type = 2;
602 // Blocks[kf].funpt = (voidg) sciblk4;
603 // Blocks[kf].type = 4;
605 // case 99: /* debugging block */
606 // Blocks[kf].funpt = (voidg) sciblk4;
607 // /*Blocks[kf].type=4;*/
612 // Blocks[kf].funpt = (voidg) sciblk4;
613 // Blocks[kf].type = 10004;
616 // sciprint(_("Undefined Function type\n"));
617 // *ierr = 1000 + kf + 1;
621 // Blocks[kf].scsptr = -i; /* set scilab function adress for sciblk */
623 //else if (i <= ntabsim)
625 // Blocks[kf].funpt = (voidg) * (tabsim[i - 1].fonc);
626 // Blocks[kf].scsptr = 0; /* this is done for being able to test if a block
627 // is a scilab block in the debugging phase when
628 // sciblk4 is called */
632 // i -= (ntabsim + 1);
633 // //TODO: see in dynamic_lin how to get funcptr from index
634 // //GetDynFunc(i, &Blocks[kf].funpt);
635 // if ( Blocks[kf].funpt == (voidf) 0)
637 // sciprint(_("Function not found\n"));
638 // *ierr = 1000 + kf + 1;
642 // Blocks[kf].scsptr = 0; /* this is done for being able to test if a block
643 // is a scilab block in the debugging phase when
644 // sciblk4 is called */
647 /* 2 : Dimension properties */
648 Blocks[kf].ztyp = ztyp[kf + 1];
649 Blocks[kf].nx = xptr[kf + 2] - xptr[kf + 1]; /* continuuous state dimension*/
650 Blocks[kf].ng = zcptr[kf + 2] - zcptr[kf + 1]; /* number of zero crossing surface*/
651 Blocks[kf].nz = zptr[kf + 2] - zptr[kf + 1]; /* number of double discrete state*/
652 Blocks[kf].noz = ozptr[kf + 2] - ozptr[kf + 1]; /* number of other discrete state*/
653 Blocks[kf].nrpar = rpptr[kf + 2] - rpptr[kf + 1]; /* size of double precision parameter vector*/
654 Blocks[kf].nipar = ipptr[kf + 2] - ipptr[kf + 1]; /* size of integer precision parameter vector*/
655 Blocks[kf].nopar = opptr[kf + 2] - opptr[kf + 1]; /* number of other parameters (matrix, data structure,..)*/
656 Blocks[kf].nin = inpptr[kf + 2] - inpptr[kf + 1]; /* number of input ports */
657 Blocks[kf].nout = outptr[kf + 2] - outptr[kf + 1]; /* number of output ports */
659 /* 3 : input port properties */
660 /* in insz, we store :
661 * - insz[0..nin-1] : first dimension of input ports
662 * - insz[nin..2*nin-1] : second dimension of input ports
663 * - insz[2*nin..3*nin-1] : type of data of input ports
665 /* allocate size and pointer arrays (number of input ports)*/
666 Blocks[kf].insz = NULL;
667 Blocks[kf].inptr = NULL;
668 if (Blocks[kf].nin != 0)
670 if ((Blocks[kf].insz = MALLOC(Blocks[kf].nin * 3 * sizeof(int))) == NULL )
676 if ((Blocks[kf].inptr = MALLOC(Blocks[kf].nin * sizeof(double*))) == NULL )
683 for (in = 0; in < Blocks[kf].nin; in++)
685 lprt = inplnk[inpptr[kf + 1] + in];
686 Blocks[kf].inptr[in] = outtbptr[lprt - 1]; /* pointer on the data*/
687 Blocks[kf].insz[in] = outtbsz[lprt - 1]; /* row dimension of the input port*/
688 Blocks[kf].insz[Blocks[kf].nin + in] = outtbsz[(lprt - 1) + nlnk]; /* column dimension of the input port*/
689 Blocks[kf].insz[2 * Blocks[kf].nin + in] = outtbtyp[lprt - 1]; /*type of data of the input port*/
691 /* 4 : output port properties */
692 /* in outsz, we store :
693 * - outsz[0..nout-1] : first dimension of output ports
694 * - outsz[nout..2*nout-1] : second dimension of output ports
695 * - outsz[2*nout..3*nout-1] : type of data of output ports
697 /* allocate size and pointer arrays (number of output ports)*/
698 Blocks[kf].outsz = NULL;
699 Blocks[kf].outptr = NULL;
700 if (Blocks[kf].nout != 0)
702 if ((Blocks[kf].outsz = MALLOC(Blocks[kf].nout * 3 * sizeof(int))) == NULL )
708 if ((Blocks[kf].outptr = MALLOC(Blocks[kf].nout * sizeof(double*))) == NULL )
716 for (out = 0; out < Blocks[kf].nout; out++) /*for each output port */
718 lprt = outlnk[outptr[kf + 1] + out];
719 Blocks[kf].outptr[out] = outtbptr[lprt - 1]; /*pointer on data */
720 Blocks[kf].outsz[out] = outtbsz[lprt - 1]; /*row dimension of output port*/
721 Blocks[kf].outsz[Blocks[kf].nout + out] = outtbsz[(lprt - 1) + nlnk]; /*column dimension of output ports*/
722 Blocks[kf].outsz[2 * Blocks[kf].nout + out] = outtbtyp[lprt - 1]; /*type of data of output port */
725 /* 5 : event output port properties */
726 Blocks[kf].evout = NULL;
727 Blocks[kf].nevout = clkptr[kf + 2] - clkptr[kf + 1];
728 if (Blocks[kf].nevout != 0)
730 if ((Blocks[kf].evout = CALLOC(Blocks[kf].nevout, sizeof(double))) == NULL )
738 /* 6 : pointer on the begining of the double discrete state array ( z) */
739 Blocks[kf].z = &(z__[zptr[kf + 1] - 1]);
741 /* 7 : type, size and pointer on the other discrete states data structures (oz) */
742 Blocks[kf].ozsz = NULL;
743 if (Blocks[kf].noz == 0)
745 Blocks[kf].ozptr = NULL;
746 Blocks[kf].oztyp = NULL;
750 Blocks[kf].ozptr = &(oz[ozptr[kf + 1] - 1]);
751 if ((Blocks[kf].ozsz = MALLOC(Blocks[kf].noz * 2 * sizeof(int))) == NULL )
757 for (i = 0; i < Blocks[kf].noz; i++)
759 Blocks[kf].ozsz[i] = ozsz[(ozptr[kf + 1] - 1) + i];
760 Blocks[kf].ozsz[i + Blocks[kf].noz] = ozsz[(ozptr[kf + 1] - 1 + noz) + i];
762 Blocks[kf].oztyp = &(oztyp[ozptr[kf + 1] - 1]);
765 /* 8 : pointer on the begining of the double parameter array ( rpar ) */
766 Blocks[kf].rpar = &(rpar[rpptr[kf + 1] - 1]);
768 /* 9 : pointer on the begining of the integer parameter array ( ipar ) */
769 Blocks[kf].ipar = &(ipar[ipptr[kf + 1] - 1]);
771 /* 10 : type, size and pointer on the other parameters data structures (opar) */
772 Blocks[kf].oparsz = NULL;
773 if (Blocks[kf].nopar == 0)
775 Blocks[kf].oparptr = NULL;
776 Blocks[kf].opartyp = NULL;
780 Blocks[kf].oparptr = &(opar[opptr[kf + 1] - 1]);
781 if ((Blocks[kf].oparsz = MALLOC(Blocks[kf].nopar * 2 * sizeof(int))) == NULL)
787 for (i = 0; i < Blocks[kf].nopar; i++)
789 Blocks[kf].oparsz[i] = oparsz[(opptr[kf + 1] - 1) + i];
790 Blocks[kf].oparsz[i + Blocks[kf].nopar] = oparsz[(opptr[kf + 1] - 1 + nopar) + i];
792 Blocks[kf].opartyp = &(opartyp[opptr[kf + 1] - 1]);
795 /* 10 : pointer on the beginning of the residual array (res) */
796 Blocks[kf].res = NULL;
797 if (Blocks[kf].nx != 0)
799 if ((Blocks[kf].res = MALLOC(Blocks[kf].nx * sizeof(double))) == NULL)
807 /* 11 : block label (label) */
809 if ((Blocks[kf].label = MALLOC(sizeof(char) * (i1 + 1))) == NULL)
815 Blocks[kf].label[i1] = '\0';
816 strcpy(Blocks[kf].label, iz[kf]);
818 /* block uid (uid) */
820 if ((Blocks[kf].uid = MALLOC(sizeof(char) * (i1 + 1))) == NULL)
826 Blocks[kf].uid[i1] = '\0';
827 strcpy(Blocks[kf].uid, uid[kf]);
829 /* 12 : block array of crossed surfaces (jroot) */
830 Blocks[kf].jroot = NULL;
831 if (Blocks[kf].ng > 0)
833 if ((Blocks[kf].jroot = CALLOC(Blocks[kf].ng, sizeof(int))) == NULL)
841 /* 13 : block work array (work) */
842 Blocks[kf].work = (void **)(((double *)work) + kf);
844 /* 14 : block modes array (mode) */
845 Blocks[kf].nmode = modptr[kf + 2] - modptr[kf + 1];
846 if (Blocks[kf].nmode != 0)
848 Blocks[kf].mode = &(mod[modptr[kf + 1] - 1]);
851 /* 15 : block xprop array (xprop) */
852 Blocks[kf].xprop = NULL;
853 if (Blocks[kf].nx != 0)
855 Blocks[kf].xprop = &(xprop[xptr[kf + 1] - 1]);
858 /* 16 : pointer on the zero crossing surface computation function of the block (g) */
860 if (Blocks[kf].ng != 0)
862 Blocks[kf].g = &(g[zcptr[kf + 1] - 1]);
865 /** all block properties are stored in the Blocks array **/
871 if ((iwa = MALLOC(sizeof(int) * (*nevts))) == NULL)
879 /* save ptr of scicos in import structure */
880 makescicosimport(x, &nx, &xptr[1], &zcptr[1], z__, &nz, &zptr[1],
881 &noz, oz, ozsz, oztyp, &ozptr[1],
882 g, &ng, mod, &nmod, &modptr[1], iz, &izptr[1],
884 &inpptr[1], &inplnk[1], &outptr[1], &outlnk[1],
885 outtbptr, outtbsz, outtbtyp,
887 &nlnk, rpar, &rpptr[1], ipar, &ipptr[1],
888 opar, oparsz, opartyp, &opptr[1],
889 &nblk, subscr, nsubs,
890 &tevts[1], &evtspt[1], nevts, pointi,
891 &iord[1], &niord, &oord[1], &noord, &zord[1], &nzord,
892 funptr, &funtyp[1], &ztyp[1],
893 &cord[1], &ncord, &ordclk[1], &nordclk, &clkptr[1],
894 &ordptr[1], &nordptr, &critev[1], iwa, Blocks,
895 t0, tf, &Atol, &rtol, &ttol, &deltat, &hmax,
898 if (*flag__ == 1) /*start*/
900 /* blocks initialization */
901 for (kf = 0; kf < nblk; ++kf)
903 *(Blocks[kf].work) = NULL;
909 kfun0 = C2F(curblk).kfun;
912 C2F(curblk).kfun = kfun0;
916 else if (*flag__ == 2) /*run*/
920 switch (C2F(cmsolver).solver)
923 case CVode_BDF_Newton:
924 case CVode_BDF_Functional:
925 case CVode_Adams_Newton:
926 case CVode_Adams_Functional:
929 case Implicit_Runge_Kutta:
933 case DDaskr_BDF_Newton:
934 case DDaskr_BDF_GMRes:
937 default: // Unknown solver number
944 kfun0 = C2F(curblk).kfun;
947 C2F(curblk).kfun = kfun0;
951 else if (*flag__ == 3) /*finish*/
956 else if (*flag__ == 4) /*linear*/
962 if ((W = MALLOC(sizeof(double) * (Max(nx, ng)))) == NULL )
970 /*---------instead of old simblk--------*/
971 /* C2F(simblk)(&nx, t0, x, W); */
973 if (ng > 0 && nmod > 0)
975 zdoit(t0, x, x + nx, W); /* updating modes as a function of state values; this was necessary in iGUI*/
977 for (jj = 0; jj < nx; jj++)
981 C2F(ierode).iero = 0;
983 if (C2F(cmsolver).solver < 100)
989 odoit(t0, x, x + nx, W);
991 C2F(ierode).iero = *ierr;
992 /*-----------------------------------------*/
993 for (i = 0; i < nx; ++i)
1000 else if (*flag__ == 5) /* initial_KINSOL= "Kinsol" */
1002 C2F(ierode).iero = 0;
1006 *ierr = C2F(ierode).iero;
1013 C2F(clearscicosimport)();
1016 /*--------------------------------------------------------------------------*/
1018 static int check_flag(void *flagvalue, char *funcname, int opt)
1020 int *errflag = NULL;
1022 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
1023 if (opt == 0 && flagvalue == NULL)
1025 sciprint(_("\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n"), funcname);
1028 /* Check if flag < 0 */
1031 errflag = (int *) flagvalue;
1034 sciprint(_("\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n"),
1035 funcname, *errflag);
1039 /* Check if function returned NULL pointer - no memory allocated */
1040 else if (opt == 2 && flagvalue == NULL)
1042 sciprint(_("\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n"), funcname);
1049 /*--------------------------------------------------------------------------*/
1050 static void cosini(double *told)
1052 static scicos_flag flag__ = 0;
1055 static int kfune = 0;
1058 SCSREAL_COP *outtbd = NULL; /*to save double of outtb*/
1059 SCSINT8_COP *outtbc = NULL; /*to save int8 of outtb*/
1060 SCSINT16_COP *outtbs = NULL; /*to save int16 of outtb*/
1061 SCSINT32_COP *outtbl = NULL; /*to save int32 of outtb*/
1062 SCSUINT8_COP *outtbuc = NULL; /*to save unsigned int8 of outtb*/
1063 SCSUINT16_COP *outtbus = NULL; /*to save unsigned int16 of outtb*/
1064 SCSUINT32_COP *outtbul = NULL; /*to save unsigned int32 of outtb*/
1065 int szouttbd = 0; /*size of arrays*/
1066 int szouttbc = 0, szouttbs = 0, szouttbl = 0;
1067 int szouttbuc = 0, szouttbus = 0, szouttbul = 0;
1068 int curouttbd = 0; /*current position in arrays*/
1069 int curouttbc = 0, curouttbs = 0, curouttbl = 0;
1070 int curouttbuc = 0, curouttbus = 0, curouttbul = 0;
1072 int ii = 0, kk = 0; /*local counters*/
1073 int sszz = 0; /*local size of element of outtb*/
1074 /*Allocation of arrays for outtb*/
1075 for (ii = 0; ii < nlnk; ii++)
1077 switch (outtbtyp[ii])
1080 szouttbd += outtbsz[ii] * outtbsz[ii + nlnk]; /*double real matrix*/
1081 outtbd = (SCSREAL_COP *) REALLOC (outtbd, szouttbd * sizeof(SCSREAL_COP));
1082 memset(outtbd, 0, szouttbd * sizeof(SCSREAL_COP));
1086 szouttbd += 2 * outtbsz[ii] * outtbsz[ii + nlnk]; /*double complex matrix*/
1087 outtbd = (SCSCOMPLEX_COP *) REALLOC (outtbd, szouttbd * sizeof(SCSCOMPLEX_COP));
1088 memset(outtbd, 0, szouttbd * sizeof(SCSCOMPLEX_COP));
1092 szouttbc += outtbsz[ii] * outtbsz[ii + nlnk]; /*int8*/
1093 outtbc = (SCSINT8_COP *) REALLOC (outtbc, szouttbc * sizeof(SCSINT8_COP));
1094 memset(outtbc, 0, szouttbc * sizeof(SCSINT8_COP));
1098 szouttbs += outtbsz[ii] * outtbsz[ii + nlnk]; /*int16*/
1099 outtbs = (SCSINT16_COP *) REALLOC (outtbs, szouttbs * sizeof(SCSINT16_COP));
1100 memset(outtbs, 0, szouttbs * sizeof(SCSINT16_COP));
1104 szouttbl += outtbsz[ii] * outtbsz[ii + nlnk]; /*int32*/
1105 outtbl = (SCSINT32_COP *) REALLOC (outtbl, szouttbl * sizeof(SCSINT32_COP));
1106 memset(outtbl, 0, szouttbl * sizeof(SCSINT32_COP));
1110 szouttbuc += outtbsz[ii] * outtbsz[ii + nlnk]; /*uint8*/
1111 outtbuc = (SCSUINT8_COP *) REALLOC (outtbuc, szouttbuc * sizeof(SCSUINT8_COP));
1112 memset(outtbuc, 0, szouttbuc * sizeof(SCSUINT8_COP));
1116 szouttbus += outtbsz[ii] * outtbsz[ii + nlnk]; /*uint16*/
1117 outtbus = (SCSUINT16_COP *) REALLOC (outtbus, szouttbus * sizeof(SCSUINT16_COP));
1118 memset(outtbus, 0, szouttbus * sizeof(SCSUINT16_COP));
1122 szouttbul += outtbsz[ii] * outtbsz[ii + nlnk]; /*uint32*/
1123 outtbul = (SCSUINT32_COP *) REALLOC (outtbul, szouttbul * sizeof(SCSUINT32_COP));
1124 memset(outtbul, 0, szouttbul * sizeof(SCSUINT32_COP));
1127 default : /* Add a message here */
1133 AJacobian_block = 0;
1138 /* initialization (flag 4) */
1139 /* loop on blocks */
1140 C2F(dset)(&ng, &c_b14, g, &c__1);
1142 for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
1145 if (Blocks[C2F(curblk).kfun - 1].nx > 0)
1147 Blocks[C2F(curblk).kfun - 1].x = &x[xptr[C2F(curblk).kfun] - 1];
1148 Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
1150 Blocks[C2F(curblk).kfun - 1].nevprt = 0;
1151 if (funtyp[C2F(curblk).kfun] >= 0) /* debug_block is not called here */
1153 /*callf(told, xd, x, x,g,&flag__);*/
1155 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1156 if (flag__ < 0 && *ierr == 0)
1159 kfune = C2F(curblk).kfun;
1161 if ((Jacobian_Flag == 1) && (AJacobian_block == 0))
1163 AJacobian_block = C2F(curblk).kfun;
1169 C2F(curblk).kfun = kfune;
1174 /* initialization (flag 6) */
1176 for (jj = 1; jj <= ncord; ++jj)
1178 C2F(curblk).kfun = cord[jj];
1179 Blocks[C2F(curblk).kfun - 1].nevprt = 0;
1180 if (funtyp[C2F(curblk).kfun] >= 0)
1182 /*callf(told, xd, x, x,g,&flag__);*/
1183 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1193 /* point-fix iterations */
1195 for (i = 1; i <= nblk + 1; ++i) /*for each block*/
1197 /* loop on blocks */
1198 for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
1200 Blocks[C2F(curblk).kfun - 1].nevprt = 0;
1201 if (funtyp[C2F(curblk).kfun] >= 0)
1203 /*callf(told, xd, x, x,g,&flag__);*/
1204 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1215 for (jj = 1; jj <= ncord; ++jj) /*for each continous block*/
1217 C2F(curblk).kfun = cord[jj];
1218 if (funtyp[C2F(curblk).kfun] >= 0)
1220 /*callf(told, xd, x, x,g,&flag__);*/
1221 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1231 /*comparison between outtb and arrays*/
1239 for (jj = 0; jj < nlnk; jj++)
1241 switch (outtbtyp[jj]) /*for each type of ports*/
1244 outtbdptr = (SCSREAL_COP *)outtbptr[jj]; /*double real matrix*/
1245 sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1246 for (kk = 0; kk < sszz; kk++)
1248 int outtbdptr_isnan = outtbdptr[kk] != outtbdptr[kk];
1249 int outtbd_isnan = (SCSREAL_COP)outtbd[curouttbd + kk] != (SCSREAL_COP)outtbd[curouttbd + kk];
1251 if (outtbdptr_isnan && outtbd_isnan)
1255 if (outtbdptr[kk] != (SCSREAL_COP)outtbd[curouttbd + kk])
1264 outtbdptr = (SCSCOMPLEX_COP *)outtbptr[jj]; /*double complex matrix*/
1265 sszz = 2 * outtbsz[jj] * outtbsz[jj + nlnk];
1266 for (kk = 0; kk < sszz; kk++)
1268 int outtbdptr_isnan = outtbdptr[kk] != outtbdptr[kk];
1269 int outtbd_isnan = (SCSCOMPLEX_COP)outtbd[curouttbd + kk] != (SCSCOMPLEX_COP)outtbd[curouttbd + kk];
1271 if (outtbdptr_isnan && outtbd_isnan)
1275 if (outtbdptr[kk] != (SCSCOMPLEX_COP)outtbd[curouttbd + kk])
1284 outtbcptr = (SCSINT8_COP *)outtbptr[jj]; /*int8*/
1285 sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1286 for (kk = 0; kk < sszz; kk++)
1288 if (outtbcptr[kk] != (SCSINT8_COP)outtbc[curouttbc + kk])
1297 outtbsptr = (SCSINT16_COP *)outtbptr[jj]; /*int16*/
1298 sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1299 for (kk = 0; kk < sszz; kk++)
1301 if (outtbsptr[kk] != (SCSINT16_COP)outtbs[curouttbs + kk])
1310 outtblptr = (SCSINT32_COP *)outtbptr[jj]; /*int32*/
1311 sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1312 for (kk = 0; kk < sszz; kk++)
1314 if (outtblptr[kk] != (SCSINT32_COP)outtbl[curouttbl + kk])
1323 outtbucptr = (SCSUINT8_COP *)outtbptr[jj]; /*uint8*/
1324 sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1325 for (kk = 0; kk < sszz; kk++)
1327 if (outtbucptr[kk] != (SCSUINT8_COP)outtbuc[curouttbuc + kk])
1336 outtbusptr = (SCSUINT16_COP *)outtbptr[jj]; /*uint16*/
1337 sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1338 for (kk = 0; kk < sszz; kk++)
1340 if (outtbusptr[kk] != (SCSUINT16_COP)outtbus[curouttbus + kk])
1349 outtbulptr = (SCSUINT32_COP *)outtbptr[jj]; /*uint32*/
1350 sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1351 for (kk = 0; kk < sszz; kk++)
1353 if (outtbulptr[kk] != (SCSUINT32_COP)outtbul[curouttbul + kk])
1361 default : /* Add a message here */
1369 /*Save data of outtb in arrays*/
1377 for (ii = 0; ii < nlnk; ii++) /*for each link*/
1379 switch (outtbtyp[ii]) /*switch to type of outtb object*/
1382 outtbdptr = (SCSREAL_COP *)outtbptr[ii]; /*double real matrix*/
1383 sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1384 C2F(dcopy)(&sszz, outtbdptr, &c__1, &outtbd[curouttbd], &c__1);
1389 outtbdptr = (SCSCOMPLEX_COP *)outtbptr[ii]; /*double complex matrix*/
1390 sszz = 2 * outtbsz[ii] * outtbsz[ii + nlnk];
1391 C2F(dcopy)(&sszz, outtbdptr, &c__1, &outtbd[curouttbd], &c__1);
1396 outtbcptr = (SCSINT8_COP *)outtbptr[ii]; /*int8*/
1397 sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1398 for (kk = 0; kk < sszz; kk++)
1400 outtbc[curouttbc + kk] = (SCSINT8_COP)outtbcptr[kk];
1406 outtbsptr = (SCSINT16_COP *)outtbptr[ii]; /*int16*/
1407 sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1408 for (kk = 0; kk < sszz; kk++)
1410 outtbs[curouttbs + kk] = (SCSINT16_COP)outtbsptr[kk];
1416 outtblptr = (SCSINT32_COP *)outtbptr[ii]; /*int32*/
1417 sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1418 for (kk = 0; kk < sszz; kk++)
1420 outtbl[curouttbl + kk] = (SCSINT32_COP)outtblptr[kk];
1426 outtbucptr = (SCSUINT8_COP *)outtbptr[ii]; /*uint8*/
1427 sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1428 for (kk = 0; kk < sszz; kk++)
1430 outtbuc[curouttbuc + kk] = (SCSUINT8_COP)outtbucptr[kk];
1436 outtbusptr = (SCSUINT16_COP *)outtbptr[ii]; /*uint16*/
1437 sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1438 for (kk = 0; kk < sszz; kk++)
1440 outtbus[curouttbus + kk] = (SCSUINT16_COP)outtbusptr[kk];
1446 outtbulptr = (SCSUINT32_COP *)outtbptr[ii]; /*uint32*/
1447 sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1448 for (kk = 0; kk < sszz; kk++)
1450 outtbul[curouttbul + kk] = (SCSUINT32_COP)outtbulptr[kk];
1455 default : /* Add a message here */
1464 /*--------------------------------------------------------------------------*/
1465 static void cossim(double *told)
1467 /* System generated locals */
1470 //** used for the [stop] button
1471 static char CommandToUnstack[1024];
1472 static int CommandLength = 0;
1473 static int SeqSync = 0;
1476 /* Local variables */
1477 static scicos_flag flag__ = 0;
1478 static int ierr1 = 0;
1479 static int j = 0, k = 0;
1480 static double t = 0.;
1482 static double rhotmp = 0., tstop = 0.;
1483 static int inxsci = 0;
1484 static int kpo = 0, kev = 0;
1485 int Discrete_Jump = 0;
1486 int *jroot = NULL, *zcros = NULL;
1487 realtype reltol = 0., abstol = 0.;
1489 void *ode_mem = NULL;
1490 int flag = 0, flagr = 0;
1492 /* Saving solver number */
1493 int solver = C2F(cmsolver).solver;
1494 /* Defining function pointers, for more readability */
1495 void(* ODEFree) (void**);
1496 int (* ODE) (void*, realtype, N_Vector, realtype*, int);
1497 int (* ODEReInit) (void*, realtype, N_Vector);
1498 int (* ODESetMaxStep) (void*, realtype);
1499 int (* ODESetStopTime) (void*, realtype);
1500 int (* ODEGetRootInfo) (void*, int*);
1501 int (* ODESStolerances) (void*, realtype, realtype);
1502 /* Generic flags for stop mode */
1503 int ODE_NORMAL = 1; /* ODE_NORMAL = CV_NORMAL = LS_NORMAL = 1 */
1504 int ODE_ONE_STEP = 2; /* ODE_ONE_STEP = CV_ONE_STEP = LS_ONE_STEP = 2 */
1507 case LSodar_Dynamic:
1508 ODEFree = &LSodarFree;
1510 ODEReInit = &LSodarReInit;
1511 ODESetMaxStep = &LSodarSetMaxStep;
1512 ODESetStopTime = &LSodarSetStopTime;
1513 ODEGetRootInfo = &LSodarGetRootInfo;
1514 ODESStolerances = &LSodarSStolerances;
1516 case CVode_BDF_Newton:
1517 case CVode_BDF_Functional:
1518 case CVode_Adams_Newton:
1519 case CVode_Adams_Functional:
1520 case Dormand_Prince:
1522 case Implicit_Runge_Kutta:
1523 ODEFree = &CVodeFree;
1525 ODEReInit = &CVodeReInit;
1526 ODESetMaxStep = &CVodeSetMaxStep;
1527 ODESetStopTime = &CVodeSetStopTime;
1528 ODEGetRootInfo = &CVodeGetRootInfo;
1529 ODESStolerances = &CVodeSStolerances;
1531 default: // Unknown solver number
1539 if ((jroot = MALLOC(sizeof(int) * ng)) == NULL )
1546 for ( jj = 0 ; jj < ng ; jj++ )
1554 if ((zcros = MALLOC(sizeof(int) * ng)) == NULL )
1565 reltol = (realtype) rtol;
1566 abstol = (realtype) Atol; /* Ith(abstol,1) = realtype) Atol;*/
1568 if (*neq > 0) /* Unfortunately CVODE does not work with NEQ==0 */
1570 y = N_VNewEmpty_Serial(*neq);
1571 if (check_flag((void *)y, "N_VNewEmpty_Serial", 0))
1589 /* Set extension of Sundials for scicos */
1590 set_sundials_with_extension(TRUE);
1594 case LSodar_Dynamic:
1595 ode_mem = LSodarCreate(neq, ng); /* Create the lsodar problem */
1597 case CVode_BDF_Newton:
1598 ode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
1600 case CVode_BDF_Functional:
1601 ode_mem = CVodeCreate(CV_BDF, CV_FUNCTIONAL);
1603 case CVode_Adams_Newton:
1604 ode_mem = CVodeCreate(CV_ADAMS, CV_NEWTON);
1606 case CVode_Adams_Functional:
1607 ode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);
1609 case Dormand_Prince:
1610 ode_mem = CVodeCreate(CV_DOPRI, CV_FUNCTIONAL);
1613 ode_mem = CVodeCreate(CV_ExpRK, CV_FUNCTIONAL);
1615 case Implicit_Runge_Kutta:
1616 ode_mem = CVodeCreate(CV_ImpRK, CV_FUNCTIONAL);
1620 /* ode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);*/
1622 if (check_flag((void *)ode_mem, "CVodeCreate", 0))
1625 N_VDestroy_Serial(y);
1631 if (solver == LSodar_Dynamic)
1633 flag = LSodarSetErrHandlerFn(ode_mem, SundialsErrHandler, NULL);
1637 flag = CVodeSetErrHandlerFn(ode_mem, SundialsErrHandler, NULL);
1639 if (check_flag(&flag, "CVodeSetErrHandlerFn", 1))
1641 *ierr = 300 + (-flag);
1646 if (solver == LSodar_Dynamic)
1648 flag = LSodarInit(ode_mem, simblklsodar, T0, y);
1652 flag = CVodeInit (ode_mem, simblk, T0, y);
1654 if (check_flag(&flag, "CVodeInit", 1))
1656 *ierr = 300 + (-flag);
1661 flag = ODESStolerances(ode_mem, reltol, abstol);
1662 if (check_flag(&flag, "CVodeSStolerances", 1))
1664 *ierr = 300 + (-flag);
1669 if (solver == LSodar_Dynamic)
1671 flag = LSodarRootInit(ode_mem, ng, grblklsodar);
1675 flag = CVodeRootInit(ode_mem, ng, grblk);
1677 if (check_flag(&flag, "CVodeRootInit", 1))
1679 *ierr = 300 + (-flag);
1684 if (solver != LSodar_Dynamic) /* Call CVDense to specify the CVDENSE dense linear solver */
1686 flag = CVDense(ode_mem, *neq);
1688 if (check_flag(&flag, "CVDense", 1))
1690 *ierr = 300 + (-flag);
1697 flag = ODESetMaxStep(ode_mem, (realtype) hmax);
1698 if (check_flag(&flag, "CVodeSetMaxStep", 1))
1700 *ierr = 300 + (-flag);
1705 /* Set the Jacobian routine to Jac (user-supplied)
1706 flag = CVDlsSetDenseJacFn(ode_mem, Jac);
1707 if (check_flag(&flag, "CVDlsSetDenseJacFn", 1)) return(1); */
1709 }/* testing if neq>0 */
1712 C2F(coshlt).halt = 0;
1715 C2F(xscion)(&inxsci);
1716 /* initialization */
1717 C2F(realtimeinit)(told, &C2F(rtfactor).scale);
1723 for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
1725 if (Blocks[C2F(curblk).kfun - 1].ng >= 1)
1727 zcros[jj] = C2F(curblk).kfun;
1731 /* . ng >= jj required */
1736 /* initialization (propagation of constant blocks outputs) */
1743 /*--discrete zero crossings----dzero--------------------*/
1744 if (ng > 0) /* storing ZC signs just after a solver call*/
1746 /*zdoit(told, g, x, x);*/
1747 zdoit(told, x, x, g);
1753 for (jj = 0; jj < ng; ++jj)
1763 /*--discrete zero crossings----dzero--------------------*/
1765 /* main loop on time */
1768 while (ismenu()) //** if the user has done something, do the actions
1771 SeqSync = GetCommand(CommandToUnstack); //** get to the action
1772 CommandLength = (int)strlen(CommandToUnstack);
1773 //syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
1775 if (C2F(coshlt).halt != 0)
1777 if (C2F(coshlt).halt == 2)
1779 *told = *tf; /* end simulation */
1781 C2F(coshlt).halt = 0;
1793 if (fabs(t - *told) < ttol)
1796 /* update output part */
1800 /* ! scheduling problem */
1807 if (xptr[nblk + 1] == 1)
1809 /* . no continuous state */
1810 if (*told + deltat + ttol > t)
1818 /* . update outputs of 'c' type blocks with no continuous state */
1821 /* . we are at the end, update continuous part before leaving */
1833 rhotmp = *tf + ttol;
1838 if (critev[kpo] == 1)
1840 rhotmp = tevts[kpo];
1855 t = Min(*told + deltat, Min(t, *tf + ttol));
1857 if (ng > 0 && hot == 0 && nmod > 0)
1859 zdoit(told, x, x, g);
1867 if (hot == 0) /* hot==0 : cold restart*/
1869 flag = ODESetStopTime(ode_mem, (realtype)tstop); /* Setting the stop time*/
1870 if (check_flag(&flag, "CVodeSetStopTime", 1))
1872 *ierr = 300 + (-flag);
1877 flag = ODEReInit(ode_mem, (realtype)(*told), y);
1878 if (check_flag(&flag, "CVodeReInit", 1))
1880 *ierr = 300 + (-flag);
1886 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1888 sciprint(_("****SUNDIALS.Cvode from: %f to %f hot= %d \n"), *told, t, hot);
1891 /*--discrete zero crossings----dzero--------------------*/
1892 /*--check for Dzeros after Mode settings or ddoit()----*/
1895 if (ng > 0 && hot == 0)
1897 zdoit(told, x, x, g);
1903 for (jj = 0; jj < ng; ++jj)
1905 if ((g[jj] >= 0.0) && (jroot[jj] == -5))
1910 else if ((g[jj] < 0.0) && (jroot[jj] == 5))
1921 /*--discrete zero crossings----dzero--------------------*/
1923 if (Discrete_Jump == 0) /* if there was a dzero, its event should be activated*/
1926 flag = ODE(ode_mem, t, y, told, ODE_NORMAL);
1936 flag = CV_ROOT_RETURN; /* in order to handle discrete jumps */
1939 /* . update outputs of 'c' type blocks if we are at the end*/
1952 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1954 sciprint(_("****SUNDIALS.Cvode reached: %f\n"), *told);
1959 else if ( flag == CV_TOO_MUCH_WORK || flag == CV_CONV_FAILURE || flag == CV_ERR_FAILURE)
1961 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1963 sciprint(_("****SUNDIALS.Cvode: too much work at time=%g (stiff region, change RTOL and ATOL)\n"), *told);
1969 *ierr = 300 + (-flag);
1978 *ierr = 300 + (-flag); /* raising errors due to internal errors, otherwise error due to flagr*/
1984 if (flag == CV_ZERO_DETACH_RETURN)
1987 }; /* new feature of sundials, detects zero-detaching */
1989 if (flag == CV_ROOT_RETURN)
1991 /* . at a least one root has been found */
1993 if (Discrete_Jump == 0)
1995 flagr = ODEGetRootInfo(ode_mem, jroot);
1996 if (check_flag(&flagr, "CVodeGetRootInfo", 1))
1998 *ierr = 300 + (-flagr);
2003 /* . at a least one root has been found */
2004 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2006 sciprint(_("root found at t=: %f\n"), *told);
2008 /* . update outputs affecting ztyp blocks ONLY FOR OLD BLOCKS */
2009 zdoit(told, x, xd, g);
2015 for (jj = 0; jj < ng; ++jj)
2017 C2F(curblk).kfun = zcros[ jj];
2018 if (C2F(curblk).kfun == -1)
2024 for (j = zcptr[C2F(curblk).kfun] - 1 ;
2025 j < zcptr[C2F(curblk).kfun + 1] - 1 ; ++j)
2036 Blocks[C2F(curblk).kfun - 1].jroot = &jroot[zcptr[C2F(curblk).kfun] - 1];
2037 if (funtyp[C2F(curblk).kfun] > 0)
2040 if (Blocks[C2F(curblk).kfun - 1].nevout > 0)
2043 if (Blocks[C2F(curblk).kfun - 1].nx > 0)
2045 Blocks[C2F(curblk).kfun - 1].x = &x[xptr[C2F(curblk).kfun] - 1];
2046 Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
2048 /* call corresponding block to determine output event (kev) */
2049 Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
2050 /*callf(told, xd, x, x,g,&flag__);*/
2051 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
2058 /* . update event agenda */
2059 for (k = 0; k < Blocks[C2F(curblk).kfun - 1].nevout; ++k)
2061 if (Blocks[C2F(curblk).kfun - 1].evout[k] >= 0.)
2063 i3 = k + clkptr[C2F(curblk).kfun] ;
2064 addevs(Blocks[C2F(curblk).kfun - 1].evout[k] + (*told), &i3, &ierr1);
2067 /* . nevts too small */
2075 /* . update state */
2076 if (Blocks[C2F(curblk).kfun - 1].nx > 0)
2078 /* . call corresponding block to update state */
2080 Blocks[C2F(curblk).kfun - 1].x = &x[xptr[C2F(curblk).kfun] - 1];
2081 Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
2082 Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
2083 /*callf(told, xd, x, x,g,&flag__);*/
2084 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
2098 /*--discrete zero crossings----dzero--------------------*/
2099 if (ng > 0) /* storing ZC signs just after a sundials call*/
2101 zdoit(told, x, x, g);
2107 for (jj = 0; jj < ng; ++jj)
2119 /*--discrete zero crossings----dzero--------------------*/
2121 C2F(realtime)(told);
2126 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2128 sciprint(_("Event: %d activated at t=%f\n"), *pointi, *told);
2129 for (kev = 0; kev < nblk; kev++)
2131 if (Blocks[kev].nmode > 0)
2133 sciprint(_("mode of block %d=%d, "), kev, Blocks[kev].mode[0]);
2136 sciprint(_("**mod**\n"));
2140 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2142 sciprint(_("End of activation\n"));
2151 /* end of main loop on time */
2156 /*--------------------------------------------------------------------------*/
2157 static void cossimdaskr(double *told)
2159 /* System generated locals */
2161 //** used for the [stop] button
2162 static char CommandToUnstack[1024];
2163 static int CommandLength = 0;
2164 static int SeqSync = 0;
2167 /* Local variables */
2168 static scicos_flag flag__ = 0;
2169 static int ierr1 = 0;
2170 static int j = 0, k = 0;
2171 static double t = 0.;
2173 static double rhotmp = 0., tstop = 0.;
2174 static int inxsci = 0;
2175 static int kpo = 0, kev = 0;
2177 int *jroot = NULL, *zcros = NULL;
2178 int *Mode_save = NULL;
2179 int Mode_change = 0;
2181 int flag = 0, flagr = 0;
2182 N_Vector yy = NULL, yp = NULL;
2183 realtype reltol = 0., abstol = 0.;
2184 int Discrete_Jump = 0;
2185 N_Vector IDx = NULL;
2186 realtype *scicos_xproperty = NULL;
2187 DlsMat TJacque = NULL;
2189 void *dae_mem = NULL;
2190 UserData data = NULL;
2191 IDAMem copy_IDA_mem = NULL;
2192 int maxnj = 0, maxnit = 0, maxnh = 0;
2193 /*-------------------- Analytical Jacobian memory allocation ----------*/
2194 int Jn = 0, Jnx = 0, Jno = 0, Jni = 0, Jactaille = 0;
2196 int cnt = 0, N_iters = 0;
2197 /* Saving solver number */
2198 int solver = C2F(cmsolver).solver;
2199 /* Flags for initial values calculation */
2200 int DAE_YA_YDP_INIT = 1;
2202 /* Defining function pointers, for more readability*/
2203 void(* DAEFree) (void**);
2204 int (* DAESolve) (void*, realtype, realtype*, N_Vector, N_Vector, int);
2205 int (* DAEReInit) (void*, realtype, N_Vector, N_Vector);
2206 int (* DAESetId) (void*, N_Vector);
2207 int (* DAECalcIC) (void*, int, realtype);
2208 int (* DAESetMaxStep) (void*, realtype);
2209 int (* DAESetUserData) (void*, void*);
2210 int (* DAESetStopTime) (void*, realtype);
2211 int (* DAEGetRootInfo) (void*, int*);
2212 int (* DAESStolerances) (void*, realtype, realtype);
2213 int (* DAEGetConsistentIC) (void*, N_Vector, N_Vector);
2214 int (* DAESetMaxNumSteps) (void*, long int);
2215 int (* DAESetMaxNumJacsIC) (void*, int);
2216 int (* DAESetMaxNumItersIC) (void*, int);
2217 int (* DAESetMaxNumStepsIC) (void*, int);
2218 int (* DAESetLineSearchOffIC) (void*, int);
2219 /* For DAEs, the generic flags for stop mode depend on the used solver */
2220 int DAE_NORMAL = 0, DAE_ONE_STEP = 0;
2221 DAE_NORMAL = (solver == IDA_BDF_Newton) ? 1 : 0; /* IDA_NORMAL = 1, DDAS_NORMAL = 0 */
2222 DAE_ONE_STEP = (solver == IDA_BDF_Newton) ? 2 : 1; /* IDA_ONE_STEP = 2, DDAS_ONE_STEP = 1 */
2225 case IDA_BDF_Newton:
2227 DAESolve = &IDASolve;
2228 DAESetId = &IDASetId;
2229 DAEReInit = &IDAReInit;
2230 DAECalcIC = &IDACalcIC;
2231 DAESetMaxStep = &IDASetMaxStep;
2232 DAESetUserData = &IDASetUserData;
2233 DAESetStopTime = &IDASetStopTime;
2234 DAEGetRootInfo = &IDAGetRootInfo;
2235 DAESStolerances = &IDASStolerances;
2236 DAESetMaxNumSteps = &IDASetMaxNumSteps;
2237 DAEGetConsistentIC = &IDAGetConsistentIC;
2238 DAESetMaxNumJacsIC = &IDASetMaxNumJacsIC;
2239 DAESetMaxNumItersIC = &IDASetMaxNumItersIC;
2240 DAESetMaxNumStepsIC = &IDASetMaxNumStepsIC;
2241 DAESetLineSearchOffIC = &IDASetLineSearchOffIC;
2243 case DDaskr_BDF_Newton:
2244 case DDaskr_BDF_GMRes:
2245 DAEFree = &DDaskrFree;
2246 DAESolve = &DDaskrSolve;
2247 DAESetId = &DDaskrSetId;
2248 DAEReInit = &DDaskrReInit;
2249 DAECalcIC = &DDaskrCalcIC;
2250 DAESetMaxStep = &DDaskrSetMaxStep;
2251 DAESetUserData = &DDaskrSetUserData;
2252 DAESetStopTime = &DDaskrSetStopTime;
2253 DAEGetRootInfo = &DDaskrGetRootInfo;
2254 DAESStolerances = &DDaskrSStolerances;
2255 DAESetMaxNumSteps = &DDaskrSetMaxNumSteps;
2256 DAEGetConsistentIC = &DDaskrGetConsistentIC;
2257 DAESetMaxNumJacsIC = &DDaskrSetMaxNumJacsIC;
2258 DAESetMaxNumItersIC = &DDaskrSetMaxNumItersIC;
2259 DAESetMaxNumStepsIC = &DDaskrSetMaxNumStepsIC;
2260 DAESetLineSearchOffIC = &DDaskrSetLineSearchOffIC;
2262 default: // Unknown solver number
2267 /* Set extension of Sundials for scicos */
2268 set_sundials_with_extension(TRUE);
2270 // CI=1.0; /* for function Get_Jacobian_ci */
2274 if ((jroot = MALLOC(sizeof(int) * ng)) == NULL )
2280 for ( jj = 0 ; jj < ng ; jj++ )
2288 if ((zcros = MALLOC(sizeof(int) * ng)) == NULL )
2302 if ((Mode_save = MALLOC(sizeof(int) * nmod)) == NULL )
2317 reltol = (realtype) rtol;
2318 abstol = (realtype) Atol; /* Ith(abstol,1) = (realtype) Atol;*/
2323 yy = N_VNewEmpty_Serial(*neq);
2324 if (check_flag((void *)yy, "N_VNew_Serial", 0))
2342 yp = N_VNewEmpty_Serial(*neq);
2343 if (check_flag((void *)yp, "N_VNew_Serial", 0))
2347 N_VDestroy_Serial(yy);
2366 IDx = N_VNew_Serial(*neq);
2367 if (check_flag((void *)IDx, "N_VNew_Serial", 0))
2372 N_VDestroy_Serial(yp);
2376 N_VDestroy_Serial(yy);
2393 /* Call the Create and Init functions to initialize DAE memory */
2395 if (solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes)
2397 dae_mem = DDaskrCreate(neq, ng, solver);
2401 dae_mem = IDACreate();
2403 if (check_flag((void *)dae_mem, "IDACreate", 0))
2407 N_VDestroy_Serial(IDx);
2411 N_VDestroy_Serial(yp);
2415 N_VDestroy_Serial(yy);
2431 if (C2F(cmsolver).solver == 100)
2433 copy_IDA_mem = (IDAMem) dae_mem;
2436 if (solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes)
2438 flag = DDaskrSetErrHandlerFn(dae_mem, SundialsErrHandler, NULL);
2442 flag = IDASetErrHandlerFn(dae_mem, SundialsErrHandler, NULL);
2444 if (check_flag(&flag, "IDASetErrHandlerFn", 1))
2446 *ierr = 200 + (-flag);
2453 N_VDestroy_Serial(IDx);
2457 N_VDestroy_Serial(yp);
2461 N_VDestroy_Serial(yy);
2478 if (solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes)
2480 flag = DDaskrInit(dae_mem, simblkddaskr, T0, yy, yp, jacpsol, psol);
2484 flag = IDAInit(dae_mem, simblkdaskr, T0, yy, yp);
2486 if (check_flag(&flag, "IDAInit", 1))
2488 *ierr = 200 + (-flag);
2495 N_VDestroy_Serial(IDx);
2499 N_VDestroy_Serial(yp);
2503 N_VDestroy_Serial(yy);
2520 flag = DAESStolerances(dae_mem, reltol, abstol);
2521 if (check_flag(&flag, "IDASStolerances", 1))
2523 *ierr = 200 + (-flag);
2530 N_VDestroy_Serial(IDx);
2534 N_VDestroy_Serial(yp);
2538 N_VDestroy_Serial(yy);
2555 if (solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes)
2557 flag = DDaskrRootInit(dae_mem, ng, grblkddaskr);
2561 flag = IDARootInit(dae_mem, ng, grblkdaskr);
2563 if (check_flag(&flag, "IDARootInit", 1))
2565 *ierr = 200 + (-flag);
2572 N_VDestroy_Serial(IDx);
2576 N_VDestroy_Serial(yp);
2580 N_VDestroy_Serial(yy);
2597 if (solver == IDA_BDF_Newton)
2599 flag = IDADense(dae_mem, *neq);
2601 if (check_flag(&flag, "IDADense", 1))
2603 *ierr = 200 + (-flag);
2604 if (*neq > 0)if (solver == IDA_BDF_Newton)
2610 N_VDestroy_Serial(IDx);
2614 N_VDestroy_Serial(yp);
2618 N_VDestroy_Serial(yy);
2636 if ((data = (UserData) MALLOC(sizeof(*data))) == NULL)
2645 N_VDestroy_Serial(IDx);
2649 N_VDestroy_Serial(yp);
2653 N_VDestroy_Serial(yy);
2669 data->dae_mem = dae_mem;
2675 data->ewt = N_VNew_Serial(*neq);
2676 if (check_flag((void *)data->ewt, "N_VNew_Serial", 0))
2678 *ierr = 200 + (-flag);
2689 N_VDestroy_Serial(IDx);
2693 N_VDestroy_Serial(yp);
2697 N_VDestroy_Serial(yy);
2715 if ((data->gwork = (double *) MALLOC(ng * sizeof(double))) == NULL)
2719 N_VDestroy_Serial(data->ewt);
2731 N_VDestroy_Serial(IDx);
2735 N_VDestroy_Serial(yp);
2739 N_VDestroy_Serial(yy);
2756 /*Jacobian_Flag=0; */
2757 if (AJacobian_block > 0) /* set by the block with A-Jac in flag-4 using Set_Jacobian_flag(1); */
2760 Jnx = Blocks[AJacobian_block - 1].nx;
2761 Jno = Blocks[AJacobian_block - 1].nout;
2762 Jni = Blocks[AJacobian_block - 1].nin;
2771 Jactaille = 3 * Jn + (Jn + Jni) * (Jn + Jno) + Jnx * (Jni + 2 * Jn + Jno) + (Jn - Jnx) * (2 * (Jn - Jnx) + Jno + Jni) + 2 * Jni * Jno;
2773 if ((data->rwork = (double *) MALLOC(Jactaille * sizeof(double))) == NULL)
2781 N_VDestroy_Serial(data->ewt);
2793 N_VDestroy_Serial(IDx);
2797 N_VDestroy_Serial(yp);
2801 N_VDestroy_Serial(yy);
2819 if (solver == IDA_BDF_Newton)
2821 flag = IDADlsSetDenseJacFn(dae_mem, Jacobians);
2823 if (check_flag(&flag, "IDADlsSetDenseJacFn", 1))
2825 *ierr = 200 + (-flag);
2830 TJacque = (DlsMat) NewDenseMat(*neq, *neq);
2832 flag = DAESetUserData(dae_mem, data);
2833 if (check_flag(&flag, "IDASetUserData", 1))
2835 *ierr = 200 + (-flag);
2842 flag = DAESetMaxStep(dae_mem, (realtype) hmax);
2843 if (check_flag(&flag, "IDASetMaxStep", 1))
2845 *ierr = 200 + (-flag);
2851 maxnj = 100; /* setting the maximum number of Jacobian evaluations during a Newton step */
2852 flag = DAESetMaxNumJacsIC(dae_mem, maxnj);
2853 if (check_flag(&flag, "IDASetMaxNumJacsIC", 1))
2855 *ierr = 200 + (-flag);
2860 maxnit = 10; /* setting the maximum number of Newton iterations in any attempt to solve CIC */
2861 if (C2F(cmsolver).solver == 102)
2863 maxnit = 15; /* By default, the Krylov max iterations should be 15 */
2865 flag = DAESetMaxNumItersIC(dae_mem, maxnit);
2866 if (check_flag(&flag, "IDASetMaxNumItersIC", 1))
2868 *ierr = 200 + (-flag);
2873 /* setting the maximum number of steps in an integration interval */
2875 flag = DAESetMaxNumSteps(dae_mem, maxnh);
2876 if (check_flag(&flag, "IDASetMaxNumSteps", 1))
2878 *ierr = 200 + (-flag);
2883 } /* testing if neq>0 */
2888 uround = uround * 0.5;
2890 while ( 1.0 + uround != 1.0);
2891 uround = uround * 2.0;
2892 SQuround = sqrt(uround);
2895 C2F(coshlt).halt = 0;
2902 C2F(xscion)(&inxsci);
2903 /* initialization */
2904 C2F(realtimeinit)(told, &C2F(rtfactor).scale);
2905 /* ATOL and RTOL are scalars */
2908 for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
2910 if (Blocks[C2F(curblk).kfun - 1].ng >= 1)
2912 zcros[jj] = C2F(curblk).kfun;
2916 /* . ng >= jj required */
2921 /* initialization (propagation of constant blocks outputs) */
2929 /*--discrete zero crossings----dzero--------------------*/
2930 if (ng > 0) /* storing ZC signs just after a solver call*/
2932 zdoit(told, x, x, g);
2938 for (jj = 0; jj < ng; ++jj)
2948 /* main loop on time */
2951 while (ismenu()) //** if the user has done something, do the actions
2954 SeqSync = GetCommand(CommandToUnstack); //** get to the action
2955 CommandLength = (int)strlen(CommandToUnstack);
2956 //syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
2958 if (C2F(coshlt).halt != 0)
2960 if (C2F(coshlt).halt == 2)
2962 *told = *tf; /* end simulation */
2964 C2F(coshlt).halt = 0;
2976 if (fabs(t - *told) < ttol)
2979 /* update output part */
2983 /* ! scheduling problem */
2990 if (xptr[nblk + 1] == 1)
2992 /* . no continuous state */
2993 if (*told + deltat + ttol > t)
3001 /* . update outputs of 'c' type blocks with no continuous state */
3004 /* . we are at the end, update continuous part before leaving */
3012 rhotmp = *tf + ttol;
3017 if (critev[kpo] == 1)
3019 rhotmp = tevts[kpo];
3030 hot = 0;/* Cold-restart the solver if the new TSTOP isn't beyong the previous one*/
3034 t = Min(*told + deltat, Min(t, *tf + ttol));
3036 if (hot == 0) /* CIC calculation when hot==0 */
3039 /* Setting the stop time*/
3040 flag = DAESetStopTime(dae_mem, (realtype)tstop);
3041 if (check_flag(&flag, "IDASetStopTime", 1))
3043 *ierr = 200 + (-flag);
3048 if (ng > 0 && nmod > 0)
3051 zdoit(told, x, xd, g);
3059 /*----------ID setting/checking------------*/
3060 N_VConst(ONE, IDx); /* Initialize id to 1's. */
3061 scicos_xproperty = NV_DATA_S(IDx);
3068 for (jj = 0; jj < *neq; jj++)
3072 scicos_xproperty[jj] = ONE;
3074 if (xprop[jj] == -1)
3076 scicos_xproperty[jj] = ZERO;
3079 /* CI=0.0;CJ=100.0; // for functions Get_Jacobian_ci and Get_Jacobian_cj
3080 Jacobians(*neq, (realtype) (*told), yy, yp, bidon, (realtype) CJ, data, TJacque, tempv1, tempv2, tempv3);
3081 for (jj=0;jj<*neq;jj++){
3082 Jacque_col=DENSE_COL(TJacque,jj);
3084 for (kk=0;kk<*neq;kk++){
3085 if ((Jacque_col[kk]-Jacque_col[kk]!=0)) {
3089 if (Jacque_col[kk]!=0){
3095 if (CI>=ZERO){ scicos_xproperty[jj]=CI;}else{fprintf(stderr,"\nWarinng! Xproperties are not match for i=%d!",jj);}
3097 /* printf("\n"); for(jj=0;jj<*neq;jj++) { printf("x%d=%g ",jj,scicos_xproperty[jj]); }*/
3099 flag = DAESetId(dae_mem, IDx);
3100 if (check_flag(&flag, "IDASetId", 1))
3102 *ierr = 200 + (-flag);
3106 // CI=1.0; // for function Get_Jacobian_ci
3107 /*--------------------------------------------*/
3108 // maxnj=100; /* setting the maximum number of Jacobian evaluation during a Newton step */
3109 // flag=DAESetMaxNumJacsIC(dae_mem, maxnj);
3110 // if (check_flag(&flag, "IDASetMaxNumJacsIC", 1)) {
3111 // *ierr=200+(-flag);
3115 // flag=DAESetLineSearchOffIC(dae_mem,FALSE); /* (def=false) */
3116 // if (check_flag(&flag, "IDASetLineSearchOffIC", 1)) {
3117 // *ierr=200+(-flag);
3121 // flag=DAESetMaxNumItersIC(dae_mem, 10);/* (def=10) setting the maximum number of Newton iterations in any attempt to solve CIC */
3122 // if (check_flag(&flag, "IDASetMaxNumItersIC", 1)) {
3123 // *ierr=200+(-flag);
3128 N_iters = 4 + nmod * 4;
3129 for (j = 0; j <= N_iters; j++)
3131 /* counter to reevaluate the
3132 modes in mode->CIC->mode->CIC-> loop
3133 do it once in the absence of mode (nmod=0) */
3134 /* updating the modes through Flag==9, Phase==1 */
3136 /* Serge Steer 29/06/2009 */
3137 while (ismenu()) //** if the user has done something, do the actions
3140 SeqSync = GetCommand(CommandToUnstack); //** get to the action
3141 CommandLength = (int)strlen(CommandToUnstack);
3142 //syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
3144 if (C2F(coshlt).halt != 0)
3146 if (C2F(coshlt).halt == 2)
3148 *told = *tf; /* end simulation */
3150 C2F(coshlt).halt = 0;
3156 flag = DAEReInit(dae_mem, (realtype)(*told), yy, yp);
3157 if (check_flag(&flag, "IDAReInit", 1))
3159 *ierr = 200 + (-flag);
3164 phase = 2; /* IDACalcIC: PHI-> yy0: if (ok) yy0_cic-> PHI*/
3165 if (C2F(cmsolver).solver == 100)
3167 copy_IDA_mem->ida_kk = 1;
3169 // the initial conditons y0 and yp0 do not satisfy the DAE
3170 flagr = DAECalcIC(dae_mem, DAE_YA_YDP_INIT, (realtype)(t));
3172 flag = DAEGetConsistentIC(dae_mem, yy, yp); /* PHI->YY */
3173 if (*ierr > 5) /* *ierr>5 => singularity in block */
3179 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3183 sciprint(_("**** SUNDIALS.IDA successfully initialized *****\n") );
3187 sciprint(_("**** SUNDIALS.IDA failed to initialize ->try again *****\n") );
3190 /*-------------------------------------*/
3191 /* saving the previous modes*/
3192 for (jj = 0; jj < nmod; ++jj)
3194 Mode_save[jj] = mod[jj];
3196 if (ng > 0 && nmod > 0)
3199 zdoit(told, x, xd, g);
3206 /*------------------------------------*/
3208 for (jj = 0; jj < nmod; ++jj)
3210 if (Mode_save[jj] != mod[jj])
3216 if (Mode_change == 0)
3220 break; /* if (flagr>=0) break; else{ *ierr=200+(-flagr); freeallx; return; }*/
3222 else if (j >= (int)( N_iters / 2))
3224 /* DAESetMaxNumStepsIC(mem,10); */ /* maxnh (def=5) */
3225 DAESetMaxNumJacsIC(dae_mem, 10); /* maxnj 100 (def=4)*/
3226 /* DAESetMaxNumItersIC(mem,100000); */ /* maxnit in IDANewtonIC (def=10) */
3227 DAESetLineSearchOffIC(dae_mem, TRUE); /* (def=false) */
3228 /* DAESetNonlinConvCoefIC(mem,1.01);*/ /* (def=0.01-0.33*/
3229 flag = DAESetMaxNumItersIC(dae_mem, 1000);
3230 if (check_flag(&flag, "IDASetMaxNumItersIC", 1))
3232 *ierr = 200 + (-flag);
3238 }/* mode-CIC counter*/
3239 if (Mode_change == 1)
3241 /* In this case, we try again by relaxing all modes and calling DAECalcIC again
3244 if (C2F(cmsolver).solver == 100)
3246 copy_IDA_mem->ida_kk = 1;
3248 flagr = DAECalcIC(dae_mem, DAE_YA_YDP_INIT, (realtype)(t));
3250 flag = DAEGetConsistentIC(dae_mem, yy, yp); /* PHI->YY */
3251 if ((flagr < 0) || (*ierr > 5)) /* *ierr>5 => singularity in block */
3258 /*-----If flagr<0 the initialization solver has not converged-----*/
3266 } /* CIC calculation when hot==0 */
3268 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3270 sciprint(_("****daskr from: %f to %f hot= %d \n"), *told, t, hot);
3273 /*--discrete zero crossings----dzero--------------------*/
3274 /*--check for Dzeros after Mode settings or ddoit()----*/
3276 if (ng > 0 && hot == 0)
3278 zdoit(told, x, xd, g);
3284 for (jj = 0; jj < ng; ++jj)
3286 if ((g[jj] >= 0.0) && ( jroot[jj] == -5))
3291 else if ((g[jj] < 0.0) && ( jroot[jj] == 5))
3303 /*--discrete zero crossings----dzero--------------------*/
3304 if (Discrete_Jump == 0) /* if there was a dzero, its event should be activated*/
3307 flagr = DAESolve(dae_mem, t, told, yy, yp, DAE_NORMAL);
3317 flagr = IDA_ROOT_RETURN; /* in order to handle discrete jumps */
3321 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3323 sciprint(_("****SUNDIALS.Ida reached: %f\n"), *told);
3328 else if ( flagr == IDA_TOO_MUCH_WORK || flagr == IDA_CONV_FAIL || flagr == IDA_ERR_FAIL)
3330 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3332 sciprint(_("**** SUNDIALS.Ida: too much work at time=%g (stiff region, change RTOL and ATOL)\n"), *told);
3338 *ierr = 200 + (-flagr);
3347 *ierr = 200 + (-flagr); /* raising errors due to internal errors, otherwise error due to flagr*/
3353 /* update outputs of 'c' type blocks if we are at the end*/
3361 if (flagr == IDA_ZERO_DETACH_RETURN)
3364 }; /* new feature of sundials, detects unmasking */
3365 if (flagr == IDA_ROOT_RETURN)
3367 /* . at least one root has been found */
3369 if (Discrete_Jump == 0)
3371 flagr = DAEGetRootInfo(dae_mem, jroot);
3372 if (check_flag(&flagr, "IDAGetRootInfo", 1))
3374 *ierr = 200 + (-flagr);
3380 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3382 sciprint(_("root found at t=: %f\n"), *told);
3384 /* . update outputs affecting ztyp blocks ONLY FOR OLD BLOCKS*/
3385 zdoit(told, x, xd, g);
3391 for (jj = 0; jj < ng; ++jj)
3393 C2F(curblk).kfun = zcros[jj];
3394 if (C2F(curblk).kfun == -1)
3399 for (j = zcptr[C2F(curblk).kfun] - 1 ;
3400 j < zcptr[C2F(curblk).kfun + 1] - 1 ; ++j)
3410 Blocks[C2F(curblk).kfun - 1].jroot = &jroot[zcptr[C2F(curblk).kfun] - 1];
3411 if (funtyp[C2F(curblk).kfun] > 0)
3413 if (Blocks[C2F(curblk).kfun - 1].nevout > 0)
3416 if (Blocks[C2F(curblk).kfun - 1].nx > 0)
3418 Blocks[C2F(curblk).kfun - 1].x = &x[xptr[C2F(curblk).kfun] - 1];
3419 Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
3421 /* call corresponding block to determine output event (kev) */
3422 Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
3423 /*callf(told, xd, x, x,g,&flag__);*/
3424 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
3431 /* update event agenda */
3432 for (k = 0; k < Blocks[C2F(curblk).kfun - 1].nevout; ++k)
3434 if (Blocks[C2F(curblk).kfun - 1].evout[k] >= 0)
3436 i3 = k + clkptr[C2F(curblk).kfun] ;
3437 addevs(Blocks[C2F(curblk).kfun - 1].evout[k] + (*told), &i3, &ierr1);
3440 /* . nevts too small */
3449 if ((Blocks[C2F(curblk).kfun - 1].nx > 0) || (*Blocks[C2F(curblk).kfun - 1].work != NULL) )
3451 /* call corresponding block to update state */
3453 if (Blocks[C2F(curblk).kfun - 1].nx > 0)
3455 Blocks[C2F(curblk).kfun - 1].x = &x[xptr[C2F(curblk).kfun] - 1];
3456 Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
3458 Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
3460 Blocks[C2F(curblk).kfun - 1].xprop = &xprop[-1 + xptr[C2F(curblk).kfun]];
3461 /*callf(told, xd, x, x,g,&flag__);*/
3462 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
3470 for (j = 0; j < *neq; j++) /* Adjust xprop for IDx */
3474 scicos_xproperty[j] = ONE;
3478 scicos_xproperty[j] = ZERO;
3486 /* Serge Steer 29/06/2009 */
3487 while (ismenu()) //** if the user has done something, do the actions
3490 SeqSync = GetCommand(CommandToUnstack); //** get to the action
3491 CommandLength = (int)strlen(CommandToUnstack);
3492 //syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
3495 if (C2F(coshlt).halt != 0)
3497 if (C2F(coshlt).halt == 2)
3499 *told = *tf; /* end simulation */
3501 C2F(coshlt).halt = 0;
3518 /*--discrete zero crossings----dzero--------------------*/
3519 if (ng > 0) /* storing ZC signs just after a ddaskr call*/
3521 zdoit(told, x, xd, g);
3527 for (jj = 0; jj < ng; ++jj)
3539 /*--discrete zero crossings----dzero--------------------*/
3541 C2F(realtime)(told);
3546 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3548 sciprint(_("Event: %d activated at t=%f\n"), *pointi, *told);
3552 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3554 sciprint(_("End of activation\n"));
3562 /* end of main loop on time */
3565 } /* cossimdaskr_ */
3566 /*--------------------------------------------------------------------------*/
3567 /* Subroutine cosend */
3568 static void cosend(double *told)
3570 /* Local variables */
3571 static scicos_flag flag__ = 0;
3573 static int kfune = 0;
3577 /* loop on blocks */
3578 for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
3581 Blocks[C2F(curblk).kfun - 1].nevprt = 0;
3582 if (funtyp[C2F(curblk).kfun] >= 0)
3584 if (Blocks[C2F(curblk).kfun - 1].nx > 0)
3586 Blocks[C2F(curblk).kfun - 1].x = &x[xptr[C2F(curblk).kfun] - 1];
3587 Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
3589 /*callf(told, xd, x, x,x,&flag__);*/
3590 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
3591 if (flag__ < 0 && *ierr == 0)
3594 kfune = C2F(curblk).kfun;
3600 C2F(curblk).kfun = kfune;
3604 /*--------------------------------------------------------------------------*/
3606 void callf(double *t, scicos_block *block, scicos_flag *flag)
3608 double* args[SZ_SIZE];
3610 double intabl[TB_SIZE];
3611 double outabl[TB_SIZE];
3613 int ii = 0, in = 0, out = 0, ki = 0, ko = 0, no = 0, ni = 0, k = 0, j = 0;
3614 int szi = 0, flagi = 0;
3615 double *ptr_d = NULL;
3617 /* function pointers type def */
3621 /* ScicosFm1 loc3;*/
3629 int solver = C2F(cmsolver).solver;
3630 int cosd = C2F(cosdebug).cosd;
3631 /*int kf = C2F(curblk).kfun;*/
3633 block_error = (int*) flag;
3635 /* debug block is never called */
3636 /*if (kf==(debug_block+1)) return;*/
3637 if (block->type == 99)
3642 /* flag 7 implicit initialization */
3643 flagi = (int) * flag;
3644 /* change flag to zero if flagi==7 for explicit block */
3645 if (flagi == 7 && block->type < 10000)
3650 /* display information for debugging mode */
3655 sciprint(_("block %d is called "), C2F(curblk).kfun);
3656 sciprint(_("with flag %d "), *flag);
3657 sciprint(_("at time %f \n"), *t);
3659 if (debug_block > -1)
3663 sciprint(_("Entering the block \n"));
3665 call_debug_scicos(block, flag, flagi, debug_block);
3668 return; /* error in debug block */
3673 C2F(scsptr).ptr = block->scsptr;
3675 /* get pointer of the function */
3678 /* continuous state */
3679 if ((solver == IDA_BDF_Newton || solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes) && block->type < 10000 && *flag == 0)
3682 block->xd = block->res;
3686 //sciprint("callf type=%d flag=%d\n",block->type,flagi);
3687 switch (block->type)
3689 /*******************/
3690 /* function type 0 */
3691 /*******************/
3694 /* This is for compatibility */
3695 /* jroot is returned in g for old type */
3696 if (block->nevprt < 0)
3698 for (j = 0; j < block->ng; ++j)
3700 block->g[j] = (double)block->jroot[j];
3704 /* concatenated entries and concatened outputs */
3705 /* catenate inputs if necessary */
3710 for (in = 0; in < block->nin; in++)
3712 szi = block->insz[in] * block->insz[in + block->nin];
3713 for (ii = 0; ii < szi; ii++)
3715 intabl[ki++] = *((double *)(block->inptr[in]) + ii);
3719 args[0] = &(intabl[0]);
3723 if (block->nin == 0)
3729 args[0] = (double *)(block->inptr[0]);
3730 ni = block->insz[0] * block->insz[1];
3734 /* catenate outputs if necessary */
3736 if (block->nout > 1)
3739 for (out = 0; out < block->nout; out++)
3741 szi = block->outsz[out] * block->outsz[out + block->nout];
3742 for (ii = 0; ii < szi; ii++)
3744 outabl[ko++] = *((double *)(block->outptr[out]) + ii);
3748 args[1] = &(outabl[0]);
3752 if (block->nout == 0)
3758 args[1] = (double *)(block->outptr[0]);
3759 no = block->outsz[0] * block->outsz[1];
3763 loc0 = (ScicosF0) loc;
3765 (*loc0)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3766 block->z, &block->nz,
3767 block->evout, &block->nevout, block->rpar, &block->nrpar,
3768 block->ipar, &block->nipar, (double *)args[0], &ni,
3769 (double *)args[1], &no);
3771 /* split output vector on each port if necessary */
3772 if (block->nout > 1)
3775 for (out = 0; out < block->nout; out++)
3777 szi = block->outsz[out] * block->outsz[out + block->nout];
3778 for (ii = 0; ii < szi; ii++)
3780 *((double *)(block->outptr[out]) + ii) = outabl[ko++];
3785 /* adjust values of output register */
3786 for (in = 0; in < block->nevout; ++in)
3788 block->evout[in] = block->evout[in] - *t;
3794 /*******************/
3795 /* function type 1 */
3796 /*******************/
3799 /* This is for compatibility */
3800 /* jroot is returned in g for old type */
3801 if (block->nevprt < 0)
3803 for (j = 0; j < block->ng; ++j)
3805 block->g[j] = (double)block->jroot[j];
3809 /* one entry for each input or output */
3810 for (in = 0 ; in < block->nin ; in++)
3812 args[in] = block->inptr[in];
3813 sz[in] = block->insz[in];
3815 for (out = 0; out < block->nout; out++)
3817 args[in + out] = block->outptr[out];
3818 sz[in + out] = block->outsz[out];
3820 /* with zero crossing */
3821 if (block->ztyp > 0)
3823 args[block->nin + block->nout] = block->g;
3824 sz[block->nin + block->nout] = block->ng;
3827 loc1 = (ScicosF) loc;
3829 (*loc1)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3830 block->z, &block->nz,
3831 block->evout, &block->nevout, block->rpar, &block->nrpar,
3832 block->ipar, &block->nipar,
3833 (double *)args[0], &sz[0],
3834 (double *)args[1], &sz[1], (double *)args[2], &sz[2],
3835 (double *)args[3], &sz[3], (double *)args[4], &sz[4],
3836 (double *)args[5], &sz[5], (double *)args[6], &sz[6],
3837 (double *)args[7], &sz[7], (double *)args[8], &sz[8],
3838 (double *)args[9], &sz[9], (double *)args[10], &sz[10],
3839 (double *)args[11], &sz[11], (double *)args[12], &sz[12],
3840 (double *)args[13], &sz[13], (double *)args[14], &sz[14],
3841 (double *)args[15], &sz[15], (double *)args[16], &sz[16],
3842 (double *)args[17], &sz[17]);
3844 /* adjust values of output register */
3845 for (in = 0; in < block->nevout; ++in)
3847 block->evout[in] = block->evout[in] - *t;
3853 /*******************/
3854 /* function type 2 */
3855 /*******************/
3858 /* This is for compatibility */
3859 /* jroot is returned in g for old type */
3860 if (block->nevprt < 0)
3862 for (j = 0; j < block->ng; ++j)
3864 block->g[j] = (double)block->jroot[j];
3868 /* no zero crossing */
3869 if (block->ztyp == 0)
3871 loc2 = (ScicosF2) loc;
3872 (*loc2)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3873 block->z, &block->nz,
3874 block->evout, &block->nevout, block->rpar, &block->nrpar,
3875 block->ipar, &block->nipar, (double **)block->inptr,
3876 block->insz, &block->nin,
3877 (double **)block->outptr, block->outsz, &block->nout);
3879 /* with zero crossing */
3882 loc2z = (ScicosF2z) loc;
3883 (*loc2z)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3884 block->z, &block->nz,
3885 block->evout, &block->nevout, block->rpar, &block->nrpar,
3886 block->ipar, &block->nipar, (double **)block->inptr,
3887 block->insz, &block->nin,
3888 (double **)block->outptr, block->outsz, &block->nout,
3889 block->g, &block->ng);
3892 /* adjust values of output register */
3893 for (in = 0; in < block->nevout; ++in)
3895 block->evout[in] = block->evout[in] - *t;
3901 /*******************/
3902 /* function type 4 */
3903 /*******************/
3906 /* get pointer of the function type 4*/
3907 loc4 = (ScicosF4) loc;
3909 (*loc4)(block, *flag);
3914 /***********************/
3915 /* function type 10001 */
3916 /***********************/
3919 /* This is for compatibility */
3920 /* jroot is returned in g for old type */
3921 if (block->nevprt < 0)
3923 for (j = 0; j < block->ng; ++j)
3925 block->g[j] = (double)block->jroot[j];
3929 /* implicit block one entry for each input or output */
3930 for (in = 0 ; in < block->nin ; in++)
3932 args[in] = block->inptr[in];
3933 sz[in] = block->insz[in];
3935 for (out = 0; out < block->nout; out++)
3937 args[in + out] = block->outptr[out];
3938 sz[in + out] = block->outsz[out];
3940 /* with zero crossing */
3941 if (block->ztyp > 0)
3943 args[block->nin + block->nout] = block->g;
3944 sz[block->nin + block->nout] = block->ng;
3947 loci1 = (ScicosFi) loc;
3948 (*loci1)(flag, &block->nevprt, t, block->res, block->xd, block->x,
3949 &block->nx, block->z, &block->nz,
3950 block->evout, &block->nevout, block->rpar, &block->nrpar,
3951 block->ipar, &block->nipar,
3952 (double *)args[0], &sz[0],
3953 (double *)args[1], &sz[1], (double *)args[2], &sz[2],
3954 (double *)args[3], &sz[3], (double *)args[4], &sz[4],
3955 (double *)args[5], &sz[5], (double *)args[6], &sz[6],
3956 (double *)args[7], &sz[7], (double *)args[8], &sz[8],
3957 (double *)args[9], &sz[9], (double *)args[10], &sz[10],
3958 (double *)args[11], &sz[11], (double *)args[12], &sz[12],
3959 (double *)args[13], &sz[13], (double *)args[14], &sz[14],
3960 (double *)args[15], &sz[15], (double *)args[16], &sz[16],
3961 (double *)args[17], &sz[17]);
3963 /* adjust values of output register */
3964 for (in = 0; in < block->nevout; ++in)
3966 block->evout[in] = block->evout[in] - *t;
3972 /***********************/
3973 /* function type 10002 */
3974 /***********************/
3977 /* This is for compatibility */
3978 /* jroot is returned in g for old type */
3979 if (block->nevprt < 0)
3981 for (j = 0; j < block->ng; ++j)
3983 block->g[j] = (double)block->jroot[j];
3987 /* implicit block, inputs and outputs given by a table of pointers */
3988 /* no zero crossing */
3989 if (block->ztyp == 0)
3991 loci2 = (ScicosFi2) loc;
3992 (*loci2)(flag, &block->nevprt, t, block->res,
3993 block->xd, block->x, &block->nx,
3994 block->z, &block->nz,
3995 block->evout, &block->nevout, block->rpar, &block->nrpar,
3996 block->ipar, &block->nipar, (double **)block->inptr,
3997 block->insz, &block->nin,
3998 (double **)block->outptr, block->outsz, &block->nout);
4000 /* with zero crossing */
4003 loci2z = (ScicosFi2z) loc;
4004 (*loci2z)(flag, &block->nevprt, t, block->res,
4005 block->xd, block->x, &block->nx,
4006 block->z, &block->nz,
4007 block->evout, &block->nevout, block->rpar, &block->nrpar,
4008 block->ipar, &block->nipar,
4009 (double **)block->inptr, block->insz, &block->nin,
4010 (double **)block->outptr, block->outsz, &block->nout,
4011 block->g, &block->ng);
4014 /* adjust values of output register */
4015 for (in = 0; in < block->nevout; ++in)
4017 block->evout[in] = block->evout[in] - *t;
4023 /***********************/
4024 /* function type 10004 */
4025 /***********************/
4028 /* get pointer of the function type 4*/
4029 loc4 = (ScicosF4) loc;
4031 (*loc4)(block, *flag);
4041 sciprint(_("Undefined Function type\n"));
4046 // sciprint("callf end flag=%d\n",*flag);
4047 /* Implicit Solver & explicit block & flag==0 */
4048 /* adjust continuous state vector after call */
4049 if ((solver == IDA_BDF_Newton || solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes) && block->type < 10000 && *flag == 0)
4054 for (k = 0; k < block->nx; k++)
4056 block->res[k] = block->res[k] - block->xd[k];
4061 for (k = 0; k < block->nx; k++)
4063 block->xd[k] = block->res[k];
4071 if (debug_block > -1)
4075 return; /* error in block */
4079 sciprint(_("Leaving block %d \n"), C2F(curblk).kfun);
4081 call_debug_scicos(block, flag, flagi, debug_block);
4082 /*call_debug_scicos(flag,kf,flagi,debug_block);*/
4086 /*--------------------------------------------------------------------------*/
4087 /* call_debug_scicos */
4088 static void call_debug_scicos(scicos_block *block, scicos_flag *flag, int flagi, int deb_blk)
4091 int solver = C2F(cmsolver).solver, k = 0;
4093 double *ptr_d = NULL;
4095 C2F(cosdebugcounter).counter = C2F(cosdebugcounter).counter + 1;
4096 C2F(scsptr).ptr = Blocks[deb_blk].scsptr;
4098 loc = Blocks[deb_blk].funpt; /* GLOBAL */
4099 loc4 = (ScicosF4) loc;
4101 /* continuous state */
4102 if ((solver == IDA_BDF_Newton || solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes) && block->type < 10000 && *flag == 0)
4105 block->xd = block->res;
4108 (*loc4)(block, *flag);
4110 /* Implicit Solver & explicit block & flag==0 */
4111 /* adjust continuous state vector after call */
4112 if ((solver == IDA_BDF_Newton || solver == DDaskr_BDF_Newton || solver == DDaskr_BDF_GMRes) && block->type < 10000 && *flag == 0)
4117 for (k = 0; k < block->nx; k++)
4119 block->res[k] = block->res[k] - block->xd[k];
4124 for (k = 0; k < block->nx; k++)
4126 block->xd[k] = block->res[k];
4133 sciprint(_("Error in the Debug block \n"));
4135 } /* call_debug_scicos */
4136 /*--------------------------------------------------------------------------*/
4138 static int simblk(realtype t, N_Vector yy, N_Vector yp, void *f_data)
4140 double tx = 0., *x = NULL, *xd = NULL;
4141 int i = 0, nantest = 0;
4144 x = (double *) NV_DATA_S(yy);
4145 xd = (double *) NV_DATA_S(yp);
4147 for (i = 0; i < *neq; i++)
4149 xd[i] = 0; /* Ã la place de "C2F(dset)(neq, &c_b14,xcdot , &c__1);"*/
4151 C2F(ierode).iero = 0;
4153 odoit(&tx, x, xd, xd);
4154 C2F(ierode).iero = *ierr;
4159 for (i = 0; i < *neq; i++) /* NaN checking */
4161 if ((xd[i] - xd[i] != 0))
4163 sciprint(_("\nWarning: The computing function #%d returns a NaN/Inf"), i);
4170 return 349; /* recoverable error; */
4174 return (abs(*ierr)); /* ierr>0 recoverable error; ierr>0 unrecoverable error; ierr=0: ok*/
4177 /*--------------------------------------------------------------------------*/
4179 static int grblk(realtype t, N_Vector yy, realtype *gout, void *g_data)
4181 double tx = 0., *x = NULL;
4182 int jj = 0, nantest = 0;
4185 x = (double *) NV_DATA_S(yy);
4187 C2F(ierode).iero = 0;
4190 zdoit(&tx, x, x, (double*) gout);
4195 for (jj = 0; jj < ng; jj++)
4196 if (gout[jj] - gout[jj] != 0)
4198 sciprint(_("\nWarning: The zero_crossing function #%d returns a NaN/Inf"), jj);
4201 } /* NaN checking */
4204 return 350; /* recoverable error; */
4207 C2F(ierode).iero = *ierr;
4211 /*--------------------------------------------------------------------------*/
4213 static void simblklsodar(int * nequations, realtype * tOld, realtype * actual, realtype * res)
4218 tx = (double) * tOld;
4220 for (i = 0; i < *nequations; ++i)
4222 res[i] = 0; /* Ã la place de "C2F(dset)(neq, &c_b14,xcdot , &c__1);"*/
4224 C2F(ierode).iero = 0;
4226 odoit(&tx, actual, res, res);
4227 C2F(ierode).iero = *ierr;
4231 for (i = 0; i < *nequations; i++) /* NaN checking */
4233 if ((res[i] - res[i] != 0))
4235 sciprint(_("\nWarning: The computing function #%d returns a NaN/Inf"), i);
4239 } /* simblklsodar */
4240 /*--------------------------------------------------------------------------*/
4242 static void grblklsodar(int * nequations, realtype * tOld, realtype * actual, int * ngc, realtype * res)
4247 tx = (double) * tOld;
4249 C2F(ierode).iero = 0;
4252 zdoit(&tx, actual, actual, res);
4256 for (jj = 0; jj < *ngc; jj++)
4258 if (res[jj] - res[jj] != 0)
4260 sciprint(_("\nWarning: The zero_crossing function #%d returns a NaN/Inf"), jj);
4261 } /* NaN checking */
4265 /*--------------------------------------------------------------------------*/
4267 static int simblkdaskr(realtype tres, N_Vector yy, N_Vector yp, N_Vector resval, void *rdata)
4270 double *xc = NULL, *xcdot = NULL, *residual = NULL;
4271 realtype alpha = 0.;
4277 int jj = 0, flag = 0, nantest = 0;
4279 data = (UserData) rdata;
4281 if (get_phase_simulation() == 1)
4283 /* Just to update mode in a very special case, i.e., when initialization using modes fails.
4284 in this case, we relax all modes and try again one more time.
4286 zdoit(&tx, NV_DATA_S(yy), NV_DATA_S(yp), NULL);
4290 flag = IDAGetCurrentStep(data->dae_mem, &hh);
4293 *ierr = 200 + (-flag);
4298 flag = IDAGetCurrentOrder(data->dae_mem, &qlast);
4301 *ierr = 200 + (-flag);
4306 for (jj = 0; jj < qlast; jj++)
4308 alpha = alpha - ONE / (jj + 1);
4311 // CJ=-alpha/hh; // For function Get_Jacobian_cj
4320 xc = (double *) NV_DATA_S(yy);
4321 xcdot = (double *) NV_DATA_S(yp);
4322 residual = (double *) NV_DATA_S(resval);
4325 C2F(dcopy)(neq, xcdot, &c__1, residual, &c__1);
4327 C2F(ierode).iero = 0;
4328 odoit(&tx, xc, xcdot, residual);
4330 C2F(ierode).iero = *ierr;
4335 for (jj = 0; jj < *neq; jj++)
4336 if (residual[jj] - residual[jj] != 0) /* NaN checking */
4338 //sciprint(_("\nWarning: The residual function #%d returns a NaN"),jj);
4344 return 257; /* recoverable error; */
4348 return (abs(*ierr)); /* ierr>0 recoverable error; ierr>0 unrecoverable error; ierr=0: ok*/
4350 /*--------------------------------------------------------------------------*/
4352 static int grblkdaskr(realtype t, N_Vector yy, N_Vector yp, realtype *gout, void *g_data)
4355 int jj = 0, nantest = 0;
4360 C2F(ierode).iero = 0;
4361 zdoit(&tx, NV_DATA_S(yy), NV_DATA_S(yp), (double *) gout);
4364 nantest = 0; /* NaN checking */
4365 for (jj = 0; jj < ng; jj++)
4367 if (gout[jj] - gout[jj] != 0)
4369 sciprint(_("\nWarning: The zero-crossing function #%d returns a NaN"), jj);
4376 return 258; /* recoverable error; */
4379 C2F(ierode).iero = *ierr;
4382 /*--------------------------------------------------------------------------*/
4384 static void simblkddaskr(realtype *tOld, realtype *actual, realtype *actualP, realtype *res, int *flag, double *dummy1, int *dummy2)
4390 if (get_phase_simulation() == 1)
4392 /* Just to update mode in a very special case, i.e., when initialization using modes fails.
4393 in this case, we relax all modes and try again one more time.
4395 zdoit(&tx, actual, actualP, NULL);
4400 tx = (double) * tOld;
4403 C2F(dcopy)(neq, actualP, &c__1, res, &c__1);
4405 C2F(ierode).iero = 0;
4406 odoit(&tx, actual, actualP, res);
4407 C2F(ierode).iero = *ierr;
4411 for (jj = 0; jj < *neq; jj++)
4412 if (res[jj] - res[jj] != 0) /* NaN checking */
4414 sciprint(_("\nWarning: The residual function #%d returns a NaN"), jj);
4415 *flag = -1; /* recoverable error; */
4424 /* *flag=-1 recoverable error; *flag=-2 unrecoverable error; *flag=0: ok*/
4427 /*--------------------------------------------------------------------------*/
4429 static void grblkddaskr(int *nequations, realtype *tOld, realtype *actual, int *ngc, realtype *res, double *dummy1, int *dummy2)
4434 tx = (double) * tOld;
4437 C2F(ierode).iero = 0;
4438 zdoit(&tx, actual, actual, res);
4439 C2F(ierode).iero = *ierr;
4444 for (jj = 0; jj < *ngc; jj++)
4446 if (res[jj] - res[jj] != 0)
4448 sciprint(_("\nWarning: The zero-crossing function #%d returns a NaN"), jj);
4455 sciprint(_("\nError: Problem in the evaluation of a root function"));
4460 /*--------------------------------------------------------------------------*/
4462 static void jacpsol(realtype *res, int *ires, int *neq, realtype *tOld, realtype *actual, realtype *actualP,
4463 realtype *rewt, realtype *savr, realtype *wk, realtype *h, realtype *cj, realtype *wp, int *iwp,
4464 int *ier, double *dummy1, int *dummy2)
4466 /* Here, we compute the system preconditioner matrix P, which is actually the jacobian matrix,
4467 so P(i,j) = dres(i)/dactual(j) + cj*dres(i)/dactualP(j), and we LU-decompose it. */
4468 int i = 0, j = 0, nrow = 0, info = 0;
4469 realtype tx = 0, del = 0, delinv = 0, ysave = 0, ypsave = 0;
4470 realtype * e = NULL;
4474 /* Work array used to evaluate res(*tOld, actual + small_increment, actualP + small_increment).
4475 savr already contains res(*tOld, actual, actualP). */
4476 e = (realtype *) calloc(*neq, sizeof(realtype));
4478 for (i = 0; i < *neq; ++i)
4480 del = max (SQuround * max(fabs(actual[i]), fabs(*h * actualP[i])), 1. / rewt[i]);
4481 del *= (*h * actualP[i] >= 0) ? 1 : -1;
4482 del = (actual[i] + del) - actual[i];
4484 ypsave = actualP[i];
4486 actualP[i] += *cj * del;
4488 C2F(ierode).iero = 0;
4489 simblkddaskr(tOld, actual, actualP, e, ires, dummy1, dummy2);
4490 C2F(ierode).iero = *ierr;
4494 for (j = 0; j < *neq; ++j)
4496 wp[nrow + j] = (e[j] - savr[j]) * delinv;
4498 if (wp[nrow + j] - wp[nrow + j] != 0)
4500 sciprint(_("\nWarning: The preconditioner evaluation function returns a NaN at index #%d."), nrow + j);
4506 actualP[i] = ypsave;
4510 sciprint(_("\nError: The preconditioner evaluation function failed."));
4517 /* Proceed to LU factorization of P. */
4518 C2F(dgefa) (wp, neq, neq, iwp, &info);
4521 sciprint(_("\nError: Failed to factor the preconditioner."));
4527 /*--------------------------------------------------------------------------*/
4529 static void psol(int *neq, realtype *tOld, realtype *actual, realtype *actualP,
4530 realtype *savr, realtype *wk, realtype *cj, realtype *wght, realtype *wp,
4531 int *iwp, realtype *b, realtype *eplin, int *ier, double *dummy1, int *dummy2)
4533 /* This function "applies" the inverse of the preconditioner to 'b' (computes P^-1*b).
4534 It is done by solving P*x = b using the linpack routine 'dgesl'. */
4537 C2F(dgesl) (wp, neq, neq, iwp, b, &job);
4540 for (i = 0; i < *neq; ++i)
4542 if (b[i] - b[i] != 0)
4544 sciprint(_("\nWarning: The preconditioner application function returns a NaN at index #%d."), i);
4545 /* Indicate a recoverable error, meaning that the step will be retried with the same step size
4546 but with a call to 'jacpsol' to update necessary data, unless the Jacobian data is current,
4547 in which case the step will be retried with a smaller step size. */
4552 /*--------------------------------------------------------------------------*/
4553 /* Subroutine addevs */
4554 static void addevs(double t, int *evtnb, int *ierr1)
4556 static int i = 0, j = 0;
4560 if (evtspt[*evtnb] != -1)
4562 if ((evtspt[*evtnb] == 0) && (*pointi == *evtnb))
4569 if (*pointi == *evtnb)
4571 *pointi = evtspt[*evtnb]; /* remove from chain */
4576 while (*evtnb != evtspt[i])
4580 evtspt[i] = evtspt[*evtnb]; /* remove old evtnb from chain */
4581 if (TCritWarning == 0)
4583 sciprint(_("\n Warning: an event is reprogrammed at t=%g by removing another"), t );
4584 sciprint(_("\n (already programmed) event. There may be an error in"));
4585 sciprint(_("\n your model. Please check your model\n"));
4588 do_cold_restart(); /* the erased event could be a critical
4589 event, so do_cold_restart is added to
4590 refresh the critical event table */
4606 if (t < tevts[*pointi])
4608 evtspt[*evtnb] = *pointi;
4620 if (t >= tevts[evtspt[i]])
4633 evtspt[*evtnb] = evtspt[i];
4637 /*--------------------------------------------------------------------------*/
4638 /* Subroutine putevs */
4639 void putevs(double *t, int *evtnb, int *ierr1)
4643 if (evtspt[*evtnb] != -1)
4658 evtspt[*evtnb] = *pointi;
4661 /*--------------------------------------------------------------------------*/
4662 /* Subroutine idoit */
4663 static void idoit(double *told)
4665 /* initialisation (propagation of constant blocks outputs) */
4666 /* Copyright INRIA */
4669 scicos_flag flag = 0;
4673 ScicosImport *scs_imp = NULL;
4676 scs_imp = getscicosimportptr();
4679 for (j = 0; j < * (scs_imp->niord); j++)
4681 kf = &scs_imp->iord[j];
4682 C2F(curblk).kfun = *kf; /* */
4683 if (scs_imp->funtyp[*kf - 1] > -1)
4685 /* continuous state */
4686 if (scs_imp->blocks[*kf - 1].nx > 0)
4688 scs_imp->blocks[*kf - 1].x = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4689 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4691 scs_imp->blocks[*kf - 1].nevprt = scs_imp->iord[j + * (scs_imp->niord)];
4692 /*callf(told, xd, x, x,x,&flag);*/
4693 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4700 if (scs_imp->blocks[*kf - 1].nevout > 0)
4702 if (scs_imp->funtyp[*kf - 1] < 0)
4704 i = synchro_nev(scs_imp, *kf, ierr);
4709 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4710 putevs(told, &i2, &ierr1);
4713 /* event conflict */
4726 /*--------------------------------------------------------------------------*/
4727 /* Subroutine doit */
4728 static void doit(double *told)
4730 /* propagation of blocks outputs on discrete activations */
4731 /* Copyright INRIA */
4734 scicos_flag flag = 0;
4737 int ii = 0, kever = 0;
4739 ScicosImport *scs_imp = NULL;
4742 scs_imp = getscicosimportptr();
4745 *pointi = evtspt[kever];
4748 nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever - 1];
4754 for (ii = scs_imp->ordptr[kever - 1]; ii <= scs_imp->ordptr[kever] - 1 ; ii++)
4756 kf = &scs_imp->ordclk[ii - 1];
4757 C2F(curblk).kfun = *kf;
4758 if (scs_imp->funtyp[*kf - 1] > -1)
4760 /* continuous state */
4761 if (scs_imp->blocks[*kf - 1].nx > 0)
4763 scs_imp->blocks[*kf - 1].x = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4764 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4766 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
4768 /*callf(told, xd, x, x,x,&flag);*/
4769 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4777 /* Initialize tvec */
4778 if (scs_imp->blocks[*kf - 1].nevout > 0)
4780 if (scs_imp->funtyp[*kf - 1] < 0)
4782 i = synchro_nev(scs_imp, *kf, ierr);
4787 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4788 putevs(told, &i2, &ierr1);
4791 /* event conflict */
4804 /*--------------------------------------------------------------------------*/
4805 /* Subroutine cdoit */
4806 static void cdoit(double *told)
4808 /* propagation of continuous blocks outputs */
4809 /* Copyright INRIA */
4812 scicos_flag flag = 0;
4816 ScicosImport *scs_imp = NULL;
4819 scs_imp = getscicosimportptr();
4822 for (j = 0; j < * (scs_imp->ncord); j++)
4824 kf = &scs_imp->cord[j];
4825 C2F(curblk).kfun = *kf;
4826 /* continuous state */
4827 if (scs_imp->blocks[*kf - 1].nx > 0)
4829 scs_imp->blocks[*kf - 1].x = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4830 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4832 scs_imp->blocks[*kf - 1].nevprt = scs_imp->cord[j + * (scs_imp->ncord)];
4833 if (scs_imp->funtyp[*kf - 1] > -1)
4836 /*callf(told, xd, x, x,x,&flag);*/
4837 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4845 /* Initialize tvec */
4846 if (scs_imp->blocks[*kf - 1].nevout > 0)
4848 if (scs_imp->funtyp[*kf - 1] < 0)
4850 i = synchro_nev(scs_imp, *kf, ierr);
4855 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4856 putevs(told, &i2, &ierr1);
4859 /* event conflict */
4872 /*--------------------------------------------------------------------------*/
4873 /* Subroutine ddoit */
4874 static void ddoit(double *told)
4876 /* update states & event out on discrete activations */
4877 /* Copyright INRIA */
4880 scicos_flag flag = 0;
4882 int i = 0, i3 = 0, ierr1 = 0;
4883 int ii = 0, keve = 0;
4885 ScicosImport *scs_imp = NULL;
4888 scs_imp = getscicosimportptr();
4898 /* update continuous and discrete states on event */
4903 for (i = 0; i < kiwa; i++)
4906 if (critev[keve] != 0)
4910 i2 = scs_imp->ordptr[keve] - 1;
4911 for (ii = scs_imp->ordptr[keve - 1]; ii <= i2; ii++)
4913 kf = &scs_imp->ordclk[ii - 1];
4914 C2F(curblk).kfun = *kf;
4915 /* continuous state */
4916 if (scs_imp->blocks[*kf - 1].nx > 0)
4918 scs_imp->blocks[*kf - 1].x = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4919 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4922 scs_imp->blocks[*kf - 1].nevprt = scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1];
4924 if (scs_imp->blocks[*kf - 1].nevout > 0)
4926 if (scs_imp->funtyp[*kf - 1] >= 0)
4928 /* initialize evout */
4929 for (j = 0; j < scs_imp->blocks[*kf - 1].nevout; j++)
4931 scs_imp->blocks[*kf - 1].evout[j] = -1;
4935 if (scs_imp->blocks[*kf - 1].nevprt > 0) /* if event has continuous origin don't call*/
4937 /*callf(told, xd, x, x ,x,&flag);*/
4938 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4946 for (j = 0; j < scs_imp->blocks[*kf - 1].nevout; j++)
4948 if (scs_imp->blocks[*kf - 1].evout[j] >= 0.)
4950 i3 = j + scs_imp->clkptr[*kf - 1] ;
4951 addevs(scs_imp->blocks[*kf - 1].evout[j] + (*told), &i3, &ierr1);
4954 /* event conflict */
4963 if (scs_imp->blocks[*kf - 1].nevprt > 0)
4965 if (scs_imp->blocks[*kf - 1].nx + scs_imp->blocks[*kf - 1].nz + scs_imp->blocks[*kf - 1].noz > 0 || \
4966 *scs_imp->blocks[*kf - 1].work != NULL)
4968 /* if a hidden state exists, must also call (for new scope eg) */
4969 /* to avoid calling non-real activations */
4971 /*callf(told, xd, x, x,x,&flag);*/
4972 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4982 if (*scs_imp->blocks[*kf - 1].work != NULL)
4985 scs_imp->blocks[*kf - 1].nevprt = 0; /* in case some hidden continuous blocks need updating */
4986 /*callf(told, xd, x, x,x,&flag);*/
4987 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4998 /*--------------------------------------------------------------------------*/
4999 /* Subroutine edoit */
5000 static void edoit(double *told, int *kiwa)
5002 /* update blocks output on discrete activations */
5003 /* Copyright INRIA */
5006 scicos_flag flag = 0;
5007 int ierr1 = 0, i = 0;
5008 int kever = 0, ii = 0;
5010 ScicosImport *scs_imp = NULL;
5014 scs_imp = getscicosimportptr();
5019 *pointi = evtspt[kever];
5022 nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever - 1];
5029 for (ii = scs_imp->ordptr[kever - 1]; ii <= scs_imp->ordptr[kever] - 1; ii++)
5031 kf = &scs_imp->ordclk[ii - 1];
5032 C2F(curblk).kfun = *kf;
5034 if (scs_imp->funtyp[*kf - 1] > -1)
5036 /* continuous state */
5037 if (scs_imp->blocks[*kf - 1].nx > 0)
5039 scs_imp->blocks[*kf - 1].x = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
5040 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5043 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5046 /*callf(told, xd, x, x,x,&flag);*/
5047 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5055 /* Initialize tvec */
5056 if (scs_imp->blocks[*kf - 1].nevout > 0)
5058 if (scs_imp->funtyp[*kf - 1] < 0)
5060 i = synchro_nev(scs_imp, *kf, ierr);
5065 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5066 putevs(told, &i2, &ierr1);
5069 /* event conflict */
5082 /*--------------------------------------------------------------------------*/
5083 /* Subroutine odoit */
5084 static void odoit(double *told, double *xt, double *xtd, double *residual)
5086 /* update blocks derivative of continuous time block */
5087 /* Copyright INRIA */
5090 scicos_flag flag = 0;
5091 int keve = 0, kiwa = 0;
5092 int ierr1 = 0, i = 0;
5095 ScicosImport *scs_imp = NULL;
5098 scs_imp = getscicosimportptr();
5102 for (jj = 0; jj < * (scs_imp->noord); jj++)
5104 kf = &scs_imp->oord[jj];
5105 C2F(curblk).kfun = *kf;
5106 /* continuous state */
5107 if (scs_imp->blocks[*kf - 1].nx > 0)
5109 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
5110 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5111 scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5114 scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[jj + * (scs_imp->noord)];
5115 if (scs_imp->funtyp[*kf - 1] > -1)
5118 /*callf(told, xtd, xt, residual,x,&flag);*/
5119 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5127 if (scs_imp->blocks[*kf - 1].nevout > 0)
5129 if (scs_imp->funtyp[*kf - 1] < 0)
5131 if (scs_imp->blocks[*kf - 1].nmode > 0)
5133 i2 = scs_imp->blocks[*kf - 1].mode[0] + scs_imp->clkptr[*kf - 1] - 1;
5137 i = synchro_nev(scs_imp, *kf, ierr);
5142 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5144 putevs(told, &i2, &ierr1);
5147 /* event conflict */
5151 ozdoit(told, xt, xtd, &kiwa);
5160 /* update states derivatives */
5161 for (ii = 0; ii < * (scs_imp->noord); ii++)
5163 kf = &scs_imp->oord[ii];
5164 C2F(curblk).kfun = *kf;
5165 if (scs_imp->blocks[*kf - 1].nx > 0 || \
5166 *scs_imp->blocks[*kf - 1].work != NULL)
5168 /* work tests if a hidden state exists, used for delay block */
5170 /* continuous state */
5171 if (scs_imp->blocks[*kf - 1].nx > 0)
5173 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
5174 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5175 scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5177 scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[ii + * (scs_imp->noord)];
5178 /*callf(told, xtd, xt, residual,xt,&flag);*/
5179 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5189 for (i = 0; i < kiwa; i++)
5192 for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1; ii++)
5194 kf = &scs_imp->ordclk[ii - 1];
5195 C2F(curblk).kfun = *kf;
5196 if (scs_imp->blocks[*kf - 1].nx > 0 || \
5197 *scs_imp->blocks[*kf - 1].work != NULL)
5199 /* work tests if a hidden state exists */
5201 /* continuous state */
5202 if (scs_imp->blocks[*kf - 1].nx > 0)
5204 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
5205 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5206 scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5208 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5209 /*callf(told, xtd, xt, residual,xt,&flag);*/
5210 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5221 /*--------------------------------------------------------------------------*/
5222 /* Subroutine reinitdoit */
5223 static void reinitdoit(double *told)
5225 /* update blocks xproperties of continuous time block */
5226 /* Copyright INRIA */
5229 scicos_flag flag = 0;
5230 int keve = 0, kiwa = 0;
5231 int ierr1 = 0, i = 0;
5234 ScicosImport *scs_imp = NULL;
5237 scs_imp = getscicosimportptr();
5241 for (jj = 0; jj < * (scs_imp->noord); jj++)
5243 kf = &scs_imp->oord[jj];
5244 C2F(curblk).kfun = *kf;
5245 /* continuous state */
5246 if (scs_imp->blocks[*kf - 1].nx > 0)
5248 scs_imp->blocks[*kf - 1].x = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
5249 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5251 scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[jj + * (scs_imp->noord)];
5252 if (scs_imp->funtyp[*kf - 1] > -1)
5255 /*callf(told, xd, x, x,x,&flag);*/
5256 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5264 if (scs_imp->blocks[*kf - 1].nevout > 0 && scs_imp->funtyp[*kf - 1] < 0)
5266 i = synchro_nev(scs_imp, *kf, ierr);
5271 if (scs_imp->blocks[*kf - 1].nmode > 0)
5273 scs_imp->blocks[*kf - 1].mode[0] = i;
5275 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5276 putevs(told, &i2, &ierr1);
5279 /* event conflict */
5292 for (ii = 0; ii < * (scs_imp->noord); ii++)
5294 kf = &scs_imp->oord[ii];
5295 C2F(curblk).kfun = *kf;
5296 if (scs_imp->blocks[*kf - 1].nx > 0)
5299 scs_imp->blocks[*kf - 1].x = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
5300 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5301 scs_imp->blocks[*kf - 1].res = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5302 scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[ii + * (scs_imp->noord)];
5303 scs_imp->blocks[*kf - 1].xprop = &scs_imp->xprop[-1 + scs_imp->xptr[*kf - 1]];
5304 /*callf(told, xd, x, xd,x,&flag);*/
5305 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5314 for (i = 0; i < kiwa; i++)
5317 for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1; ii++)
5319 kf = &scs_imp->ordclk[ii - 1];
5320 C2F(curblk).kfun = *kf;
5321 if (scs_imp->blocks[*kf - 1].nx > 0)
5324 scs_imp->blocks[*kf - 1].x = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
5325 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5326 scs_imp->blocks[*kf - 1].res = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5327 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5328 scs_imp->blocks[*kf - 1].xprop = &scs_imp->xprop[-1 + scs_imp->xptr[*kf - 1]];
5329 /*callf(told, xd, x, xd,x,&flag);*/
5330 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5340 /*--------------------------------------------------------------------------*/
5341 /* Subroutine ozdoit */
5342 static void ozdoit(double *told, double *xt, double *xtd, int *kiwa)
5344 /* update blocks output of continuous time block on discrete activations */
5345 /* Copyright INRIA */
5348 scicos_flag flag = 0;
5350 int ierr1 = 0, i = 0;
5351 int ii = 0, kever = 0;
5353 ScicosImport *scs_imp = NULL;
5356 scs_imp = getscicosimportptr();
5360 *pointi = evtspt[kever];
5363 nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever - 1];
5371 for (ii = scs_imp->ordptr[kever - 1]; ii <= scs_imp->ordptr[kever] - 1; ii++)
5373 kf = &scs_imp->ordclk[ii - 1];
5374 C2F(curblk).kfun = *kf;
5375 if (scs_imp->funtyp[*kf - 1] > -1)
5377 /* continuous state */
5378 if (scs_imp->blocks[*kf - 1].nx > 0)
5380 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
5381 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5383 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5385 /*callf(told, xtd, xt, xt,x,&flag);*/
5386 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5394 /* Initialize tvec */
5395 if (scs_imp->blocks[*kf - 1].nevout > 0)
5397 if (scs_imp->funtyp[*kf - 1] < 0)
5399 if (phase == 1 || scs_imp->blocks[*kf - 1].nmode == 0)
5401 i = synchro_nev(scs_imp, *kf, ierr);
5409 i = scs_imp->blocks[*kf - 1].mode[0];
5411 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5412 putevs(told, &i2, &ierr1);
5415 /* event conflict */
5419 ozdoit(told, xt, xtd, kiwa);
5424 /*--------------------------------------------------------------------------*/
5425 /* Subroutine zdoit */
5426 static void zdoit(double *told, double *xt, double *xtd, double *g)
5428 /* update blocks zcross of continuous time block */
5429 /* Copyright INRIA */
5431 scicos_flag flag = 0;
5432 int keve = 0, kiwa = 0;
5433 int ierr1 = 0, i = 0, j = 0;
5436 ScicosImport *scs_imp = NULL;
5439 scs_imp = getscicosimportptr();
5442 for (i = 0; i < * (scs_imp->ng); i++)
5448 for (jj = 0; jj < * (scs_imp->nzord); jj++)
5450 kf = &scs_imp->zord[jj];
5451 C2F(curblk).kfun = *kf;
5452 /* continuous state */
5453 if (scs_imp->blocks[*kf - 1].nx > 0)
5455 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
5456 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5458 scs_imp->blocks[*kf - 1].nevprt = scs_imp->zord[jj + * (scs_imp->nzord)];
5460 if (scs_imp->funtyp[*kf - 1] > -1)
5463 /*callf(told, xtd, xt, xt,xt,&flag);*/
5464 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5472 /* Initialize tvec */
5473 if (scs_imp->blocks[*kf - 1].nevout > 0)
5475 if (scs_imp->funtyp[*kf - 1] < 0)
5477 if (phase == 1 || scs_imp->blocks[*kf - 1].nmode == 0)
5479 i = synchro_nev(scs_imp, *kf, ierr);
5487 i = scs_imp->blocks[*kf - 1].mode[0];
5489 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5490 putevs(told, &i2, &ierr1);
5493 /* event conflict */
5497 ozdoit(told, xt, xtd, &kiwa);
5506 /* update zero crossing surfaces */
5507 for (ii = 0; ii < * (scs_imp->nzord); ii++)
5509 kf = &scs_imp->zord[ii];
5510 C2F(curblk).kfun = *kf;
5511 if (scs_imp->blocks[*kf - 1].ng > 0)
5513 /* update g array ptr */
5514 scs_imp->blocks[*kf - 1].g = &g[scs_imp->zcptr[*kf - 1] - 1];
5515 if (scs_imp->funtyp[*kf - 1] > 0)
5518 /* continuous state */
5519 if (scs_imp->blocks[*kf - 1].nx > 0)
5521 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
5522 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5524 scs_imp->blocks[*kf - 1].nevprt = scs_imp->zord[ii + * (scs_imp->nzord)];
5525 /*callf(told, xtd, xt, xtd,g,&flag);*/
5526 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5535 j = synchro_g_nev(scs_imp, g, *kf, ierr);
5540 if ( (phase == 1) && (scs_imp->blocks[*kf - 1].nmode > 0) )
5542 scs_imp->blocks[*kf - 1].mode[0] = j;
5546 // scs_imp->blocks[*kf-1].g = &scs_imp->g[scs_imp->zcptr[*kf]-1];
5551 for (i = 0; i < kiwa; i++)
5554 for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1; ii++)
5556 kf = &scs_imp->ordclk[ii - 1];
5557 C2F(curblk).kfun = *kf;
5558 if (scs_imp->blocks[*kf - 1].ng > 0)
5560 /* update g array ptr */
5561 scs_imp->blocks[*kf - 1].g = &g[scs_imp->zcptr[*kf - 1] - 1];
5562 if (scs_imp->funtyp[*kf - 1] > 0)
5565 /* continuous state */
5566 if (scs_imp->blocks[*kf - 1].nx > 0)
5568 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
5569 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5571 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5572 /*callf(told, xtd, xt, xtd,g,&flag);*/
5573 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5582 j = synchro_g_nev(scs_imp, g, *kf, ierr);
5587 if ((phase == 1) && (scs_imp->blocks[*kf - 1].nmode > 0))
5589 scs_imp->blocks[*kf - 1].mode[0] = j;
5593 //scs_imp->blocks[*kf-1].g = &scs_imp->g[scs_imp->zcptr[*kf]-1];
5598 /*--------------------------------------------------------------------------*/
5599 /* Subroutine Jdoit */
5600 void Jdoit(double *told, double *xt, double *xtd, double *residual, int *job)
5602 /* update blocks jacobian of continuous time block */
5603 /* Copyright INRIA */
5606 scicos_flag flag = 0;
5607 int keve = 0, kiwa = 0;
5608 int ierr1 = 0, i = 0;
5611 ScicosImport *scs_imp = NULL;
5614 scs_imp = getscicosimportptr();
5618 for (jj = 0; jj < * (scs_imp->noord); jj++)
5620 kf = &scs_imp->oord[jj];
5621 C2F(curblk).kfun = *kf;
5622 scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[jj + * (scs_imp->noord)];
5623 if (scs_imp->funtyp[*kf - 1] > -1)
5626 /* applying desired output */
5627 if ((*job == 2) && (scs_imp->oord[jj] == AJacobian_block))
5631 /* continuous state */
5632 if (scs_imp->blocks[*kf - 1].nx > 0)
5634 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
5635 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5636 scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5639 /*callf(told, xtd, xt, residual,x,&flag);*/
5640 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5648 if (scs_imp->blocks[*kf - 1].nevout > 0)
5650 if (scs_imp->funtyp[*kf - 1] < 0)
5652 if (scs_imp->blocks[*kf - 1].nmode > 0)
5654 i2 = scs_imp->blocks[*kf - 1].mode[0] + scs_imp->clkptr[*kf - 1] - 1;
5658 i = synchro_nev(scs_imp, *kf, ierr);
5663 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5665 putevs(told, &i2, &ierr1);
5668 /* event conflict */
5672 ozdoit(told, xt, xtd, &kiwa);
5677 /* update states derivatives */
5678 for (ii = 0; ii < * (scs_imp->noord); ii++)
5680 kf = &scs_imp->oord[ii];
5681 C2F(curblk).kfun = *kf;
5682 if (scs_imp->blocks[*kf - 1].nx > 0 || \
5683 *scs_imp->blocks[*kf - 1].work != NULL)
5685 /* work tests if a hidden state exists, used for delay block */
5687 if (((*job == 1) && (scs_imp->oord[ii] == AJacobian_block)) || (*job != 1))
5693 /* continuous state */
5694 if (scs_imp->blocks[*kf - 1].nx > 0)
5696 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
5697 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5698 scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5700 scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[ii + * (scs_imp->noord)];
5701 /*callf(told, xtd, xt, residual,xt,&flag);*/
5702 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5712 for (i = 0; i < kiwa; i++)
5715 for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1; ii++)
5717 kf = &scs_imp->ordclk[ii - 1];
5718 C2F(curblk).kfun = *kf;
5719 if (scs_imp->blocks[*kf - 1].nx > 0 || \
5720 *scs_imp->blocks[*kf - 1].work != NULL)
5722 /* work tests if a hidden state exists */
5724 if (((*job == 1) && (scs_imp->oord[ii - 1] == AJacobian_block)) || (*job != 1))
5730 /* continuous state */
5731 if (scs_imp->blocks[*kf - 1].nx > 0)
5733 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
5734 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5735 scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5737 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5738 /*callf(told, xtd, xt, residual,xt,&flag);*/
5739 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5750 /*--------------------------------------------------------------------------*/
5751 /* Subroutine synchro_nev */
5752 static int synchro_nev(ScicosImport *scs_imp, int kf, int *ierr)
5754 /* synchro blocks computation */
5755 /* Copyright INRIA */
5756 SCSREAL_COP *outtbdptr = NULL; /*to store double of outtb*/
5757 SCSINT8_COP *outtbcptr = NULL; /*to store int8 of outtb*/
5758 SCSINT16_COP *outtbsptr = NULL; /*to store int16 of outtb*/
5759 SCSINT32_COP *outtblptr = NULL; /*to store int32 of outtb*/
5760 SCSUINT8_COP *outtbucptr = NULL; /*to store unsigned int8 of outtb */
5761 SCSUINT16_COP *outtbusptr = NULL; /*to store unsigned int16 of outtb */
5762 SCSUINT32_COP *outtbulptr = NULL; /*to store unsigned int32 of outtb */
5765 int i = 0; /* return 0 by default */
5767 /* variable for param */
5769 void **outtbptr = NULL;
5775 outtbtyp = scs_imp->outtbtyp;
5776 outtbptr = scs_imp->outtbptr;
5777 funtyp = scs_imp->funtyp;
5778 inplnk = scs_imp->inplnk;
5779 inpptr = scs_imp->inpptr;
5781 /* if-then-else blk */
5782 if (funtyp[kf - 1] == -1)
5784 switch (outtbtyp[-1 + inplnk[inpptr[kf - 1] - 1]])
5787 outtbdptr = (SCSREAL_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5788 cond = (*outtbdptr <= 0.);
5792 outtbdptr = (SCSCOMPLEX_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5793 cond = (*outtbdptr <= 0.);
5797 outtbcptr = (SCSINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5798 cond = (*outtbcptr <= 0);
5802 outtbsptr = (SCSINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5803 cond = (*outtbsptr <= 0);
5807 outtblptr = (SCSINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5808 cond = (*outtblptr <= 0);