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 /*--------------------------------------------------------------------------*/
41 /* Sundials includes */
42 #include <cvode/cvode.h> /* prototypes for CVODES fcts. and consts. */
43 #include <cvode/cvode_dense.h> /* prototype for CVDense */
44 #include <cvode/cvode_direct.h> /* prototypes for various DlsMat operations */
46 #include <ida/ida_dense.h>
47 #include <ida/ida_direct.h>
48 #include <nvector/nvector_serial.h> /* serial N_Vector types, fcts., and macros */
49 #include <sundials/sundials_dense.h> /* prototypes for various DlsMat operations */
50 #include <sundials/sundials_direct.h> /* definitions of DlsMat and DENSE_ELEM */
51 #include <sundials/sundials_types.h> /* definition of type realtype */
52 #include <sundials/sundials_math.h>
53 #include <kinsol/kinsol.h>
54 #include <kinsol/kinsol_dense.h>
55 #include <kinsol/kinsol_direct.h>
56 #include <sundials/sundials_extension.h> /* uses extension for scicos */
59 #include "machine.h" /* C2F */
60 #include "dynamic_link.h"
61 #include "scicos-def.h"
62 #include "stack-def.h"
67 #include "core_math.h"
68 #include "storeCommand.h"
71 #include "math_graphics.h"
79 #include "dynlib_scicos.h"
81 #include "lsodar.h" /* prototypes for lsodar fcts. and consts. */
83 #if defined(linux) && defined(__i386__)
84 #include "setPrecisionFPU.h"
87 #include "localization.h"
88 #include "charEncoding.h"
89 /*--------------------------------------------------------------------------*/
96 double *gwork; /* just added for a very special use: a
97 space passing to grblkdakr for zero crossing surfaces
98 when updating mode variables during initialization */
101 SCICOS_IMPEXP SCSPTR_struct C2F(scsptr);
102 /*--------------------------------------------------------------------------*/
105 if (*neq>0) {if (C2F(cmsolver).solver) CVodeFree(&ode_mem); else LSodarFree(&ode_mem);} \
106 if (*neq>0) N_VDestroy_Serial(y); \
107 if ( ng>0 ) FREE(jroot); \
108 if ( ng>0 ) FREE(zcros);
111 /* TJacque allocated by sundials */
113 if (*neq>0) free(TJacque); \
114 if (*neq>0) FREE(data->rwork); \
115 if (( ng>0 )&& (*neq>0)) FREE(data->gwork); \
116 if (*neq>0) N_VDestroy_Serial(data->ewt); \
117 if (*neq>0) FREE(data); \
118 if (*neq>0) {if (C2F(cmsolver).solver == 100) IDAFree(&dae_mem);} \
119 if (*neq>0) N_VDestroy_Serial(IDx); \
120 if (*neq>0) N_VDestroy_Serial(yp); \
121 if (*neq>0) N_VDestroy_Serial(yy); \
122 if ( ng>0 ) FREE(jroot); \
123 if ( ng>0 ) FREE(zcros); \
124 if (nmod>0) FREE(Mode_save);
126 #define freeouttbptr \
137 N_VDestroy_Serial(y); \
138 N_VDestroy_Serial(fscale); \
139 N_VDestroy_Serial(yscale); \
143 #define ONE RCONST(1.0)
144 #define ZERO RCONST(0.0)
145 #define T0 RCONST(0.0)
146 /*--------------------------------------------------------------------------*/
147 /* Table of constant values */
148 static int c__90 = 90;
149 static int c__91 = 91;
151 static double c_b14 = 0.;
154 int TCritWarning = 0;
156 static double *t0 = NULL, *tf = NULL;
157 static double *x = NULL, *tevts = NULL, *xd = NULL, *g = NULL;
158 static double Atol = 0., rtol = 0., ttol = 0., deltat = 0., hmax = 0.;
159 static int *ierr = NULL;
160 static int *pointi = NULL;
161 static int *xptr = NULL, *modptr = NULL, *evtspt = NULL;
162 static int *funtyp = NULL, *inpptr = NULL, *outptr = NULL, *inplnk = NULL, *outlnk = NULL;
163 static int *clkptr = NULL, *ordptr = NULL, *ordclk = NULL, *cord = NULL, *iord = NULL, *oord = NULL, *zord = NULL, *critev = NULL, *zcptr = NULL;
164 static int nblk = 0, nordptr = 0, nlnk = 0, nx = 0, ng = 0, ncord = 0, noord = 0, nzord = 0;
165 static int nordclk = 0, niord = 0, nmod = 0;
166 static int nelem = 0;
167 static int *mod = NULL;
168 static int *neq = NULL;
169 static int *xprop = NULL; /* xproperty */
170 static int debug_block = 0;
171 static int *iwa = NULL;
174 /* declaration of ptr for typed port */
175 static void **outtbptr = NULL; /*pointer array of object of outtb*/
176 static int *outtbsz = NULL; /*size of object of outtb*/
177 static int *outtbtyp = NULL; /*type of object of outtb*/
179 SCSREAL_COP *outtbdptr = NULL; /*to store double of outtb*/
180 SCSINT8_COP *outtbcptr = NULL; /*to store int8 of outtb*/
181 SCSINT16_COP *outtbsptr = NULL; /*to store int16 of outtb*/
182 SCSINT32_COP *outtblptr = NULL; /*to store int32 of outtb*/
183 SCSUINT8_COP *outtbucptr = NULL; /*to store unsigned int8 of outtb */
184 SCSUINT16_COP *outtbusptr = NULL; /*to store unsigned int16 of outtb */
185 SCSUINT32_COP *outtbulptr = NULL; /*to store unsigned int32 of outtb */
187 static outtb_el *outtb_elem = NULL;
189 static scicos_block *Blocks = NULL;
191 /* pass to external variable for code generation */
192 /* reserved variable name */
193 int *block_error = NULL;
194 double scicos_time = 0.;
196 int Jacobian_Flag = 0;
197 // double CI = 0., CJ = 0.; // doubles returned by Get_Jacobian_ci and Get_Jacobian_cj respectively
198 double CJJ = 0.; // returned by Get_Jacobian_parameter
199 double SQuround = 0.;
201 static int AJacobian_block = 0;
204 /* Variable declaration moved to scicos.c because it was in the scicos-def.h therefore
205 * multiple declaration of the variable and linkers were complaining about duplicate
208 SCICOS_IMPEXP COSDEBUGCOUNTER_struct C2F(cosdebugcounter);
209 SCICOS_IMPEXP RTFACTOR_struct C2F(rtfactor);
210 SCICOS_IMPEXP SOLVER_struct C2F(cmsolver);
211 SCICOS_IMPEXP CURBLK_struct C2F(curblk);
212 SCICOS_IMPEXP COSDEBUG_struct C2F(cosdebug);
213 SCICOS_IMPEXP COSHLT_struct C2F(coshlt);
214 SCICOS_IMPEXP DBCOS_struct C2F(dbcos);
215 SCICOS_IMPEXP COSTOL_struct C2F(costol);
216 SCICOS_IMPEXP COSERR_struct coserr;
217 /*--------------------------------------------------------------------------*/
218 static void FREE_blocks();
219 static void cosini(double *told);
220 static void cosend(double *told);
221 static void cossim(double *told);
222 static void cossimdaskr(double *told);
223 static void doit(double *told);
224 static void idoit(double *told);
225 static void zdoit(double *told, double *xt, double *xtd, double *g);
226 static void cdoit(double *told);
227 static void ozdoit(double *told, double *xt, double *xtd, int *kiwa);
228 static void odoit(double *told, double *xt, double *xtd, double *residual);
229 static void ddoit(double *told);
230 static void edoit(double *told, int *kiwa);
231 static void reinitdoit(double *told);
232 static int CallKinsol(double *told);
233 static int simblk(realtype t, N_Vector yy, N_Vector yp, void *f_data);
234 static int grblkdaskr(realtype t, N_Vector yy, N_Vector yp, realtype *gout, void *g_data);
235 static int grblk(realtype t, N_Vector yy, realtype *gout, void *g_data);
236 static int simblklsodar(int * nequations, realtype * tOld, realtype * actual, realtype * res);
237 static int grblklsodar(int * nequations, realtype * tOld, realtype * actual, int * ngroot, realtype * res);
238 static void addevs(double t, int *evtnb, int *ierr1);
239 static int synchro_g_nev(ScicosImport *scs_imp, double *g, int kf, int *ierr);
240 static void Multp(double *A, double *B, double *R, int ra, int rb, int ca, int cb);
241 static int read_id(ezxml_t *elements, char *id, double *value);
242 static int simblkdaskr(realtype tres, N_Vector yy, N_Vector yp, N_Vector resval, void *rdata);
243 static void SundialsErrHandler(int error_code, const char *module, const char *function, char *msg, void *user_data);
244 static int Jacobians(long int Neq, realtype tt, realtype cj, N_Vector yy,
245 N_Vector yp, N_Vector resvec, DlsMat Jacque, void *jdata,
246 N_Vector tempv1, N_Vector tempv2, N_Vector tempv3);
247 static void call_debug_scicos(scicos_block *block, scicos_flag *flag, int flagi, int deb_blk);
248 static int synchro_nev(ScicosImport *scs_imp, int kf, int *ierr);
249 /*--------------------------------------------------------------------------*/
250 extern int C2F(dset)(int *n, double *dx, double *dy, int *incy);
251 extern int C2F(dcopy)(int *, double *, int *, double *, int *);
252 extern int C2F(msgs)();
253 extern int C2F(getscsmax)();
254 extern int C2F(makescicosimport)();
255 extern int C2F(clearscicosimport)();
256 /*--------------------------------------------------------------------------*/
257 void putevs(double *t, int *evtnb, int *ierr1);
258 void Jdoit(double *told, double *xt, double *xtd, double *residual, int *job);
259 int simblkKinsol(N_Vector yy, N_Vector resval, void *rdata);
260 /*--------------------------------------------------------------------------*/
261 int C2F(scicos)(double *x_in, int *xptr_in, double *z__,
262 void **work, int *zptr, int *modptr_in,
263 void **oz, int *ozsz, int *oztyp, int *ozptr,
264 int *iz, int *izptr, double *t0_in,
265 double *tf_in, double *tevts_in, int *evtspt_in,
266 int *nevts, int *pointi_in, void **outtbptr_in,
267 int *outtbsz_in, int *outtbtyp_in,
268 outtb_el *outtb_elem_in, int *nelem1, int *nlnk1,
269 int *funptr, int *funtyp_in, int *inpptr_in,
270 int *outptr_in, int *inplnk_in, int *outlnk_in,
271 double *rpar, int *rpptr, int *ipar, int *ipptr,
272 void **opar, int *oparsz, int *opartyp, int *opptr,
273 int *clkptr_in, int *ordptr_in, int *nordptr1,
274 int *ordclk_in, int *cord_in, int *ncord1,
275 int *iord_in, int *niord1, int *oord_in,
276 int *noord1, int *zord_in, int *nzord1,
277 int *critev_in, int *nblk1, int *ztyp,
278 int *zcptr_in, int *subscr, int *nsubs,
279 double *simpar, int *flag__, int *ierr_out)
281 int i1, kf, lprt, in, out, job = 1;
284 static int mxtb = 0, ierr0 = 0, kfun0 = 0, i = 0, j = 0, k = 0, jj = 0;
285 static int ni = 0, no = 0;
286 static int nz = 0, noz = 0, nopar = 0;
289 // Set FPU Flag to Extended for scicos simulation
290 // in order to override Java setting it to Double.
291 #if defined(linux) && defined(__i386__)
294 /* Copyright INRIA */
295 /* iz,izptr are used to pass block labels */
302 /* Parameter adjustments */
306 modptr = modptr_in - 1;
310 evtspt = evtspt_in - 1;
311 tevts = tevts_in - 1;
312 outtbptr = outtbptr_in;
313 outtbsz = outtbsz_in;
314 outtbtyp = outtbtyp_in;
315 outtb_elem = outtb_elem_in;
316 funtyp = funtyp_in - 1;
317 inpptr = inpptr_in - 1;
318 outptr = outptr_in - 1;
319 inplnk = inplnk_in - 1;
320 outlnk = outlnk_in - 1;
324 clkptr = clkptr_in - 1;
325 ordptr = ordptr_in - 1;
326 ordclk = ordclk_in - 1;
332 critev = critev_in - 1;
334 zcptr = zcptr_in - 1;
342 C2F(rtfactor).scale = simpar[5];
343 C2F(cmsolver).solver = (int) simpar[6];
356 nordclk = ordptr[nordptr] - 1; /* number of rows in ordclk is ordptr(nclkp1)-1 */
357 ng = zcptr[nblk + 1] - 1; /* computes number of zero crossing surfaces */
358 nmod = modptr[nblk + 1] - 1; /* computes number of modes */
359 nz = zptr[nblk + 1] - 1; /* number of discrete real states */
360 noz = ozptr[nblk + 1] - 1; /* number of discrete object states */
361 nopar = opptr[nblk + 1] - 1; /* number of object parameters */
362 nx = xptr[nblk + 1] - 1; /* number of object parameters */
365 xd = &x[xptr[nblk + 1] - 1];
367 /* check for hard coded maxsize */
368 for (i = 1; i <= nblk; ++i)
370 if (funtyp[i] < 10000)
376 funtyp[i] = funtyp[i] % 1000 + 10000;
378 ni = inpptr[i + 1] - inpptr[i];
379 no = outptr[i + 1] - outptr[i];
384 /* hard coded maxsize in callf.c */
385 C2F(msgs)(&c__90, &c__0);
386 C2F(curblk).kfun = i;
391 else if (funtyp[i] == 2 || funtyp[i] == 3)
393 /* hard coded maxsize in scicos.h */
394 if (ni + no > SZ_SIZE)
396 C2F(msgs)(&c__90, &c__0);
397 C2F(curblk).kfun = i;
407 for (j = 1; j <= ni; ++j)
409 k = inplnk[inpptr[i] - 1 + j];
410 mxtb = mxtb + (outtbsz[k - 1] * outtbsz[(k - 1) + nlnk]);
415 for (j = 1; j <= no; ++j)
417 k = outlnk[outptr[i] - 1 + j];
418 mxtb = mxtb + (outtbsz[k - 1] * outtbsz[(k - 1) + nlnk]);
423 C2F(msgs)(&c__91, &c__0);
424 C2F(curblk).kfun = i;
431 if (nx > 0) /* xprop */
433 if ((xprop = MALLOC(sizeof(int) * nx)) == NULL )
439 for (i = 0; i < nx; i++) /* initialize */
443 if (nmod > 0) /* mod */
445 if ((mod = MALLOC(sizeof(int) * nmod)) == NULL )
455 if (ng > 0) /* g becomes global */
457 if ((g = MALLOC(sizeof(double) * ng)) == NULL )
472 debug_block = -1; /* no debug block for start */
473 C2F(cosdebugcounter).counter = 0;
475 /** Create Block's array **/
476 if ((Blocks = MALLOC(sizeof(scicos_block) * nblk)) == NULL )
494 /** Setting blocks properties for each entry in Block's array **/
496 /* 1 : type and pointer on simulation function */
497 for (kf = 0; kf < nblk; ++kf) /*for each block */
499 C2F(curblk).kfun = kf + 1;
501 Blocks[kf].type = funtyp[kf + 1];
504 switch (funtyp[kf + 1])
507 Blocks[kf].funpt = F2C(sciblk);
510 sciprint(_("type 1 function not allowed for scilab blocks\n"));
511 *ierr = 1000 + kf + 1;
515 sciprint(_("type 2 function not allowed for scilab blocks\n"));
516 *ierr = 1000 + kf + 1;
520 Blocks[kf].funpt = sciblk2;
524 Blocks[kf].funpt = sciblk4;
527 case 99: /* debugging block */
528 Blocks[kf].funpt = sciblk4;
529 /*Blocks[kf].type=4;*/
534 Blocks[kf].funpt = sciblk4;
535 Blocks[kf].type = 10004;
538 sciprint(_("Undefined Function type\n"));
539 *ierr = 1000 + kf + 1;
543 Blocks[kf].scsptr = -i; /* set scilab function adress for sciblk */
545 else if (i <= ntabsim)
547 Blocks[kf].funpt = *(tabsim[i - 1].fonc);
548 Blocks[kf].scsptr = 0; /* this is done for being able to test if a block
549 is a scilab block in the debugging phase when
555 GetDynFunc(i, &Blocks[kf].funpt);
556 if ( Blocks[kf].funpt == (voidf) 0)
558 sciprint(_("Function not found\n"));
559 *ierr = 1000 + kf + 1;
563 Blocks[kf].scsptr = 0; /* this is done for being able to test if a block
564 is a scilab block in the debugging phase when
568 /* 2 : Dimension properties */
569 Blocks[kf].ztyp = ztyp[kf + 1];
570 Blocks[kf].nx = xptr[kf + 2] - xptr[kf + 1]; /* continuuous state dimension*/
571 Blocks[kf].ng = zcptr[kf + 2] - zcptr[kf + 1]; /* number of zero crossing surface*/
572 Blocks[kf].nz = zptr[kf + 2] - zptr[kf + 1]; /* number of double discrete state*/
573 Blocks[kf].noz = ozptr[kf + 2] - ozptr[kf + 1]; /* number of other discrete state*/
574 Blocks[kf].nrpar = rpptr[kf + 2] - rpptr[kf + 1]; /* size of double precision parameter vector*/
575 Blocks[kf].nipar = ipptr[kf + 2] - ipptr[kf + 1]; /* size of integer precision parameter vector*/
576 Blocks[kf].nopar = opptr[kf + 2] - opptr[kf + 1]; /* number of other parameters (matrix, data structure,..)*/
577 Blocks[kf].nin = inpptr[kf + 2] - inpptr[kf + 1]; /* number of input ports */
578 Blocks[kf].nout = outptr[kf + 2] - outptr[kf + 1]; /* number of output ports */
580 /* 3 : input port properties */
581 /* in insz, we store :
582 * - insz[0..nin-1] : first dimension of input ports
583 * - insz[nin..2*nin-1] : second dimension of input ports
584 * - insz[2*nin..3*nin-1] : type of data of input ports
586 /* allocate size and pointer arrays (number of input ports)*/
587 Blocks[kf].insz = NULL;
588 Blocks[kf].inptr = NULL;
589 if (Blocks[kf].nin != 0)
591 if ((Blocks[kf].insz = MALLOC(Blocks[kf].nin * 3 * sizeof(int))) == NULL )
597 if ((Blocks[kf].inptr = MALLOC(Blocks[kf].nin * sizeof(double*))) == NULL )
604 for (in = 0; in < Blocks[kf].nin; in++)
606 lprt = inplnk[inpptr[kf + 1] + in];
607 Blocks[kf].inptr[in] = outtbptr[lprt - 1]; /* pointer on the data*/
608 Blocks[kf].insz[in] = outtbsz[lprt - 1]; /* row dimension of the input port*/
609 Blocks[kf].insz[Blocks[kf].nin + in] = outtbsz[(lprt - 1) + nlnk]; /* column dimension of the input port*/
610 Blocks[kf].insz[2 * Blocks[kf].nin + in] = outtbtyp[lprt - 1]; /*type of data of the input port*/
612 /* 4 : output port properties */
613 /* in outsz, we store :
614 * - outsz[0..nout-1] : first dimension of output ports
615 * - outsz[nout..2*nout-1] : second dimension of output ports
616 * - outsz[2*nout..3*nout-1] : type of data of output ports
618 /* allocate size and pointer arrays (number of output ports)*/
619 Blocks[kf].outsz = NULL;
620 Blocks[kf].outptr = NULL;
621 if (Blocks[kf].nout != 0)
623 if ((Blocks[kf].outsz = MALLOC(Blocks[kf].nout * 3 * sizeof(int))) == NULL )
629 if ((Blocks[kf].outptr = MALLOC(Blocks[kf].nout * sizeof(double*))) == NULL )
637 for (out = 0; out < Blocks[kf].nout; out++) /*for each output port */
639 lprt = outlnk[outptr[kf + 1] + out];
640 Blocks[kf].outptr[out] = outtbptr[lprt - 1]; /*pointer on data */
641 Blocks[kf].outsz[out] = outtbsz[lprt - 1]; /*row dimension of output port*/
642 Blocks[kf].outsz[Blocks[kf].nout + out] = outtbsz[(lprt - 1) + nlnk]; /*column dimension of output ports*/
643 Blocks[kf].outsz[2 * Blocks[kf].nout + out] = outtbtyp[lprt - 1]; /*type of data of output port */
646 /* 5 : event output port properties */
647 Blocks[kf].evout = NULL;
648 Blocks[kf].nevout = clkptr[kf + 2] - clkptr[kf + 1];
649 if (Blocks[kf].nevout != 0)
651 if ((Blocks[kf].evout = CALLOC(Blocks[kf].nevout, sizeof(double))) == NULL )
659 /* 6 : pointer on the begining of the double discrete state array ( z) */
660 Blocks[kf].z = &(z__[zptr[kf + 1] - 1]);
662 /* 7 : type, size and pointer on the other discrete states data structures (oz) */
663 Blocks[kf].ozsz = NULL;
664 if (Blocks[kf].noz == 0)
666 Blocks[kf].ozptr = NULL;
667 Blocks[kf].oztyp = NULL;
671 Blocks[kf].ozptr = &(oz[ozptr[kf + 1] - 1]);
672 if ((Blocks[kf].ozsz = MALLOC(Blocks[kf].noz * 2 * sizeof(int))) == NULL )
678 for (i = 0; i < Blocks[kf].noz; i++)
680 Blocks[kf].ozsz[i] = ozsz[(ozptr[kf + 1] - 1) + i];
681 Blocks[kf].ozsz[i + Blocks[kf].noz] = ozsz[(ozptr[kf + 1] - 1 + noz) + i];
683 Blocks[kf].oztyp = &(oztyp[ozptr[kf + 1] - 1]);
686 /* 8 : pointer on the begining of the double parameter array ( rpar ) */
687 Blocks[kf].rpar = &(rpar[rpptr[kf + 1] - 1]);
689 /* 9 : pointer on the begining of the integer parameter array ( ipar ) */
690 Blocks[kf].ipar = &(ipar[ipptr[kf + 1] - 1]);
692 /* 10 : type, size and pointer on the other parameters data structures (opar) */
693 Blocks[kf].oparsz = NULL;
694 if (Blocks[kf].nopar == 0)
696 Blocks[kf].oparptr = NULL;
697 Blocks[kf].opartyp = NULL;
701 Blocks[kf].oparptr = &(opar[opptr[kf + 1] - 1]);
702 if ((Blocks[kf].oparsz = MALLOC(Blocks[kf].nopar * 2 * sizeof(int))) == NULL)
708 for (i = 0; i < Blocks[kf].nopar; i++)
710 Blocks[kf].oparsz[i] = oparsz[(opptr[kf + 1] - 1) + i];
711 Blocks[kf].oparsz[i + Blocks[kf].nopar] = oparsz[(opptr[kf + 1] - 1 + nopar) + i];
713 Blocks[kf].opartyp = &(opartyp[opptr[kf + 1] - 1]);
716 /* 10 : pointer on the beginning of the residual array (res) */
717 Blocks[kf].res = NULL;
718 if (Blocks[kf].nx != 0)
720 if ((Blocks[kf].res = MALLOC(Blocks[kf].nx * sizeof(double))) == NULL)
728 /* 11 : block label (label) */
729 i1 = izptr[kf + 2] - izptr[kf + 1];
730 if ((Blocks[kf].label = MALLOC(sizeof(char) * (i1 + 1))) == NULL)
736 Blocks[kf].label[i1] = '\0';
737 C2F(cvstr)(&i1, &(iz[izptr[kf + 1] - 1]), Blocks[kf].label, &job, i1);
739 /* 12 : block array of crossed surfaces (jroot) */
740 Blocks[kf].jroot = NULL;
741 if (Blocks[kf].ng != 0)
743 if ((Blocks[kf].jroot = CALLOC(Blocks[kf].ng, sizeof(int))) == NULL)
751 /* 13 : block work array (work) */
752 Blocks[kf].work = (void **)(((double *)work) + kf);
754 /* 14 : block modes array (mode) */
755 Blocks[kf].nmode = modptr[kf + 2] - modptr[kf + 1];
756 if (Blocks[kf].nmode != 0)
758 Blocks[kf].mode = &(mod[modptr[kf + 1] - 1]);
761 /* 15 : block xprop array (xprop) */
762 Blocks[kf].xprop = NULL;
763 if (Blocks[kf].nx != 0)
765 Blocks[kf].xprop = &(xprop[xptr[kf + 1] - 1]);
768 /* 16 : pointer on the zero crossing surface computation function of the block (g) */
770 if (Blocks[kf].ng != 0)
772 Blocks[kf].g = &(g[zcptr[kf + 1] - 1]);
775 /** all block properties are stored in the Blocks array **/
781 if ((iwa = MALLOC(sizeof(int) * (*nevts))) == NULL)
789 /* save ptr of scicos in import structure */
790 C2F(makescicosimport)(x, &nx, &xptr[1], &zcptr[1], z__, &nz, &zptr[1],
791 &noz, oz, ozsz, oztyp, &ozptr[1],
792 g, &ng, mod, &nmod, &modptr[1], iz, &izptr[1],
793 &inpptr[1], &inplnk[1], &outptr[1], &outlnk[1],
794 outtbptr, outtbsz, outtbtyp,
796 &nlnk, rpar, &rpptr[1], ipar, &ipptr[1],
797 opar, oparsz, opartyp, &opptr[1],
798 &nblk, subscr, nsubs,
799 &tevts[1], &evtspt[1], nevts, pointi,
800 &iord[1], &niord, &oord[1], &noord, &zord[1], &nzord,
801 funptr, &funtyp[1], &ztyp[1],
802 &cord[1], &ncord, &ordclk[1], &nordclk, &clkptr[1],
803 &ordptr[1], &nordptr, &critev[1], iwa, Blocks,
804 t0, tf, &Atol, &rtol, &ttol, &deltat, &hmax,
807 if (*flag__ == 1) /*start*/
809 /* blocks initialization */
810 for (kf = 0; kf < nblk; ++kf)
812 *(Blocks[kf].work) = NULL;
818 kfun0 = C2F(curblk).kfun;
821 C2F(curblk).kfun = kfun0;
825 else if (*flag__ == 2) /*run*/
829 switch (C2F(cmsolver).solver)
831 case 0: // LSodar - Dynamic / Dynamic
832 case 1: // CVode - BDF / Newton
833 case 2: // CVode - BDF / Functional
834 case 3: // CVode - Adams / Newton
835 case 4: // CVode - Adams / Functional
836 case 5: // Dormand-Prince
837 case 6: // Runge-Kutta
838 case 7: // Implicit Runge-Kutta - RK / Fixed-Point
844 default: // Unknown solver number
851 kfun0 = C2F(curblk).kfun;
854 C2F(curblk).kfun = kfun0;
858 else if (*flag__ == 3) /*finish*/
863 else if (*flag__ == 4) /*linear*/
869 if ((W = MALLOC(sizeof(double) * (Max(nx, ng)))) == NULL )
877 /*---------instead of old simblk--------*/
878 /* C2F(simblk)(&nx, t0, x, W); */
880 if (ng > 0 && nmod > 0)
882 zdoit(t0, x, x + nx, W); /* updating modes as a function of state values; this was necessary in iGUI*/
884 for (jj = 0; jj < nx; jj++)
888 C2F(ierode).iero = 0;
890 if (C2F(cmsolver).solver < 100)
896 odoit(t0, x, x + nx, W);
898 C2F(ierode).iero = *ierr;
899 /*-----------------------------------------*/
900 for (i = 0; i < nx; ++i)
907 else if (*flag__ == 5) /* initial_KINSOL= "Kinsol" */
909 C2F(ierode).iero = 0;
913 *ierr = C2F(ierode).iero;
920 C2F(clearscicosimport)();
923 /*--------------------------------------------------------------------------*/
925 static int check_flag(void *flagvalue, char *funcname, int opt)
929 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
930 if (opt == 0 && flagvalue == NULL)
932 sciprint(_("\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n"), funcname);
935 /* Check if flag < 0 */
938 errflag = (int *) flagvalue;
941 sciprint(_("\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n"),
946 /* Check if function returned NULL pointer - no memory allocated */
947 else if (opt == 2 && flagvalue == NULL)
949 sciprint(_("\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n"), funcname);
956 /*--------------------------------------------------------------------------*/
957 static void cosini(double *told)
959 static scicos_flag flag__ = 0;
962 static int kfune = 0;
965 SCSREAL_COP *outtbd = NULL; /*to save double of outtb*/
966 SCSINT8_COP *outtbc = NULL; /*to save int8 of outtb*/
967 SCSINT16_COP *outtbs = NULL; /*to save int16 of outtb*/
968 SCSINT32_COP *outtbl = NULL; /*to save int32 of outtb*/
969 SCSUINT8_COP *outtbuc = NULL; /*to save unsigned int8 of outtb*/
970 SCSUINT16_COP *outtbus = NULL; /*to save unsigned int16 of outtb*/
971 SCSUINT32_COP *outtbul = NULL; /*to save unsigned int32 of outtb*/
972 int szouttbd = 0; /*size of arrays*/
973 int szouttbc = 0, szouttbs = 0, szouttbl = 0;
974 int szouttbuc = 0, szouttbus = 0, szouttbul = 0;
975 int curouttbd = 0; /*current position in arrays*/
976 int curouttbc = 0, curouttbs = 0, curouttbl = 0;
977 int curouttbuc = 0, curouttbus = 0, curouttbul = 0;
979 int ii = 0, kk = 0; /*local counters*/
980 int sszz = 0; /*local size of element of outtb*/
981 /*Allocation of arrays for outtb*/
982 for (ii = 0; ii < nlnk; ii++)
984 switch (outtbtyp[ii])
987 szouttbd += outtbsz[ii] * outtbsz[ii + nlnk]; /*double real matrix*/
988 outtbd = (SCSREAL_COP *) REALLOC (outtbd, szouttbd * sizeof(SCSREAL_COP));
992 szouttbd += 2 * outtbsz[ii] * outtbsz[ii + nlnk]; /*double complex matrix*/
993 outtbd = (SCSCOMPLEX_COP *) REALLOC (outtbd, szouttbd * sizeof(SCSCOMPLEX_COP));
997 szouttbc += outtbsz[ii] * outtbsz[ii + nlnk]; /*int8*/
998 outtbc = (SCSINT8_COP *) REALLOC (outtbc, szouttbc * sizeof(SCSINT8_COP));
1002 szouttbs += outtbsz[ii] * outtbsz[ii + nlnk]; /*int16*/
1003 outtbs = (SCSINT16_COP *) REALLOC (outtbs, szouttbs * sizeof(SCSINT16_COP));
1007 szouttbl += outtbsz[ii] * outtbsz[ii + nlnk]; /*int32*/
1008 outtbl = (SCSINT32_COP *) REALLOC (outtbl, szouttbl * sizeof(SCSINT32_COP));
1012 szouttbuc += outtbsz[ii] * outtbsz[ii + nlnk]; /*uint8*/
1013 outtbuc = (SCSUINT8_COP *) REALLOC (outtbuc, szouttbuc * sizeof(SCSUINT8_COP));
1017 szouttbus += outtbsz[ii] * outtbsz[ii + nlnk]; /*uint16*/
1018 outtbus = (SCSUINT16_COP *) REALLOC (outtbus, szouttbus * sizeof(SCSUINT16_COP));
1022 szouttbul += outtbsz[ii] * outtbsz[ii + nlnk]; /*uint32*/
1023 outtbul = (SCSUINT32_COP *) REALLOC (outtbul, szouttbul * sizeof(SCSUINT32_COP));
1026 default : /* Add a message here */
1032 AJacobian_block = 0;
1037 /* initialization (flag 4) */
1038 /* loop on blocks */
1039 C2F(dset)(&ng, &c_b14, g, &c__1);
1041 for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
1044 if (Blocks[C2F(curblk).kfun - 1].nx > 0)
1046 Blocks[C2F(curblk).kfun - 1].x = &x[xptr[C2F(curblk).kfun] - 1];
1047 Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
1049 Blocks[C2F(curblk).kfun - 1].nevprt = 0;
1050 if (funtyp[C2F(curblk).kfun] >= 0) /* debug_block is not called here */
1052 /*callf(told, xd, x, x,g,&flag__);*/
1054 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1055 if (flag__ < 0 && *ierr == 0)
1058 kfune = C2F(curblk).kfun;
1060 if ((Jacobian_Flag == 1) && (AJacobian_block == 0))
1062 AJacobian_block = C2F(curblk).kfun;
1068 C2F(curblk).kfun = kfune;
1073 /* initialization (flag 6) */
1075 for (jj = 1; jj <= ncord; ++jj)
1077 C2F(curblk).kfun = cord[jj];
1078 Blocks[C2F(curblk).kfun - 1].nevprt = 0;
1079 if (funtyp[C2F(curblk).kfun] >= 0)
1081 /*callf(told, xd, x, x,g,&flag__);*/
1082 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1092 /* point-fix iterations */
1094 for (i = 1; i <= nblk + 1; ++i) /*for each block*/
1096 /* loop on blocks */
1097 for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
1099 Blocks[C2F(curblk).kfun - 1].nevprt = 0;
1100 if (funtyp[C2F(curblk).kfun] >= 0)
1102 /*callf(told, xd, x, x,g,&flag__);*/
1103 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1114 for (jj = 1; jj <= ncord; ++jj) /*for each continous block*/
1116 C2F(curblk).kfun = cord[jj];
1117 if (funtyp[C2F(curblk).kfun] >= 0)
1119 /*callf(told, xd, x, x,g,&flag__);*/
1120 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1130 /*comparison between outtb and arrays*/
1138 for (jj = 0; jj < nlnk; jj++)
1140 switch (outtbtyp[jj]) /*for each type of ports*/
1143 outtbdptr = (SCSREAL_COP *)outtbptr[jj]; /*double real matrix*/
1144 sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1145 for (kk = 0; kk < sszz; kk++)
1147 if (outtbdptr[kk] != (SCSREAL_COP)outtbd[curouttbd + kk])
1156 outtbdptr = (SCSCOMPLEX_COP *)outtbptr[jj]; /*double complex matrix*/
1157 sszz = 2 * outtbsz[jj] * outtbsz[jj + nlnk];
1158 for (kk = 0; kk < sszz; kk++)
1160 if (outtbdptr[kk] != (SCSCOMPLEX_COP)outtbd[curouttbd + kk])
1169 outtbcptr = (SCSINT8_COP *)outtbptr[jj]; /*int8*/
1170 sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1171 for (kk = 0; kk < sszz; kk++)
1173 if (outtbcptr[kk] != (SCSINT8_COP)outtbc[curouttbc + kk])
1182 outtbsptr = (SCSINT16_COP *)outtbptr[jj]; /*int16*/
1183 sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1184 for (kk = 0; kk < sszz; kk++)
1186 if (outtbsptr[kk] != (SCSINT16_COP)outtbs[curouttbs + kk])
1195 outtblptr = (SCSINT32_COP *)outtbptr[jj]; /*int32*/
1196 sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1197 for (kk = 0; kk < sszz; kk++)
1199 if (outtblptr[kk] != (SCSINT32_COP)outtbl[curouttbl + kk])
1208 outtbucptr = (SCSUINT8_COP *)outtbptr[jj]; /*uint8*/
1209 sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1210 for (kk = 0; kk < sszz; kk++)
1212 if (outtbucptr[kk] != (SCSUINT8_COP)outtbuc[curouttbuc + kk])
1221 outtbusptr = (SCSUINT16_COP *)outtbptr[jj]; /*uint16*/
1222 sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1223 for (kk = 0; kk < sszz; kk++)
1225 if (outtbusptr[kk] != (SCSUINT16_COP)outtbus[curouttbus + kk])
1234 outtbulptr = (SCSUINT32_COP *)outtbptr[jj]; /*uint32*/
1235 sszz = outtbsz[jj] * outtbsz[jj + nlnk];
1236 for (kk = 0; kk < sszz; kk++)
1238 if (outtbulptr[kk] != (SCSUINT32_COP)outtbul[curouttbul + kk])
1246 default : /* Add a message here */
1254 /*Save data of outtb in arrays*/
1262 for (ii = 0; ii < nlnk; ii++) /*for each link*/
1264 switch (outtbtyp[ii]) /*switch to type of outtb object*/
1267 outtbdptr = (SCSREAL_COP *)outtbptr[ii]; /*double real matrix*/
1268 sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1269 C2F(dcopy)(&sszz, outtbdptr, &c__1, &outtbd[curouttbd], &c__1);
1274 outtbdptr = (SCSCOMPLEX_COP *)outtbptr[ii]; /*double complex matrix*/
1275 sszz = 2 * outtbsz[ii] * outtbsz[ii + nlnk];
1276 C2F(dcopy)(&sszz, outtbdptr, &c__1, &outtbd[curouttbd], &c__1);
1281 outtbcptr = (SCSINT8_COP *)outtbptr[ii]; /*int8*/
1282 sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1283 for (kk = 0; kk < sszz; kk++)
1285 outtbc[curouttbc + kk] = (SCSINT8_COP)outtbcptr[kk];
1291 outtbsptr = (SCSINT16_COP *)outtbptr[ii]; /*int16*/
1292 sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1293 for (kk = 0; kk < sszz; kk++)
1295 outtbs[curouttbs + kk] = (SCSINT16_COP)outtbsptr[kk];
1301 outtblptr = (SCSINT32_COP *)outtbptr[ii]; /*int32*/
1302 sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1303 for (kk = 0; kk < sszz; kk++)
1305 outtbl[curouttbl + kk] = (SCSINT32_COP)outtblptr[kk];
1311 outtbucptr = (SCSUINT8_COP *)outtbptr[ii]; /*uint8*/
1312 sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1313 for (kk = 0; kk < sszz; kk++)
1315 outtbuc[curouttbuc + kk] = (SCSUINT8_COP)outtbucptr[kk];
1321 outtbusptr = (SCSUINT16_COP *)outtbptr[ii]; /*uint16*/
1322 sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1323 for (kk = 0; kk < sszz; kk++)
1325 outtbus[curouttbus + kk] = (SCSUINT16_COP)outtbusptr[kk];
1331 outtbulptr = (SCSUINT32_COP *)outtbptr[ii]; /*uint32*/
1332 sszz = outtbsz[ii] * outtbsz[ii + nlnk];
1333 for (kk = 0; kk < sszz; kk++)
1335 outtbul[curouttbul + kk] = (SCSUINT32_COP)outtbulptr[kk];
1340 default : /* Add a message here */
1349 /*--------------------------------------------------------------------------*/
1350 static void cossim(double *told)
1352 /* System generated locals */
1355 //** used for the [stop] button
1356 static char CommandToUnstack[1024];
1357 static int CommandLength = 0;
1358 static int SeqSync = 0;
1361 /* Local variables */
1362 static scicos_flag flag__ = 0;
1363 static int ierr1 = 0;
1364 static int j = 0, k = 0;
1365 static double t = 0.;
1367 static double rhotmp = 0., tstop = 0.;
1368 static int inxsci = 0;
1369 static int kpo = 0, kev = 0;
1370 int Discrete_Jump = 0;
1371 int *jroot = NULL, *zcros = NULL;
1372 realtype reltol = 0., abstol = 0.;
1374 void *ode_mem = NULL;
1375 int flag = 0, flagr = 0;
1377 // Defining function pointers, for more readability
1378 int (* ODE) (void*, realtype, N_Vector, realtype*, int);
1379 int (* ODEReInit) (void*, realtype, N_Vector);
1380 int (* ODESetMaxStep) (void*, realtype);
1381 int (* ODESetStopTime) (void*, realtype);
1382 int (* ODEGetRootInfo) (void*, int*);
1383 int (* ODESStolerances) (void*, realtype, realtype);
1384 /* Generic flags for stop mode */
1385 int ODE_NORMAL = 1; /* ODE_NORMAL = CV_NORMAL = LS_NORMAL = 1 */
1386 int ODE_ONE_STEP = 2; /* ODE_ONE_STEP = CV_ONE_STEP = LS_ONE_STEP = 2 */
1387 switch (C2F(cmsolver).solver)
1391 ODEReInit = &LSodarReInit;
1392 ODESetMaxStep = &LSodarSetMaxStep;
1393 ODESetStopTime = &LSodarSetStopTime;
1394 ODEGetRootInfo = &LSodarGetRootInfo;
1395 ODESStolerances = &LSodarSStolerances;
1397 case 1: // CVode BDF / Newton
1398 case 2: // CVode BDF / Functional
1399 case 3: // CVode Adams / Newton
1400 case 4: // CVode Adams / Functional
1401 case 5: // Dormand-Prince
1402 case 6: // Runge-Kutta
1403 case 7: // Implicit Runge-Kutta
1405 ODEReInit = &CVodeReInit;
1406 ODESetMaxStep = &CVodeSetMaxStep;
1407 ODESetStopTime = &CVodeSetStopTime;
1408 ODEGetRootInfo = &CVodeGetRootInfo;
1409 ODESStolerances = &CVodeSStolerances;
1411 default: // Unknown solver number
1419 if ((jroot = MALLOC(sizeof(int) * ng)) == NULL )
1426 for ( jj = 0 ; jj < ng ; jj++ )
1434 if ((zcros = MALLOC(sizeof(int) * ng)) == NULL )
1445 reltol = (realtype) rtol;
1446 abstol = (realtype) Atol; /* Ith(abstol,1) = realtype) Atol;*/
1448 if (*neq > 0) /* Unfortunately CVODE does not work with NEQ==0 */
1450 y = N_VNewEmpty_Serial(*neq);
1451 if (check_flag((void *)y, "N_VNewEmpty_Serial", 0))
1469 /* Set extension of Sundials for scicos */
1470 set_sundials_with_extension(TRUE);
1472 switch (C2F(cmsolver).solver)
1475 ode_mem = LSodarCreate(neq, ng); /* Create the lsodar problem */
1478 ode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
1481 ode_mem = CVodeCreate(CV_BDF, CV_FUNCTIONAL);
1484 ode_mem = CVodeCreate(CV_ADAMS, CV_NEWTON);
1487 ode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);
1490 ode_mem = CVodeCreate(CV_DOPRI, CV_FUNCTIONAL);
1493 ode_mem = CVodeCreate(CV_ExpRK, CV_FUNCTIONAL);
1496 ode_mem = CVodeCreate(CV_ImpRK, CV_FUNCTIONAL);
1500 /* ode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);*/
1502 if (check_flag((void *)ode_mem, "CVodeCreate", 0))
1505 N_VDestroy_Serial(y);
1511 if (C2F(cmsolver).solver == 0)
1513 flag = LSodarSetErrHandlerFn(ode_mem, SundialsErrHandler, NULL);
1517 flag = CVodeSetErrHandlerFn(ode_mem, SundialsErrHandler, NULL);
1519 if (check_flag(&flag, "CVodeSetErrHandlerFn", 1))
1521 *ierr = 300 + (-flag);
1526 if (C2F(cmsolver).solver == 0)
1528 flag = LSodarInit (ode_mem, simblklsodar, T0, y);
1532 flag = CVodeInit (ode_mem, simblk, T0, y);
1534 if (check_flag(&flag, "CVodeInit", 1))
1536 *ierr = 300 + (-flag);
1541 flag = ODESStolerances(ode_mem, reltol, abstol);
1542 if (check_flag(&flag, "CVodeSStolerances", 1))
1544 *ierr = 300 + (-flag);
1549 if (C2F(cmsolver).solver == 0)
1551 flag = LSodarRootInit(ode_mem, ng, grblklsodar);
1555 flag = CVodeRootInit(ode_mem, ng, grblk);
1557 if (check_flag(&flag, "CVodeRootInit", 1))
1559 *ierr = 300 + (-flag);
1564 if (C2F(cmsolver).solver != 0)
1565 /* Call CVDense to specify the CVDENSE dense linear solver */
1567 flag = CVDense(ode_mem, *neq);
1569 if (check_flag(&flag, "CVDense", 1))
1571 *ierr = 300 + (-flag);
1578 flag = ODESetMaxStep(ode_mem, (realtype) hmax);
1579 if (check_flag(&flag, "CVodeSetMaxStep", 1))
1581 *ierr = 300 + (-flag);
1586 /* Set the Jacobian routine to Jac (user-supplied)
1587 flag = CVDlsSetDenseJacFn(ode_mem, Jac);
1588 if (check_flag(&flag, "CVDlsSetDenseJacFn", 1)) return(1); */
1590 }/* testing if neq>0 */
1593 C2F(coshlt).halt = 0;
1596 C2F(xscion)(&inxsci);
1597 /* initialization */
1598 C2F(realtimeinit)(told, &C2F(rtfactor).scale);
1604 for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
1606 if (Blocks[C2F(curblk).kfun - 1].ng >= 1)
1608 zcros[jj] = C2F(curblk).kfun;
1612 /* . ng >= jj required */
1617 /* initialization (propagation of constant blocks outputs) */
1624 /*--discrete zero crossings----dzero--------------------*/
1625 if (ng > 0) /* storing ZC signs just after a solver call*/
1627 /*zdoit(told, g, x, x);*/
1628 zdoit(told, x, x, g);
1634 for (jj = 0; jj < ng; ++jj)
1644 /*--discrete zero crossings----dzero--------------------*/
1646 /* main loop on time */
1649 while (ismenu()) //** if the user has done something, do the actions
1652 SeqSync = GetCommand(CommandToUnstack); //** get to the action
1653 CommandLength = (int)strlen(CommandToUnstack);
1654 syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
1656 if (C2F(coshlt).halt != 0)
1658 if (C2F(coshlt).halt == 2)
1660 *told = *tf; /* end simulation */
1662 C2F(coshlt).halt = 0;
1674 if (fabs(t - *told) < ttol)
1677 /* update output part */
1681 /* ! scheduling problem */
1688 if (xptr[nblk + 1] == 1)
1690 /* . no continuous state */
1691 if (*told + deltat + ttol > t)
1699 /* . update outputs of 'c' type blocks with no continuous state */
1702 /* . we are at the end, update continuous part before leaving */
1714 rhotmp = *tf + ttol;
1719 if (critev[kpo] == 1)
1721 rhotmp = tevts[kpo];
1736 t = Min(*told + deltat, Min(t, *tf + ttol));
1738 if (ng > 0 && hot == 0 && nmod > 0)
1740 zdoit(told, x, x, g);
1748 if (hot == 0) /* hot==0 : cold restart*/
1750 flag = ODESetStopTime(ode_mem, (realtype)tstop); /* Setting the stop time*/
1751 if (check_flag(&flag, "CVodeSetStopTime", 1))
1753 *ierr = 300 + (-flag);
1758 flag = ODEReInit(ode_mem, (realtype)(*told), y);
1759 if (check_flag(&flag, "CVodeReInit", 1))
1761 *ierr = 300 + (-flag);
1767 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1769 sciprint(_("****SUNDIALS.Cvode from: %f to %f hot= %d \n"), *told, t, hot);
1772 /*--discrete zero crossings----dzero--------------------*/
1773 /*--check for Dzeros after Mode settings or ddoit()----*/
1776 if (ng > 0 && hot == 0)
1778 zdoit(told, x, x, g);
1784 for (jj = 0; jj < ng; ++jj)
1786 if ((g[jj] >= 0.0) && (jroot[jj] == -5))
1791 else if ((g[jj] < 0.0) && (jroot[jj] == 5))
1802 /*--discrete zero crossings----dzero--------------------*/
1804 if (Discrete_Jump == 0) /* if there was a dzero, its event should be activated*/
1807 flag = ODE(ode_mem, t, y, told, ODE_NORMAL);
1817 flag = CV_ROOT_RETURN; /* in order to handle discrete jumps */
1820 /* . update outputs of 'c' type blocks if we are at the end*/
1833 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1835 sciprint(_("****SUNDIALS.Cvode reached: %f\n"), *told);
1840 else if ( flag == CV_TOO_MUCH_WORK || flag == CV_CONV_FAILURE || flag == CV_ERR_FAILURE)
1842 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1844 sciprint(_("****SUNDIALS.Cvode: too much work at time=%g (stiff region, change RTOL and ATOL)\n"), *told);
1850 *ierr = 300 + (-flag);
1859 *ierr = 300 + (-flag); /* raising errors due to internal errors, otherwise error due to flagr*/
1865 if (flag == CV_ZERO_DETACH_RETURN)
1868 }; /* new feature of sundials, detects zero-detaching */
1870 if (flag == CV_ROOT_RETURN)
1872 /* . at a least one root has been found */
1874 if (Discrete_Jump == 0)
1876 flagr = ODEGetRootInfo(ode_mem, jroot);
1877 if (check_flag(&flagr, "CVodeGetRootInfo", 1))
1879 *ierr = 300 + (-flagr);
1884 /* . at a least one root has been found */
1885 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
1887 sciprint(_("root found at t=: %f\n"), *told);
1889 /* . update outputs affecting ztyp blocks ONLY FOR OLD BLOCKS */
1890 zdoit(told, x, xd, g);
1896 for (jj = 0; jj < ng; ++jj)
1898 C2F(curblk).kfun = zcros[ jj];
1899 if (C2F(curblk).kfun == -1)
1905 for (j = zcptr[C2F(curblk).kfun] - 1 ;
1906 j < zcptr[C2F(curblk).kfun + 1] - 1 ; ++j)
1917 Blocks[C2F(curblk).kfun - 1].jroot = &jroot[zcptr[C2F(curblk).kfun] - 1];
1918 if (funtyp[C2F(curblk).kfun] > 0)
1921 if (Blocks[C2F(curblk).kfun - 1].nevout > 0)
1924 if (Blocks[C2F(curblk).kfun - 1].nx > 0)
1926 Blocks[C2F(curblk).kfun - 1].x = &x[xptr[C2F(curblk).kfun] - 1];
1927 Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
1929 /* call corresponding block to determine output event (kev) */
1930 Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
1931 /*callf(told, xd, x, x,g,&flag__);*/
1932 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1939 /* . update event agenda */
1940 for (k = 0; k < Blocks[C2F(curblk).kfun - 1].nevout; ++k)
1942 if (Blocks[C2F(curblk).kfun - 1].evout[k] >= 0.)
1944 i3 = k + clkptr[C2F(curblk).kfun] ;
1945 addevs(Blocks[C2F(curblk).kfun - 1].evout[k] + (*told), &i3, &ierr1);
1948 /* . nevts too small */
1956 /* . update state */
1957 if (Blocks[C2F(curblk).kfun - 1].nx > 0)
1959 /* . call corresponding block to update state */
1961 Blocks[C2F(curblk).kfun - 1].x = &x[xptr[C2F(curblk).kfun] - 1];
1962 Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
1963 Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
1964 /*callf(told, xd, x, x,g,&flag__);*/
1965 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
1979 /*--discrete zero crossings----dzero--------------------*/
1980 if (ng > 0) /* storing ZC signs just after a sundials call*/
1982 zdoit(told, x, x, g);
1988 for (jj = 0; jj < ng; ++jj)
2000 /*--discrete zero crossings----dzero--------------------*/
2002 C2F(realtime)(told);
2007 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2009 sciprint(_("Event: %d activated at t=%f\n"), *pointi, *told);
2010 for (kev = 0; kev < nblk; kev++)
2012 if (Blocks[kev].nmode > 0)
2014 sciprint(_("mode of block %d=%d, "), kev, Blocks[kev].mode[0]);
2017 sciprint(_("**mod**\n"));
2021 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
2023 sciprint(_("End of activation\n"));
2032 /* end of main loop on time */
2037 /*--------------------------------------------------------------------------*/
2038 static void cossimdaskr(double *told)
2040 /* System generated locals */
2042 //** used for the [stop] button
2043 static char CommandToUnstack[1024];
2044 static int CommandLength = 0;
2045 static int SeqSync = 0;
2048 /* Local variables */
2049 static scicos_flag flag__ = 0;
2050 static int ierr1 = 0;
2051 static int j = 0, k = 0;
2052 static double t = 0.;
2054 static double rhotmp = 0., tstop = 0.;
2055 static int inxsci = 0;
2056 static int kpo = 0, kev = 0;
2058 int *jroot = NULL, *zcros = NULL;
2059 int *Mode_save = NULL;
2060 int Mode_change = 0;
2062 int flag = 0, flagr = 0;
2063 N_Vector yy = NULL, yp = NULL;
2064 realtype reltol = 0., abstol = 0.;
2065 int Discrete_Jump = 0;
2066 N_Vector IDx = NULL;
2067 realtype *scicos_xproperty = NULL;
2068 DlsMat TJacque = NULL;
2070 void *dae_mem = NULL;
2071 UserData data = NULL;
2072 IDAMem copy_IDA_mem = NULL;
2073 int maxnj = 0, maxnit = 0;
2074 /*-------------------- Analytical Jacobian memory allocation ----------*/
2075 int Jn = 0, Jnx = 0, Jno = 0, Jni = 0, Jactaille = 0;
2077 int cnt = 0, N_iters = 0;
2079 // Defining function pointers, for more readability
2080 int (* DAESolve) (void*, realtype, realtype*, N_Vector, N_Vector, int);
2081 int (* DAEReInit) (void*, realtype, N_Vector, N_Vector);
2082 int (* DAESetMaxStep) (void*, realtype);
2083 int (* DAESetUserData) (void*, void*);
2084 int (* DAESetStopTime) (void*, realtype);
2085 int (* DAEGetRootInfo) (void*, int*);
2086 int (* DAESStolerances) (void*, realtype, realtype);
2087 /* Generic flag for stop mode */
2088 int DAE_NORMAL = 1; /* IDA_NORMAL = 1 */
2089 int DAE_ONE_STEP = 2; /* IDA_ONE_STEP = 2 */
2090 switch (C2F(cmsolver).solver)
2093 DAESolve = &IDASolve;
2094 DAEReInit = &IDAReInit;
2095 DAESetMaxStep = &IDASetMaxStep;
2096 DAESetUserData = &IDASetUserData;
2097 DAESetStopTime = &IDASetStopTime;
2098 DAEGetRootInfo = &IDAGetRootInfo;
2099 DAESStolerances = &IDASStolerances;
2101 default: // Unknown solver number
2107 /* Set extension of Sundials for scicos */
2108 set_sundials_with_extension(TRUE);
2110 // CI=1.0; /* for function Get_Jacobian_ci */
2114 if ((jroot = MALLOC(sizeof(int) * ng)) == NULL )
2120 for ( jj = 0 ; jj < ng ; jj++ )
2128 if ((zcros = MALLOC(sizeof(int) * ng)) == NULL )
2142 if ((Mode_save = MALLOC(sizeof(int) * nmod)) == NULL )
2157 reltol = (realtype) rtol;
2158 abstol = (realtype) Atol; /* Ith(abstol,1) = (realtype) Atol;*/
2163 yy = N_VNewEmpty_Serial(*neq);
2164 if (check_flag((void *)yy, "N_VNew_Serial", 0))
2182 yp = N_VNewEmpty_Serial(*neq);
2183 if (check_flag((void *)yp, "N_VNew_Serial", 0))
2187 N_VDestroy_Serial(yy);
2206 IDx = N_VNew_Serial(*neq);
2207 if (check_flag((void *)IDx, "N_VNew_Serial", 0))
2212 N_VDestroy_Serial(yp);
2216 N_VDestroy_Serial(yy);
2233 /* Call the Create and Init functions to initialize DAE memory */
2235 dae_mem = IDACreate();
2236 if (check_flag((void *)dae_mem, "IDACreate", 0))
2240 N_VDestroy_Serial(IDx);
2244 N_VDestroy_Serial(yp);
2248 N_VDestroy_Serial(yy);
2264 copy_IDA_mem = (IDAMem) dae_mem;
2266 if (C2F(cmsolver).solver == 100)
2268 flag = IDASetErrHandlerFn(dae_mem, SundialsErrHandler, NULL);
2270 if (check_flag(&flag, "IDASetErrHandlerFn", 1))
2272 *ierr = 200 + (-flag);
2273 if (*neq > 0)if (C2F(cmsolver).solver == 100)
2279 N_VDestroy_Serial(IDx);
2283 N_VDestroy_Serial(yp);
2287 N_VDestroy_Serial(yy);
2304 if (C2F(cmsolver).solver == 100)
2306 flag = IDAInit(dae_mem, simblkdaskr, T0, yy, yp);
2308 if (check_flag(&flag, "IDAInit", 1))
2310 *ierr = 200 + (-flag);
2311 if (*neq > 0)if (C2F(cmsolver).solver == 100)
2317 N_VDestroy_Serial(IDx);
2321 N_VDestroy_Serial(yp);
2325 N_VDestroy_Serial(yy);
2342 flag = DAESStolerances(dae_mem, reltol, abstol);
2343 if (check_flag(&flag, "IDASStolerances", 1))
2345 *ierr = 200 + (-flag);
2346 if (*neq > 0)if (C2F(cmsolver).solver == 100)
2352 N_VDestroy_Serial(IDx);
2356 N_VDestroy_Serial(yp);
2360 N_VDestroy_Serial(yy);
2377 if (C2F(cmsolver).solver == 100)
2379 flag = IDARootInit(dae_mem, ng, grblkdaskr);
2381 if (check_flag(&flag, "IDARootInit", 1))
2383 *ierr = 200 + (-flag);
2384 if (*neq > 0)if (C2F(cmsolver).solver == 100)
2390 N_VDestroy_Serial(IDx);
2394 N_VDestroy_Serial(yp);
2398 N_VDestroy_Serial(yy);
2415 if (C2F(cmsolver).solver == 100)
2417 flag = IDADense(dae_mem, *neq);
2419 if (check_flag(&flag, "IDADense", 1))
2421 *ierr = 200 + (-flag);
2422 if (*neq > 0)if (C2F(cmsolver).solver == 100)
2428 N_VDestroy_Serial(IDx);
2432 N_VDestroy_Serial(yp);
2436 N_VDestroy_Serial(yy);
2454 if ((data = (UserData) MALLOC(sizeof(*data))) == NULL)
2457 if (*neq > 0)if (C2F(cmsolver).solver == 100)
2463 N_VDestroy_Serial(IDx);
2467 N_VDestroy_Serial(yp);
2471 N_VDestroy_Serial(yy);
2487 data->dae_mem = dae_mem;
2493 data->ewt = N_VNew_Serial(*neq);
2494 if (check_flag((void *)data->ewt, "N_VNew_Serial", 0))
2496 *ierr = 200 + (-flag);
2501 if (*neq > 0)if (C2F(cmsolver).solver == 100)
2507 N_VDestroy_Serial(IDx);
2511 N_VDestroy_Serial(yp);
2515 N_VDestroy_Serial(yy);
2533 if ((data->gwork = (double *) MALLOC(ng * sizeof(double))) == NULL)
2537 N_VDestroy_Serial(data->ewt);
2543 if (*neq > 0)if (C2F(cmsolver).solver == 100)
2549 N_VDestroy_Serial(IDx);
2553 N_VDestroy_Serial(yp);
2557 N_VDestroy_Serial(yy);
2574 /*Jacobian_Flag=0; */
2575 if (AJacobian_block > 0) /* set by the block with A-Jac in flag-4 using Set_Jacobian_flag(1); */
2578 Jnx = Blocks[AJacobian_block - 1].nx;
2579 Jno = Blocks[AJacobian_block - 1].nout;
2580 Jni = Blocks[AJacobian_block - 1].nin;
2589 Jactaille = 3 * Jn + (Jn + Jni) * (Jn + Jno) + Jnx * (Jni + 2 * Jn + Jno) + (Jn - Jnx) * (2 * (Jn - Jnx) + Jno + Jni) + 2 * Jni * Jno;
2591 if ((data->rwork = (double *) MALLOC(Jactaille * sizeof(double))) == NULL)
2599 N_VDestroy_Serial(data->ewt);
2605 if (*neq > 0)if (C2F(cmsolver).solver == 100)
2611 N_VDestroy_Serial(IDx);
2615 N_VDestroy_Serial(yp);
2619 N_VDestroy_Serial(yy);
2637 if (C2F(cmsolver).solver == 100)
2639 flag = IDADlsSetDenseJacFn(dae_mem, Jacobians);
2641 if (check_flag(&flag, "IDADlsSetDenseJacFn", 1))
2643 *ierr = 200 + (-flag);
2648 TJacque = (DlsMat) NewDenseMat(*neq, *neq);
2650 flag = DAESetUserData(dae_mem, data);
2651 if (check_flag(&flag, "IDASetUserData", 1))
2653 *ierr = 200 + (-flag);
2660 flag = DAESetMaxStep(dae_mem, (realtype) hmax);
2661 if (check_flag(&flag, "IDASetMaxStep", 1))
2663 *ierr = 200 + (-flag);
2669 maxnj = 100; /* setting the maximum number of Jacobian evaluations during a Newton step */
2670 if (C2F(cmsolver).solver == 100)
2672 flag = IDASetMaxNumJacsIC(dae_mem, maxnj);
2674 if (check_flag(&flag, "IDASetMaxNumJacsIC", 1))
2676 *ierr = 200 + (-flag);
2681 maxnit = 10; /* setting the maximum number of Newton iterations in any attempt to solve CIC */
2682 if (C2F(cmsolver).solver == 100)
2684 flag = IDASetMaxNumItersIC(dae_mem, maxnit);
2686 if (check_flag(&flag, "IDASetMaxNumItersIC", 1))
2688 *ierr = 200 + (-flag);
2693 /* setting the maximum number of steps in an integration interval */
2694 if (C2F(cmsolver).solver == 100)
2696 flag = IDASetMaxNumSteps(dae_mem, 2000);
2698 if (check_flag(&flag, "IDASetMaxNumSteps", 1))
2700 *ierr = 200 + (-flag);
2705 } /* testing if neq>0 */
2710 uround = uround * 0.5;
2712 while ( 1.0 + uround != 1.0);
2713 uround = uround * 2.0;
2714 SQuround = sqrt(uround);
2717 C2F(coshlt).halt = 0;
2724 C2F(xscion)(&inxsci);
2725 /* initialization */
2726 C2F(realtimeinit)(told, &C2F(rtfactor).scale);
2727 /* ATOL and RTOL are scalars */
2730 for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
2732 if (Blocks[C2F(curblk).kfun - 1].ng >= 1)
2734 zcros[jj] = C2F(curblk).kfun;
2738 /* . ng >= jj required */
2743 /* initialization (propagation of constant blocks outputs) */
2751 /*--discrete zero crossings----dzero--------------------*/
2752 if (ng > 0) /* storing ZC signs just after a solver call*/
2754 zdoit(told, x, x, g);
2760 for (jj = 0; jj < ng; ++jj)
2770 /* main loop on time */
2773 while (ismenu()) //** if the user has done something, do the actions
2776 SeqSync = GetCommand(CommandToUnstack); //** get to the action
2777 CommandLength = (int)strlen(CommandToUnstack);
2778 syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
2780 if (C2F(coshlt).halt != 0)
2782 if (C2F(coshlt).halt == 2)
2784 *told = *tf; /* end simulation */
2786 C2F(coshlt).halt = 0;
2798 if (fabs(t - *told) < ttol)
2801 /* update output part */
2805 /* ! scheduling problem */
2812 if (xptr[nblk + 1] == 1)
2814 /* . no continuous state */
2815 if (*told + deltat + ttol > t)
2823 /* . update outputs of 'c' type blocks with no continuous state */
2826 /* . we are at the end, update continuous part before leaving */
2834 rhotmp = *tf + ttol;
2839 if (critev[kpo] == 1)
2841 rhotmp = tevts[kpo];
2852 hot = 0;/* Cold-restart the solver if the new TSTOP isn't beyong the previous one*/
2856 t = Min(*told + deltat, Min(t, *tf + ttol));
2858 if (hot == 0) /* CIC calculation when hot==0 */
2861 /* Setting the stop time*/
2862 flag = DAESetStopTime(dae_mem, (realtype)tstop);
2863 if (check_flag(&flag, "IDASetStopTime", 1))
2865 *ierr = 200 + (-flag);
2870 if (ng > 0 && nmod > 0)
2873 zdoit(told, x, xd, g);
2881 /*----------ID setting/checking------------*/
2882 N_VConst(ONE, IDx); /* Initialize id to 1's. */
2883 scicos_xproperty = NV_DATA_S(IDx);
2890 for (jj = 0; jj < *neq; jj++)
2894 scicos_xproperty[jj] = ONE;
2896 if (xprop[jj] == -1)
2898 scicos_xproperty[jj] = ZERO;
2901 /* CI=0.0;CJ=100.0; // for functions Get_Jacobian_ci and Get_Jacobian_cj
2902 Jacobians(*neq, (realtype) (*told), yy, yp, bidon, (realtype) CJ, data, TJacque, tempv1, tempv2, tempv3);
2903 for (jj=0;jj<*neq;jj++){
2904 Jacque_col=DENSE_COL(TJacque,jj);
2906 for (kk=0;kk<*neq;kk++){
2907 if ((Jacque_col[kk]-Jacque_col[kk]!=0)) {
2911 if (Jacque_col[kk]!=0){
2917 if (CI>=ZERO){ scicos_xproperty[jj]=CI;}else{fprintf(stderr,"\nWarinng! Xproperties are not match for i=%d!",jj);}
2919 /* printf("\n"); for(jj=0;jj<*neq;jj++) { printf("x%d=%g ",jj,scicos_xproperty[jj]); }*/
2921 if (C2F(cmsolver).solver == 100)
2923 flag = IDASetId(dae_mem, IDx);
2925 if (check_flag(&flag, "IDASetId", 1))
2927 *ierr = 200 + (-flag);
2931 // CI=1.0; // for function Get_Jacobian_ci
2932 /*--------------------------------------------*/
2933 // maxnj=100; /* setting the maximum number of Jacobian evaluation during a Newton step */
2934 // flag=IDASetMaxNumJacsIC(dae_mem, maxnj);
2935 // if (check_flag(&flag, "IDASetMaxNumItersIC", 1)) {
2936 // *ierr=200+(-flag);
2940 // flag=IDASetLineSearchOffIC(dae_mem,FALSE); /* (def=false) */
2941 // if (check_flag(&flag, "IDASetLineSearchOffIC", 1)) {
2942 // *ierr=200+(-flag);
2946 // flag=IDASetMaxNumItersIC(dae_mem, 10);/* (def=10) setting the maximum number of Newton iterations in any attempt to solve CIC */
2947 // if (check_flag(&flag, "IDASetMaxNumItersIC", 1)) {
2948 // *ierr=200+(-flag);
2953 N_iters = 4 + nmod * 4;
2954 for (j = 0; j <= N_iters; j++)
2956 /* counter to reevaluate the
2957 modes in mode->CIC->mode->CIC-> loop
2958 do it once in the absence of mode (nmod=0) */
2959 /* updating the modes through Flag==9, Phase==1 */
2961 /* Serge Steer 29/06/2009 */
2962 while (ismenu()) //** if the user has done something, do the actions
2965 SeqSync = GetCommand(CommandToUnstack); //** get to the action
2966 CommandLength = (int)strlen(CommandToUnstack);
2967 syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
2969 if (C2F(coshlt).halt != 0)
2971 if (C2F(coshlt).halt == 2)
2973 *told = *tf; /* end simulation */
2975 C2F(coshlt).halt = 0;
2981 flag = DAEReInit(dae_mem, (realtype)(*told), yy, yp);
2982 if (check_flag(&flag, "IDAReInit", 1))
2984 *ierr = 200 + (-flag);
2989 phase = 2; /* IDACalcIC: PHI-> yy0: if (ok) yy0_cic-> PHI*/
2990 copy_IDA_mem->ida_kk = 1;
2992 // the initial conditons y0 and yp0 do not satisfy the DAE
2993 if (C2F(cmsolver).solver == 100)
2995 flagr = IDACalcIC(dae_mem, IDA_YA_YDP_INIT, (realtype)(t));
2999 if (C2F(cmsolver).solver == 100)
3001 flag = IDAGetConsistentIC(dae_mem, yy, yp); /* PHI->YY */
3004 if (*ierr > 5) /* *ierr>5 => singularity in block */
3010 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3014 sciprint(_("**** SUNDIALS.IDA successfully initialized *****\n") );
3018 sciprint(_("**** SUNDIALS.IDA failed to initialize ->try again *****\n") );
3021 /*-------------------------------------*/
3022 /* saving the previous modes*/
3023 for (jj = 0; jj < nmod; ++jj)
3025 Mode_save[jj] = mod[jj];
3027 if (ng > 0 && nmod > 0)
3030 zdoit(told, x, xd, g);
3037 /*------------------------------------*/
3039 for (jj = 0; jj < nmod; ++jj)
3041 if (Mode_save[jj] != mod[jj])
3047 if (Mode_change == 0)
3051 break; /* if (flagr>=0) break; else{ *ierr=200+(-flagr); freeallx; return; }*/
3053 else if (j >= (int)( N_iters / 2))
3055 /* IDASetMaxNumStepsIC(mem,10); */ /* maxnh (def=5) */
3056 if (C2F(cmsolver).solver == 100)
3058 IDASetMaxNumJacsIC(dae_mem, 10); /* maxnj 100 (def=4)*/
3060 /* IDASetMaxNumItersIC(mem,100000); */ /* maxnit in IDANewtonIC (def=10) */
3061 if (C2F(cmsolver).solver == 100)
3063 IDASetLineSearchOffIC(dae_mem, TRUE); /* (def=false) */
3065 /* IDASetNonlinConvCoefIC(mem,1.01);*/ /* (def=0.01-0.33*/
3066 if (C2F(cmsolver).solver == 100)
3068 flag = IDASetMaxNumItersIC(dae_mem, 1000);
3070 if (check_flag(&flag, "IDASetMaxNumItersIC", 1))
3072 *ierr = 200 + (-flag);
3078 }/* mode-CIC counter*/
3079 if (Mode_change == 1)
3081 /* In this case, we try again by relaxing all modes and calling IDA_calc again
3084 copy_IDA_mem->ida_kk = 1;
3085 if (C2F(cmsolver).solver == 100)
3087 flagr = IDACalcIC(dae_mem, IDA_YA_YDP_INIT, (realtype)(t));
3090 if (C2F(cmsolver).solver == 100)
3092 flag = IDAGetConsistentIC(dae_mem, yy, yp); /* PHI->YY */
3094 if ((flagr < 0) || (*ierr > 5)) /* *ierr>5 => singularity in block */
3101 /*-----If flagr<0 the initialization solver has not converged-----*/
3109 } /* CIC calculation when hot==0 */
3111 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3113 sciprint(_("****daskr from: %f to %f hot= %d \n"), *told, t, hot);
3116 /*--discrete zero crossings----dzero--------------------*/
3117 /*--check for Dzeros after Mode settings or ddoit()----*/
3119 if (ng > 0 && hot == 0)
3121 zdoit(told, x, xd, g);
3127 for (jj = 0; jj < ng; ++jj)
3129 if ((g[jj] >= 0.0) && ( jroot[jj] == -5))
3134 else if ((g[jj] < 0.0) && ( jroot[jj] == 5))
3146 /*--discrete zero crossings----dzero--------------------*/
3147 if (Discrete_Jump == 0) /* if there was a dzero, its event should be activated*/
3150 flagr = DAESolve(dae_mem, t, told, yy, yp, DAE_NORMAL);
3160 flagr = IDA_ROOT_RETURN; /* in order to handle discrete jumps */
3164 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3166 sciprint(_("****SUNDIALS.Ida reached: %f\n"), *told);
3171 else if ( flagr == IDA_TOO_MUCH_WORK || flagr == IDA_CONV_FAIL || flagr == IDA_ERR_FAIL)
3173 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3175 sciprint(_("**** SUNDIALS.Ida: too much work at time=%g (stiff region, change RTOL and ATOL)\n"), *told);
3181 *ierr = 200 + (-flagr);
3190 *ierr = 200 + (-flagr); /* raising errors due to internal errors, otherwise error due to flagr*/
3196 /* update outputs of 'c' type blocks if we are at the end*/
3204 if (flagr == IDA_ZERO_DETACH_RETURN)
3207 }; /* new feature of sundials, detects unmasking */
3208 if (flagr == IDA_ROOT_RETURN)
3210 /* . at a least one root has been found */
3212 if (Discrete_Jump == 0)
3214 flagr = DAEGetRootInfo(dae_mem, jroot);
3215 if (check_flag(&flagr, "IDAGetRootInfo", 1))
3217 *ierr = 200 + (-flagr);
3223 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3225 sciprint(_("root found at t=: %f\n"), *told);
3227 /* . update outputs affecting ztyp blocks ONLY FOR OLD BLOCKS*/
3228 zdoit(told, x, xd, g);
3234 for (jj = 0; jj < ng; ++jj)
3236 C2F(curblk).kfun = zcros[jj];
3237 if (C2F(curblk).kfun == -1)
3242 for (j = zcptr[C2F(curblk).kfun] - 1 ;
3243 j < zcptr[C2F(curblk).kfun + 1] - 1 ; ++j)
3253 Blocks[C2F(curblk).kfun - 1].jroot = &jroot[zcptr[C2F(curblk).kfun] - 1];
3254 if (funtyp[C2F(curblk).kfun] > 0)
3256 if (Blocks[C2F(curblk).kfun - 1].nevout > 0)
3259 if (Blocks[C2F(curblk).kfun - 1].nx > 0)
3261 Blocks[C2F(curblk).kfun - 1].x = &x[xptr[C2F(curblk).kfun] - 1];
3262 Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
3264 /* call corresponding block to determine output event (kev) */
3265 Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
3266 /*callf(told, xd, x, x,g,&flag__);*/
3267 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
3274 /* update event agenda */
3275 for (k = 0; k < Blocks[C2F(curblk).kfun - 1].nevout; ++k)
3277 if (Blocks[C2F(curblk).kfun - 1].evout[k] >= 0)
3279 i3 = k + clkptr[C2F(curblk).kfun] ;
3280 addevs(Blocks[C2F(curblk).kfun - 1].evout[k] + (*told), &i3, &ierr1);
3283 /* . nevts too small */
3292 if ((Blocks[C2F(curblk).kfun - 1].nx > 0) || (*Blocks[C2F(curblk).kfun - 1].work != NULL) )
3294 /* call corresponding block to update state */
3296 if (Blocks[C2F(curblk).kfun - 1].nx > 0)
3298 Blocks[C2F(curblk).kfun - 1].x = &x[xptr[C2F(curblk).kfun] - 1];
3299 Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
3301 Blocks[C2F(curblk).kfun - 1].nevprt = -kev;
3303 Blocks[C2F(curblk).kfun - 1].xprop = &xprop[-1 + xptr[C2F(curblk).kfun]];
3304 /*callf(told, xd, x, x,g,&flag__);*/
3305 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
3313 for (j = 0; j < *neq; j++) /* Adjust xprop for IDx */
3317 scicos_xproperty[j] = ONE;
3321 scicos_xproperty[j] = ZERO;
3329 /* Serge Steer 29/06/2009 */
3330 while (ismenu()) //** if the user has done something, do the actions
3333 SeqSync = GetCommand(CommandToUnstack); //** get to the action
3334 CommandLength = (int)strlen(CommandToUnstack);
3335 syncexec(CommandToUnstack, &CommandLength, &ierr2, &one, CommandLength); //** execute it
3338 if (C2F(coshlt).halt != 0)
3340 if (C2F(coshlt).halt == 2)
3342 *told = *tf; /* end simulation */
3344 C2F(coshlt).halt = 0;
3361 /*--discrete zero crossings----dzero--------------------*/
3362 if (ng > 0) /* storing ZC signs just after a ddaskr call*/
3364 zdoit(told, x, xd, g);
3370 for (jj = 0; jj < ng; ++jj)
3382 /*--discrete zero crossings----dzero--------------------*/
3384 C2F(realtime)(told);
3389 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3391 sciprint(_("Event: %d activated at t=%f\n"), *pointi, *told);
3395 if ((C2F(cosdebug).cosd >= 1) && (C2F(cosdebug).cosd != 3))
3397 sciprint(_("End of activation"));
3405 /* end of main loop on time */
3408 } /* cossimdaskr_ */
3409 /*--------------------------------------------------------------------------*/
3410 /* Subroutine cosend */
3411 static void cosend(double *told)
3413 /* Local variables */
3414 static scicos_flag flag__ = 0;
3416 static int kfune = 0;
3420 /* loop on blocks */
3421 for (C2F(curblk).kfun = 1; C2F(curblk).kfun <= nblk; ++C2F(curblk).kfun)
3424 Blocks[C2F(curblk).kfun - 1].nevprt = 0;
3425 if (funtyp[C2F(curblk).kfun] >= 0)
3427 if (Blocks[C2F(curblk).kfun - 1].nx > 0)
3429 Blocks[C2F(curblk).kfun - 1].x = &x[xptr[C2F(curblk).kfun] - 1];
3430 Blocks[C2F(curblk).kfun - 1].xd = &xd[xptr[C2F(curblk).kfun] - 1];
3432 /*callf(told, xd, x, x,x,&flag__);*/
3433 callf(told, &Blocks[C2F(curblk).kfun - 1], &flag__);
3434 if (flag__ < 0 && *ierr == 0)
3437 kfune = C2F(curblk).kfun;
3443 C2F(curblk).kfun = kfune;
3447 /*--------------------------------------------------------------------------*/
3449 void callf(double *t, scicos_block *block, scicos_flag *flag)
3451 double* args[SZ_SIZE];
3453 double intabl[TB_SIZE];
3454 double outabl[TB_SIZE];
3456 int ii = 0, in = 0, out = 0, ki = 0, ko = 0, no = 0, ni = 0, k = 0, j = 0;
3457 int szi = 0, flagi = 0;
3458 double *ptr_d = NULL;
3460 /* function pointers type def */
3464 /* ScicosFm1 loc3;*/
3472 int solver = C2F(cmsolver).solver;
3473 int cosd = C2F(cosdebug).cosd;
3474 /*int kf = C2F(curblk).kfun;*/
3476 block_error = (int*) flag;
3478 /* debug block is never called */
3479 /*if (kf==(debug_block+1)) return;*/
3480 if (block->type == 99)
3485 /* flag 7 implicit initialization */
3486 flagi = (int) * flag;
3487 /* change flag to zero if flagi==7 for explicit block */
3488 if (flagi == 7 && block->type < 10000)
3493 /* display information for debugging mode */
3498 sciprint(_("block %d is called "), C2F(curblk).kfun);
3499 sciprint(_("with flag %d "), *flag);
3500 sciprint(_("at time %f \n"), *t);
3502 if (debug_block > -1)
3506 sciprint(_("Entering the block \n"));
3508 call_debug_scicos(block, flag, flagi, debug_block);
3511 return; /* error in debug block */
3516 C2F(scsptr).ptr = block->scsptr;
3518 /* get pointer of the function */
3521 /* continuous state */
3522 if ((solver == 100) && block->type < 10000 && *flag == 0)
3525 block->xd = block->res;
3529 //sciprint("callf type=%d flag=%d\n",block->type,flagi);
3530 switch (block->type)
3532 /*******************/
3533 /* function type 0 */
3534 /*******************/
3537 /* This is for compatibility */
3538 /* jroot is returned in g for old type */
3539 if (block->nevprt < 0)
3541 for (j = 0; j < block->ng; ++j)
3543 block->g[j] = (double)block->jroot[j];
3547 /* concatenated entries and concatened outputs */
3548 /* catenate inputs if necessary */
3553 for (in = 0; in < block->nin; in++)
3555 szi = block->insz[in] * block->insz[in + block->nin];
3556 for (ii = 0; ii < szi; ii++)
3558 intabl[ki++] = *((double *)(block->inptr[in]) + ii);
3562 args[0] = &(intabl[0]);
3566 if (block->nin == 0)
3572 args[0] = (double *)(block->inptr[0]);
3573 ni = block->insz[0] * block->insz[1];
3577 /* catenate outputs if necessary */
3579 if (block->nout > 1)
3582 for (out = 0; out < block->nout; out++)
3584 szi = block->outsz[out] * block->outsz[out + block->nout];
3585 for (ii = 0; ii < szi; ii++)
3587 outabl[ko++] = *((double *)(block->outptr[out]) + ii);
3591 args[1] = &(outabl[0]);
3595 if (block->nout == 0)
3601 args[1] = (double *)(block->outptr[0]);
3602 no = block->outsz[0] * block->outsz[1];
3606 loc0 = (ScicosF0) loc;
3608 (*loc0)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3609 block->z, &block->nz,
3610 block->evout, &block->nevout, block->rpar, &block->nrpar,
3611 block->ipar, &block->nipar, (double *)args[0], &ni,
3612 (double *)args[1], &no);
3614 /* split output vector on each port if necessary */
3615 if (block->nout > 1)
3618 for (out = 0; out < block->nout; out++)
3620 szi = block->outsz[out] * block->outsz[out + block->nout];
3621 for (ii = 0; ii < szi; ii++)
3623 *((double *)(block->outptr[out]) + ii) = outabl[ko++];
3628 /* adjust values of output register */
3629 for (in = 0; in < block->nevout; ++in)
3631 block->evout[in] = block->evout[in] - *t;
3637 /*******************/
3638 /* function type 1 */
3639 /*******************/
3642 /* This is for compatibility */
3643 /* jroot is returned in g for old type */
3644 if (block->nevprt < 0)
3646 for (j = 0; j < block->ng; ++j)
3648 block->g[j] = (double)block->jroot[j];
3652 /* one entry for each input or output */
3653 for (in = 0 ; in < block->nin ; in++)
3655 args[in] = block->inptr[in];
3656 sz[in] = block->insz[in];
3658 for (out = 0; out < block->nout; out++)
3660 args[in + out] = block->outptr[out];
3661 sz[in + out] = block->outsz[out];
3663 /* with zero crossing */
3664 if (block->ztyp > 0)
3666 args[block->nin + block->nout] = block->g;
3667 sz[block->nin + block->nout] = block->ng;
3670 loc1 = (ScicosF) loc;
3672 (*loc1)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3673 block->z, &block->nz,
3674 block->evout, &block->nevout, block->rpar, &block->nrpar,
3675 block->ipar, &block->nipar,
3676 (double *)args[0], &sz[0],
3677 (double *)args[1], &sz[1], (double *)args[2], &sz[2],
3678 (double *)args[3], &sz[3], (double *)args[4], &sz[4],
3679 (double *)args[5], &sz[5], (double *)args[6], &sz[6],
3680 (double *)args[7], &sz[7], (double *)args[8], &sz[8],
3681 (double *)args[9], &sz[9], (double *)args[10], &sz[10],
3682 (double *)args[11], &sz[11], (double *)args[12], &sz[12],
3683 (double *)args[13], &sz[13], (double *)args[14], &sz[14],
3684 (double *)args[15], &sz[15], (double *)args[16], &sz[16],
3685 (double *)args[17], &sz[17]);
3687 /* adjust values of output register */
3688 for (in = 0; in < block->nevout; ++in)
3690 block->evout[in] = block->evout[in] - *t;
3696 /*******************/
3697 /* function type 2 */
3698 /*******************/
3701 /* This is for compatibility */
3702 /* jroot is returned in g for old type */
3703 if (block->nevprt < 0)
3705 for (j = 0; j < block->ng; ++j)
3707 block->g[j] = (double)block->jroot[j];
3711 /* no zero crossing */
3712 if (block->ztyp == 0)
3714 loc2 = (ScicosF2) loc;
3715 (*loc2)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3716 block->z, &block->nz,
3717 block->evout, &block->nevout, block->rpar, &block->nrpar,
3718 block->ipar, &block->nipar, (double **)block->inptr,
3719 block->insz, &block->nin,
3720 (double **)block->outptr, block->outsz, &block->nout);
3722 /* with zero crossing */
3725 loc2z = (ScicosF2z) loc;
3726 (*loc2z)(flag, &block->nevprt, t, block->xd, block->x, &block->nx,
3727 block->z, &block->nz,
3728 block->evout, &block->nevout, block->rpar, &block->nrpar,
3729 block->ipar, &block->nipar, (double **)block->inptr,
3730 block->insz, &block->nin,
3731 (double **)block->outptr, block->outsz, &block->nout,
3732 block->g, &block->ng);
3735 /* adjust values of output register */
3736 for (in = 0; in < block->nevout; ++in)
3738 block->evout[in] = block->evout[in] - *t;
3744 /*******************/
3745 /* function type 4 */
3746 /*******************/
3749 /* get pointer of the function type 4*/
3750 loc4 = (ScicosF4) loc;
3752 (*loc4)(block, *flag);
3757 /***********************/
3758 /* function type 10001 */
3759 /***********************/
3762 /* This is for compatibility */
3763 /* jroot is returned in g for old type */
3764 if (block->nevprt < 0)
3766 for (j = 0; j < block->ng; ++j)
3768 block->g[j] = (double)block->jroot[j];
3772 /* implicit block one entry for each input or output */
3773 for (in = 0 ; in < block->nin ; in++)
3775 args[in] = block->inptr[in];
3776 sz[in] = block->insz[in];
3778 for (out = 0; out < block->nout; out++)
3780 args[in + out] = block->outptr[out];
3781 sz[in + out] = block->outsz[out];
3783 /* with zero crossing */
3784 if (block->ztyp > 0)
3786 args[block->nin + block->nout] = block->g;
3787 sz[block->nin + block->nout] = block->ng;
3790 loci1 = (ScicosFi) loc;
3791 (*loci1)(flag, &block->nevprt, t, block->res, block->xd, block->x,
3792 &block->nx, block->z, &block->nz,
3793 block->evout, &block->nevout, block->rpar, &block->nrpar,
3794 block->ipar, &block->nipar,
3795 (double *)args[0], &sz[0],
3796 (double *)args[1], &sz[1], (double *)args[2], &sz[2],
3797 (double *)args[3], &sz[3], (double *)args[4], &sz[4],
3798 (double *)args[5], &sz[5], (double *)args[6], &sz[6],
3799 (double *)args[7], &sz[7], (double *)args[8], &sz[8],
3800 (double *)args[9], &sz[9], (double *)args[10], &sz[10],
3801 (double *)args[11], &sz[11], (double *)args[12], &sz[12],
3802 (double *)args[13], &sz[13], (double *)args[14], &sz[14],
3803 (double *)args[15], &sz[15], (double *)args[16], &sz[16],
3804 (double *)args[17], &sz[17]);
3806 /* adjust values of output register */
3807 for (in = 0; in < block->nevout; ++in)
3809 block->evout[in] = block->evout[in] - *t;
3815 /***********************/
3816 /* function type 10002 */
3817 /***********************/
3820 /* This is for compatibility */
3821 /* jroot is returned in g for old type */
3822 if (block->nevprt < 0)
3824 for (j = 0; j < block->ng; ++j)
3826 block->g[j] = (double)block->jroot[j];
3830 /* implicit block, inputs and outputs given by a table of pointers */
3831 /* no zero crossing */
3832 if (block->ztyp == 0)
3834 loci2 = (ScicosFi2) loc;
3835 (*loci2)(flag, &block->nevprt, t, block->res,
3836 block->xd, block->x, &block->nx,
3837 block->z, &block->nz,
3838 block->evout, &block->nevout, block->rpar, &block->nrpar,
3839 block->ipar, &block->nipar, (double **)block->inptr,
3840 block->insz, &block->nin,
3841 (double **)block->outptr, block->outsz, &block->nout);
3843 /* with zero crossing */
3846 loci2z = (ScicosFi2z) loc;
3847 (*loci2z)(flag, &block->nevprt, t, block->res,
3848 block->xd, block->x, &block->nx,
3849 block->z, &block->nz,
3850 block->evout, &block->nevout, block->rpar, &block->nrpar,
3851 block->ipar, &block->nipar,
3852 (double **)block->inptr, block->insz, &block->nin,
3853 (double **)block->outptr, block->outsz, &block->nout,
3854 block->g, &block->ng);
3857 /* adjust values of output register */
3858 for (in = 0; in < block->nevout; ++in)
3860 block->evout[in] = block->evout[in] - *t;
3866 /***********************/
3867 /* function type 10004 */
3868 /***********************/
3871 /* get pointer of the function type 4*/
3872 loc4 = (ScicosF4) loc;
3874 (*loc4)(block, *flag);
3884 sciprint(_("Undefined Function type\n"));
3889 // sciprint("callf end flag=%d\n",*flag);
3890 /* Implicit Solver & explicit block & flag==0 */
3891 /* adjust continuous state vector after call */
3892 if ((solver == 100) && block->type < 10000 && *flag == 0)
3897 for (k = 0; k < block->nx; k++)
3899 block->res[k] = block->res[k] - block->xd[k];
3904 for (k = 0; k < block->nx; k++)
3906 block->xd[k] = block->res[k];
3914 if (debug_block > -1)
3918 return; /* error in block */
3922 sciprint(_("Leaving block %d \n"), C2F(curblk).kfun);
3924 call_debug_scicos(block, flag, flagi, debug_block);
3925 /*call_debug_scicos(flag,kf,flagi,debug_block);*/
3929 /*--------------------------------------------------------------------------*/
3930 /* call_debug_scicos */
3931 static void call_debug_scicos(scicos_block *block, scicos_flag *flag, int flagi, int deb_blk)
3934 int solver = C2F(cmsolver).solver, k = 0;
3936 double *ptr_d = NULL;
3938 C2F(cosdebugcounter).counter = C2F(cosdebugcounter).counter + 1;
3939 C2F(scsptr).ptr = Blocks[deb_blk].scsptr;
3941 loc = Blocks[deb_blk].funpt; /* GLOBAL */
3942 loc4 = (ScicosF4) loc;
3944 /* continuous state */
3945 if ((solver == 100) && block->type < 10000 && *flag == 0)
3948 block->xd = block->res;
3951 (*loc4)(block, *flag);
3953 /* Implicit Solver & explicit block & flag==0 */
3954 /* adjust continuous state vector after call */
3955 if ((solver == 100) && block->type < 10000 && *flag == 0)
3960 for (k = 0; k < block->nx; k++)
3962 block->res[k] = block->res[k] - block->xd[k];
3967 for (k = 0; k < block->nx; k++)
3969 block->xd[k] = block->res[k];
3976 sciprint(_("Error in the Debug block \n"));
3978 } /* call_debug_scicos */
3979 /*--------------------------------------------------------------------------*/
3981 static int simblk(realtype t, N_Vector yy, N_Vector yp, void *f_data)
3983 double tx = 0., *x = NULL, *xd = NULL;
3984 int i = 0, nantest = 0;
3987 x = (double *) NV_DATA_S(yy);
3988 xd = (double *) NV_DATA_S(yp);
3990 for (i = 0; i < *neq; i++)
3992 xd[i] = 0; /* Ã la place de "C2F(dset)(neq, &c_b14,xcdot , &c__1);"*/
3994 C2F(ierode).iero = 0;
3996 odoit(&tx, x, xd, xd);
3997 C2F(ierode).iero = *ierr;
4002 for (i = 0; i < *neq; i++) /* NaN checking */
4004 if ((xd[i] - xd[i] != 0))
4006 sciprint(_("\nWarning: The computing function #%d returns a NaN/Inf"), i);
4013 return 349; /* recoverable error; */
4017 return (abs(*ierr)); /* ierr>0 recoverable error; ierr>0 unrecoverable error; ierr=0: ok*/
4020 /*--------------------------------------------------------------------------*/
4022 static int grblk(realtype t, N_Vector yy, realtype *gout, void *g_data)
4024 double tx = 0., *x = NULL;
4025 int jj = 0, nantest = 0;
4028 x = (double *) NV_DATA_S(yy);
4030 C2F(ierode).iero = 0;
4033 zdoit(&tx, x, x, (double*) gout);
4038 for (jj = 0; jj < ng; jj++)
4039 if (gout[jj] - gout[jj] != 0)
4041 sciprint(_("\nWarning: The zero_crossing function #%d returns a NaN/Inf"), jj);
4044 } /* NaN checking */
4047 return 350; /* recoverable error; */
4050 C2F(ierode).iero = *ierr;
4054 /*--------------------------------------------------------------------------*/
4056 static int simblklsodar(int * nequations, realtype * tOld, realtype * actual, realtype * res)
4059 int i = 0, nantest = 0;
4061 tx = (double) * tOld;
4063 for (i = 0; i < *nequations; ++i)
4065 res[i] = 0; /* Ã la place de "C2F(dset)(neq, &c_b14,xcdot , &c__1);"*/
4067 C2F(ierode).iero = 0;
4069 odoit(&tx, actual, res, res);
4070 C2F(ierode).iero = *ierr;
4075 for (i = 0; i < *nequations; i++) /* NaN checking */
4077 if ((res[i] - res[i] != 0))
4079 sciprint(_("\nWarning: The computing function #%d returns a NaN/Inf"), i);
4086 return 349; /* recoverable error; */
4090 return (abs(*ierr)); /* ierr>0 recoverable error; ierr>0 unrecoverable error; ierr=0: ok*/
4092 } /* simblklsodar */
4093 /*--------------------------------------------------------------------------*/
4095 static int grblklsodar(int * nequations, realtype * tOld, realtype * actual, int * ngc, realtype * res)
4098 int jj = 0, nantest = 0;
4100 tx = (double) * tOld;
4102 C2F(ierode).iero = 0;
4105 zdoit(&tx, actual, actual, res);
4110 for (jj = 0; jj < *ngc; jj++)
4111 if (res[jj] - res[jj] != 0)
4113 sciprint(_("\nWarning: The zero_crossing function #%d returns a NaN/Inf"), jj);
4116 } /* NaN checking */
4119 return 350; /* recoverable error; */
4122 C2F(ierode).iero = *ierr;
4126 /*--------------------------------------------------------------------------*/
4128 static int simblkdaskr(realtype tres, N_Vector yy, N_Vector yp, N_Vector resval, void *rdata)
4131 double *xc = NULL, *xcdot = NULL, *residual = NULL;
4132 realtype alpha = 0.;
4138 int jj = 0, flag = 0, nantest = 0;
4140 data = (UserData) rdata;
4142 if (get_phase_simulation() == 1)
4144 /* Just to update mode in a very special case, i.e., when initialization using modes fails.
4145 in this case, we relax all modes and try again one more time.
4147 zdoit(&tx, NV_DATA_S(yy), NV_DATA_S(yp), (double *)data->gwork);
4151 flag = IDAGetCurrentStep(data->dae_mem, &hh);
4154 *ierr = 200 + (-flag);
4159 flag = IDAGetCurrentOrder(data->dae_mem, &qlast);
4162 *ierr = 200 + (-flag);
4167 for (jj = 0; jj < qlast; jj++)
4169 alpha = alpha - ONE / (jj + 1);
4172 // CJ=-alpha/hh; // For function Get_Jacobian_cj
4181 xc = (double *) NV_DATA_S(yy);
4182 xcdot = (double *) NV_DATA_S(yp);
4183 residual = (double *) NV_DATA_S(resval);
4186 C2F(dcopy)(neq, xcdot, &c__1, residual, &c__1);
4188 C2F(ierode).iero = 0;
4189 odoit(&tx, xc, xcdot, residual);
4191 C2F(ierode).iero = *ierr;
4196 for (jj = 0; jj < *neq; jj++)
4197 if (residual[jj] - residual[jj] != 0) /* NaN checking */
4199 //sciprint(_("\nWarning: The residual function #%d returns a NaN"),jj);
4205 return 257; /* recoverable error; */
4209 return (abs(*ierr)); /* ierr>0 recoverable error; ierr>0 unrecoverable error; ierr=0: ok*/
4211 /*--------------------------------------------------------------------------*/
4213 static int grblkdaskr(realtype t, N_Vector yy, N_Vector yp, realtype *gout, void *g_data)
4216 int jj = 0, nantest = 0;
4221 C2F(ierode).iero = 0;
4222 zdoit(&tx, NV_DATA_S(yy), NV_DATA_S(yp), (double *)gout);
4225 nantest = 0; /* NaN checking */
4226 for (jj = 0; jj < ng; jj++)
4228 if (gout[jj] - gout[jj] != 0)
4230 sciprint(_("\nWarning: The zero-crossing function #%d returns a NaN"), jj);
4237 return 258; /* recoverable error; */
4240 C2F(ierode).iero = *ierr;
4243 /*--------------------------------------------------------------------------*/
4244 /* Subroutine addevs */
4245 static void addevs(double t, int *evtnb, int *ierr1)
4247 static int i = 0, j = 0;
4251 if (evtspt[*evtnb] != -1)
4253 if ((evtspt[*evtnb] == 0) && (*pointi == *evtnb))
4260 if (*pointi == *evtnb)
4262 *pointi = evtspt[*evtnb]; /* remove from chain */
4267 while (*evtnb != evtspt[i])
4271 evtspt[i] = evtspt[*evtnb]; /* remove old evtnb from chain */
4272 if (TCritWarning == 0)
4274 sciprint(_("\n Warning: an event is reprogrammed at t=%g by removing another"), t );
4275 sciprint(_("\n (already programmed) event. There may be an error in"));
4276 sciprint(_("\n your model. Please check your model\n"));
4279 do_cold_restart(); /* the erased event could be a critical
4280 event, so do_cold_restart is added to
4281 refresh the critical event table */
4297 if (t < tevts[*pointi])
4299 evtspt[*evtnb] = *pointi;
4311 if (t >= tevts[evtspt[i]])
4324 evtspt[*evtnb] = evtspt[i];
4328 /*--------------------------------------------------------------------------*/
4329 /* Subroutine putevs */
4330 void putevs(double *t, int *evtnb, int *ierr1)
4334 if (evtspt[*evtnb] != -1)
4349 evtspt[*evtnb] = *pointi;
4352 /*--------------------------------------------------------------------------*/
4353 /* Subroutine idoit */
4354 static void idoit(double *told)
4356 /* initialisation (propagation of constant blocks outputs) */
4357 /* Copyright INRIA */
4360 scicos_flag flag = 0;
4364 ScicosImport *scs_imp = NULL;
4367 scs_imp = getscicosimportptr();
4370 for (j = 0; j < * (scs_imp->niord); j++)
4372 kf = &scs_imp->iord[j];
4373 C2F(curblk).kfun = *kf; /* */
4374 if (scs_imp->funtyp[*kf - 1] > -1)
4376 /* continuous state */
4377 if (scs_imp->blocks[*kf - 1].nx > 0)
4379 scs_imp->blocks[*kf - 1].x = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4380 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4382 scs_imp->blocks[*kf - 1].nevprt = scs_imp->iord[j + * (scs_imp->niord)];
4383 /*callf(told, xd, x, x,x,&flag);*/
4384 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4391 if (scs_imp->blocks[*kf - 1].nevout > 0)
4393 if (scs_imp->funtyp[*kf - 1] < 0)
4395 i = synchro_nev(scs_imp, *kf, ierr);
4400 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4401 putevs(told, &i2, &ierr1);
4404 /* event conflict */
4417 /*--------------------------------------------------------------------------*/
4418 /* Subroutine doit */
4419 static void doit(double *told)
4421 /* propagation of blocks outputs on discrete activations */
4422 /* Copyright INRIA */
4425 scicos_flag flag = 0;
4428 int ii = 0, kever = 0;
4430 ScicosImport *scs_imp = NULL;
4433 scs_imp = getscicosimportptr();
4436 *pointi = evtspt[kever];
4439 nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever - 1];
4445 for (ii = scs_imp->ordptr[kever - 1]; ii <= scs_imp->ordptr[kever] - 1 ; ii++)
4447 kf = &scs_imp->ordclk[ii - 1];
4448 C2F(curblk).kfun = *kf;
4449 if (scs_imp->funtyp[*kf - 1] > -1)
4451 /* continuous state */
4452 if (scs_imp->blocks[*kf - 1].nx > 0)
4454 scs_imp->blocks[*kf - 1].x = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4455 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4457 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
4459 /*callf(told, xd, x, x,x,&flag);*/
4460 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4468 /* Initialize tvec */
4469 if (scs_imp->blocks[*kf - 1].nevout > 0)
4471 if (scs_imp->funtyp[*kf - 1] < 0)
4473 i = synchro_nev(scs_imp, *kf, ierr);
4478 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4479 putevs(told, &i2, &ierr1);
4482 /* event conflict */
4495 /*--------------------------------------------------------------------------*/
4496 /* Subroutine cdoit */
4497 static void cdoit(double *told)
4499 /* propagation of continuous blocks outputs */
4500 /* Copyright INRIA */
4503 scicos_flag flag = 0;
4507 ScicosImport *scs_imp = NULL;
4510 scs_imp = getscicosimportptr();
4513 for (j = 0; j < * (scs_imp->ncord); j++)
4515 kf = &scs_imp->cord[j];
4516 C2F(curblk).kfun = *kf;
4517 /* continuous state */
4518 if (scs_imp->blocks[*kf - 1].nx > 0)
4520 scs_imp->blocks[*kf - 1].x = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4521 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4523 scs_imp->blocks[*kf - 1].nevprt = scs_imp->cord[j + * (scs_imp->ncord)];
4524 if (scs_imp->funtyp[*kf - 1] > -1)
4527 /*callf(told, xd, x, x,x,&flag);*/
4528 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4536 /* Initialize tvec */
4537 if (scs_imp->blocks[*kf - 1].nevout > 0)
4539 if (scs_imp->funtyp[*kf - 1] < 0)
4541 i = synchro_nev(scs_imp, *kf, ierr);
4546 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4547 putevs(told, &i2, &ierr1);
4550 /* event conflict */
4563 /*--------------------------------------------------------------------------*/
4564 /* Subroutine ddoit */
4565 static void ddoit(double *told)
4567 /* update states & event out on discrete activations */
4568 /* Copyright INRIA */
4571 scicos_flag flag = 0;
4573 int i = 0, i3 = 0, ierr1 = 0;
4574 int ii = 0, keve = 0;
4576 ScicosImport *scs_imp = NULL;
4579 scs_imp = getscicosimportptr();
4589 /* update continuous and discrete states on event */
4594 for (i = 0; i < kiwa; i++)
4597 if (critev[keve] != 0)
4601 i2 = scs_imp->ordptr[keve] - 1;
4602 for (ii = scs_imp->ordptr[keve - 1]; ii <= i2; ii++)
4604 kf = &scs_imp->ordclk[ii - 1];
4605 C2F(curblk).kfun = *kf;
4606 /* continuous state */
4607 if (scs_imp->blocks[*kf - 1].nx > 0)
4609 scs_imp->blocks[*kf - 1].x = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4610 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4613 scs_imp->blocks[*kf - 1].nevprt = scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1];
4615 if (scs_imp->blocks[*kf - 1].nevout > 0)
4617 if (scs_imp->funtyp[*kf - 1] >= 0)
4619 /* initialize evout */
4620 for (j = 0; j < scs_imp->blocks[*kf - 1].nevout; j++)
4622 scs_imp->blocks[*kf - 1].evout[j] = -1;
4626 if (scs_imp->blocks[*kf - 1].nevprt > 0) /* if event has continuous origin don't call*/
4628 /*callf(told, xd, x, x ,x,&flag);*/
4629 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4637 for (j = 0; j < scs_imp->blocks[*kf - 1].nevout; j++)
4639 if (scs_imp->blocks[*kf - 1].evout[j] >= 0.)
4641 i3 = j + scs_imp->clkptr[*kf - 1] ;
4642 addevs(scs_imp->blocks[*kf - 1].evout[j] + (*told), &i3, &ierr1);
4645 /* event conflict */
4654 if (scs_imp->blocks[*kf - 1].nevprt > 0)
4656 if (scs_imp->blocks[*kf - 1].nx + scs_imp->blocks[*kf - 1].nz + scs_imp->blocks[*kf - 1].noz > 0 || \
4657 *scs_imp->blocks[*kf - 1].work != NULL)
4659 /* if a hidden state exists, must also call (for new scope eg) */
4660 /* to avoid calling non-real activations */
4662 /*callf(told, xd, x, x,x,&flag);*/
4663 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4673 if (*scs_imp->blocks[*kf - 1].work != NULL)
4676 scs_imp->blocks[*kf - 1].nevprt = 0; /* in case some hidden continuous blocks need updating */
4677 /*callf(told, xd, x, x,x,&flag);*/
4678 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4689 /*--------------------------------------------------------------------------*/
4690 /* Subroutine edoit */
4691 static void edoit(double *told, int *kiwa)
4693 /* update blocks output on discrete activations */
4694 /* Copyright INRIA */
4697 scicos_flag flag = 0;
4698 int ierr1 = 0, i = 0;
4699 int kever = 0, ii = 0;
4701 ScicosImport *scs_imp = NULL;
4705 scs_imp = getscicosimportptr();
4710 *pointi = evtspt[kever];
4713 nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever - 1];
4720 for (ii = scs_imp->ordptr[kever - 1]; ii <= scs_imp->ordptr[kever] - 1; ii++)
4722 kf = &scs_imp->ordclk[ii - 1];
4723 C2F(curblk).kfun = *kf;
4725 if (scs_imp->funtyp[*kf - 1] > -1)
4727 /* continuous state */
4728 if (scs_imp->blocks[*kf - 1].nx > 0)
4730 scs_imp->blocks[*kf - 1].x = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4731 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4734 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
4737 /*callf(told, xd, x, x,x,&flag);*/
4738 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4746 /* Initialize tvec */
4747 if (scs_imp->blocks[*kf - 1].nevout > 0)
4749 if (scs_imp->funtyp[*kf - 1] < 0)
4751 i = synchro_nev(scs_imp, *kf, ierr);
4756 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4757 putevs(told, &i2, &ierr1);
4760 /* event conflict */
4773 /*--------------------------------------------------------------------------*/
4774 /* Subroutine odoit */
4775 static void odoit(double *told, double *xt, double *xtd, double *residual)
4777 /* update blocks derivative of continuous time block */
4778 /* Copyright INRIA */
4781 scicos_flag flag = 0;
4782 int keve = 0, kiwa = 0;
4783 int ierr1 = 0, i = 0;
4786 ScicosImport *scs_imp = NULL;
4789 scs_imp = getscicosimportptr();
4793 for (jj = 0; jj < * (scs_imp->noord); jj++)
4795 kf = &scs_imp->oord[jj];
4796 C2F(curblk).kfun = *kf;
4797 /* continuous state */
4798 if (scs_imp->blocks[*kf - 1].nx > 0)
4800 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
4801 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
4802 scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
4805 scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[jj + * (scs_imp->noord)];
4806 if (scs_imp->funtyp[*kf - 1] > -1)
4809 /*callf(told, xtd, xt, residual,x,&flag);*/
4810 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4818 if (scs_imp->blocks[*kf - 1].nevout > 0)
4820 if (scs_imp->funtyp[*kf - 1] < 0)
4822 if (scs_imp->blocks[*kf - 1].nmode > 0)
4824 i2 = scs_imp->blocks[*kf - 1].mode[0] + scs_imp->clkptr[*kf - 1] - 1;
4828 i = synchro_nev(scs_imp, *kf, ierr);
4833 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4835 putevs(told, &i2, &ierr1);
4838 /* event conflict */
4842 ozdoit(told, xt, xtd, &kiwa);
4851 /* update states derivatives */
4852 for (ii = 0; ii < * (scs_imp->noord); ii++)
4854 kf = &scs_imp->oord[ii];
4855 C2F(curblk).kfun = *kf;
4856 if (scs_imp->blocks[*kf - 1].nx > 0 || \
4857 *scs_imp->blocks[*kf - 1].work != NULL)
4859 /* work tests if a hidden state exists, used for delay block */
4861 /* continuous state */
4862 if (scs_imp->blocks[*kf - 1].nx > 0)
4864 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
4865 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
4866 scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
4868 scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[ii + * (scs_imp->noord)];
4869 /*callf(told, xtd, xt, residual,xt,&flag);*/
4870 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4880 for (i = 0; i < kiwa; i++)
4883 for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1; ii++)
4885 kf = &scs_imp->ordclk[ii - 1];
4886 C2F(curblk).kfun = *kf;
4887 if (scs_imp->blocks[*kf - 1].nx > 0 || \
4888 *scs_imp->blocks[*kf - 1].work != NULL)
4890 /* work tests if a hidden state exists */
4892 /* continuous state */
4893 if (scs_imp->blocks[*kf - 1].nx > 0)
4895 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
4896 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
4897 scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
4899 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
4900 /*callf(told, xtd, xt, residual,xt,&flag);*/
4901 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4912 /*--------------------------------------------------------------------------*/
4913 /* Subroutine reinitdoit */
4914 static void reinitdoit(double *told)
4916 /* update blocks xproperties of continuous time block */
4917 /* Copyright INRIA */
4920 scicos_flag flag = 0;
4921 int keve = 0, kiwa = 0;
4922 int ierr1 = 0, i = 0;
4925 ScicosImport *scs_imp = NULL;
4928 scs_imp = getscicosimportptr();
4932 for (jj = 0; jj < * (scs_imp->noord); jj++)
4934 kf = &scs_imp->oord[jj];
4935 C2F(curblk).kfun = *kf;
4936 /* continuous state */
4937 if (scs_imp->blocks[*kf - 1].nx > 0)
4939 scs_imp->blocks[*kf - 1].x = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4940 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4942 scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[jj + * (scs_imp->noord)];
4943 if (scs_imp->funtyp[*kf - 1] > -1)
4946 /*callf(told, xd, x, x,x,&flag);*/
4947 callf(told, &scs_imp->blocks[*kf - 1], &flag);
4955 if (scs_imp->blocks[*kf - 1].nevout > 0 && scs_imp->funtyp[*kf - 1] < 0)
4957 i = synchro_nev(scs_imp, *kf, ierr);
4962 if (scs_imp->blocks[*kf - 1].nmode > 0)
4964 scs_imp->blocks[*kf - 1].mode[0] = i;
4966 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
4967 putevs(told, &i2, &ierr1);
4970 /* event conflict */
4983 for (ii = 0; ii < * (scs_imp->noord); ii++)
4985 kf = &scs_imp->oord[ii];
4986 C2F(curblk).kfun = *kf;
4987 if (scs_imp->blocks[*kf - 1].nx > 0)
4990 scs_imp->blocks[*kf - 1].x = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
4991 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4992 scs_imp->blocks[*kf - 1].res = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
4993 scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[ii + * (scs_imp->noord)];
4994 scs_imp->blocks[*kf - 1].xprop = &scs_imp->xprop[-1 + scs_imp->xptr[*kf - 1]];
4995 /*callf(told, xd, x, xd,x,&flag);*/
4996 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5005 for (i = 0; i < kiwa; i++)
5008 for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1; ii++)
5010 kf = &scs_imp->ordclk[ii - 1];
5011 C2F(curblk).kfun = *kf;
5012 if (scs_imp->blocks[*kf - 1].nx > 0)
5015 scs_imp->blocks[*kf - 1].x = &scs_imp->x[scs_imp->xptr[*kf - 1] - 1];
5016 scs_imp->blocks[*kf - 1].xd = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5017 scs_imp->blocks[*kf - 1].res = &scs_imp->xd[scs_imp->xptr[*kf - 1] - 1];
5018 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5019 scs_imp->blocks[*kf - 1].xprop = &scs_imp->xprop[-1 + scs_imp->xptr[*kf - 1]];
5020 /*callf(told, xd, x, xd,x,&flag);*/
5021 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5031 /*--------------------------------------------------------------------------*/
5032 /* Subroutine ozdoit */
5033 static void ozdoit(double *told, double *xt, double *xtd, int *kiwa)
5035 /* update blocks output of continuous time block on discrete activations */
5036 /* Copyright INRIA */
5039 scicos_flag flag = 0;
5041 int ierr1 = 0, i = 0;
5042 int ii = 0, kever = 0;
5044 ScicosImport *scs_imp = NULL;
5047 scs_imp = getscicosimportptr();
5051 *pointi = evtspt[kever];
5054 nord = scs_imp->ordptr[kever] - scs_imp->ordptr[kever - 1];
5062 for (ii = scs_imp->ordptr[kever - 1]; ii <= scs_imp->ordptr[kever] - 1; ii++)
5064 kf = &scs_imp->ordclk[ii - 1];
5065 C2F(curblk).kfun = *kf;
5066 if (scs_imp->funtyp[*kf - 1] > -1)
5068 /* continuous state */
5069 if (scs_imp->blocks[*kf - 1].nx > 0)
5071 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
5072 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5074 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5076 /*callf(told, xtd, xt, xt,x,&flag);*/
5077 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5085 /* Initialize tvec */
5086 if (scs_imp->blocks[*kf - 1].nevout > 0)
5088 if (scs_imp->funtyp[*kf - 1] < 0)
5090 if (phase == 1 || scs_imp->blocks[*kf - 1].nmode == 0)
5092 i = synchro_nev(scs_imp, *kf, ierr);
5100 i = scs_imp->blocks[*kf - 1].mode[0];
5102 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5103 putevs(told, &i2, &ierr1);
5106 /* event conflict */
5110 ozdoit(told, xt, xtd, kiwa);
5115 /*--------------------------------------------------------------------------*/
5116 /* Subroutine zdoit */
5117 static void zdoit(double *told, double *xt, double *xtd, double *g)
5119 /* update blocks zcross of continuous time block */
5120 /* Copyright INRIA */
5122 scicos_flag flag = 0;
5123 int keve = 0, kiwa = 0;
5124 int ierr1 = 0, i = 0, j = 0;
5127 ScicosImport *scs_imp = NULL;
5130 scs_imp = getscicosimportptr();
5133 for (i = 0; i < * (scs_imp->ng); i++)
5139 for (jj = 0; jj < * (scs_imp->nzord); jj++)
5141 kf = &scs_imp->zord[jj];
5142 C2F(curblk).kfun = *kf;
5143 /* continuous state */
5144 if (scs_imp->blocks[*kf - 1].nx > 0)
5146 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
5147 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5149 scs_imp->blocks[*kf - 1].nevprt = scs_imp->zord[jj + * (scs_imp->nzord)];
5151 if (scs_imp->funtyp[*kf - 1] > -1)
5154 /*callf(told, xtd, xt, xt,xt,&flag);*/
5155 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5163 /* Initialize tvec */
5164 if (scs_imp->blocks[*kf - 1].nevout > 0)
5166 if (scs_imp->funtyp[*kf - 1] < 0)
5168 if (phase == 1 || scs_imp->blocks[*kf - 1].nmode == 0)
5170 i = synchro_nev(scs_imp, *kf, ierr);
5178 i = scs_imp->blocks[*kf - 1].mode[0];
5180 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5181 putevs(told, &i2, &ierr1);
5184 /* event conflict */
5188 ozdoit(told, xt, xtd, &kiwa);
5197 /* update zero crossing surfaces */
5198 for (ii = 0; ii < * (scs_imp->nzord); ii++)
5200 kf = &scs_imp->zord[ii];
5201 C2F(curblk).kfun = *kf;
5202 if (scs_imp->blocks[*kf - 1].ng > 0)
5204 /* update g array ptr */
5205 scs_imp->blocks[*kf - 1].g = &g[scs_imp->zcptr[*kf - 1] - 1];
5206 if (scs_imp->funtyp[*kf - 1] > 0)
5209 /* continuous state */
5210 if (scs_imp->blocks[*kf - 1].nx > 0)
5212 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
5213 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5215 scs_imp->blocks[*kf - 1].nevprt = scs_imp->zord[ii + * (scs_imp->nzord)];
5216 /*callf(told, xtd, xt, xtd,g,&flag);*/
5217 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5226 j = synchro_g_nev(scs_imp, g, *kf, ierr);
5231 if ( (phase == 1) && (scs_imp->blocks[*kf - 1].nmode > 0) )
5233 scs_imp->blocks[*kf - 1].mode[0] = j;
5237 // scs_imp->blocks[*kf-1].g = &scs_imp->g[scs_imp->zcptr[*kf]-1];
5242 for (i = 0; i < kiwa; i++)
5245 for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1; ii++)
5247 kf = &scs_imp->ordclk[ii - 1];
5248 C2F(curblk).kfun = *kf;
5249 if (scs_imp->blocks[*kf - 1].ng > 0)
5251 /* update g array ptr */
5252 scs_imp->blocks[*kf - 1].g = &g[scs_imp->zcptr[*kf - 1] - 1];
5253 if (scs_imp->funtyp[*kf - 1] > 0)
5256 /* continuous state */
5257 if (scs_imp->blocks[*kf - 1].nx > 0)
5259 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
5260 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5262 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5263 /*callf(told, xtd, xt, xtd,g,&flag);*/
5264 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5273 j = synchro_g_nev(scs_imp, g, *kf, ierr);
5278 if ((phase == 1) && (scs_imp->blocks[*kf - 1].nmode > 0))
5280 scs_imp->blocks[*kf - 1].mode[0] = j;
5284 //scs_imp->blocks[*kf-1].g = &scs_imp->g[scs_imp->zcptr[*kf]-1];
5289 /*--------------------------------------------------------------------------*/
5290 /* Subroutine Jdoit */
5291 void Jdoit(double *told, double *xt, double *xtd, double *residual, int *job)
5293 /* update blocks jacobian of continuous time block */
5294 /* Copyright INRIA */
5297 scicos_flag flag = 0;
5298 int keve = 0, kiwa = 0;
5299 int ierr1 = 0, i = 0;
5302 ScicosImport *scs_imp = NULL;
5305 scs_imp = getscicosimportptr();
5309 for (jj = 0; jj < * (scs_imp->noord); jj++)
5311 kf = &scs_imp->oord[jj];
5312 C2F(curblk).kfun = *kf;
5313 scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[jj + * (scs_imp->noord)];
5314 if (scs_imp->funtyp[*kf - 1] > -1)
5317 /* applying desired output */
5318 if ((*job == 2) && (scs_imp->oord[jj] == AJacobian_block))
5322 /* continuous state */
5323 if (scs_imp->blocks[*kf - 1].nx > 0)
5325 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
5326 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5327 scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5330 /*callf(told, xtd, xt, residual,x,&flag);*/
5331 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5339 if (scs_imp->blocks[*kf - 1].nevout > 0)
5341 if (scs_imp->funtyp[*kf - 1] < 0)
5343 if (scs_imp->blocks[*kf - 1].nmode > 0)
5345 i2 = scs_imp->blocks[*kf - 1].mode[0] + scs_imp->clkptr[*kf - 1] - 1;
5349 i = synchro_nev(scs_imp, *kf, ierr);
5354 i2 = i + scs_imp->clkptr[*kf - 1] - 1;
5356 putevs(told, &i2, &ierr1);
5359 /* event conflict */
5363 ozdoit(told, xt, xtd, &kiwa);
5368 /* update states derivatives */
5369 for (ii = 0; ii < * (scs_imp->noord); ii++)
5371 kf = &scs_imp->oord[ii];
5372 C2F(curblk).kfun = *kf;
5373 if (scs_imp->blocks[*kf - 1].nx > 0 || \
5374 *scs_imp->blocks[*kf - 1].work != NULL)
5376 /* work tests if a hidden state exists, used for delay block */
5378 if (((*job == 1) && (scs_imp->oord[ii] == AJacobian_block)) || (*job != 1))
5384 /* continuous state */
5385 if (scs_imp->blocks[*kf - 1].nx > 0)
5387 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
5388 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5389 scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5391 scs_imp->blocks[*kf - 1].nevprt = scs_imp->oord[ii + * (scs_imp->noord)];
5392 /*callf(told, xtd, xt, residual,xt,&flag);*/
5393 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5403 for (i = 0; i < kiwa; i++)
5406 for (ii = scs_imp->ordptr[keve - 1]; ii <= scs_imp->ordptr[keve] - 1; ii++)
5408 kf = &scs_imp->ordclk[ii - 1];
5409 C2F(curblk).kfun = *kf;
5410 if (scs_imp->blocks[*kf - 1].nx > 0 || \
5411 *scs_imp->blocks[*kf - 1].work != NULL)
5413 /* work tests if a hidden state exists */
5415 if (((*job == 1) && (scs_imp->oord[ii - 1] == AJacobian_block)) || (*job != 1))
5421 /* continuous state */
5422 if (scs_imp->blocks[*kf - 1].nx > 0)
5424 scs_imp->blocks[*kf - 1].x = &xt[scs_imp->xptr[*kf - 1] - 1];
5425 scs_imp->blocks[*kf - 1].xd = &xtd[scs_imp->xptr[*kf - 1] - 1];
5426 scs_imp->blocks[*kf - 1].res = &residual[scs_imp->xptr[*kf - 1] - 1];
5428 scs_imp->blocks[*kf - 1].nevprt = abs(scs_imp->ordclk[ii + * (scs_imp->nordclk) - 1]);
5429 /*callf(told, xtd, xt, residual,xt,&flag);*/
5430 callf(told, &scs_imp->blocks[*kf - 1], &flag);
5441 /*--------------------------------------------------------------------------*/
5442 /* Subroutine synchro_nev */
5443 static int synchro_nev(ScicosImport *scs_imp, int kf, int *ierr)
5445 /* synchro blocks computation */
5446 /* Copyright INRIA */
5447 SCSREAL_COP *outtbdptr = NULL; /*to store double of outtb*/
5448 SCSINT8_COP *outtbcptr = NULL; /*to store int8 of outtb*/
5449 SCSINT16_COP *outtbsptr = NULL; /*to store int16 of outtb*/
5450 SCSINT32_COP *outtblptr = NULL; /*to store int32 of outtb*/
5451 SCSUINT8_COP *outtbucptr = NULL; /*to store unsigned int8 of outtb */
5452 SCSUINT16_COP *outtbusptr = NULL; /*to store unsigned int16 of outtb */
5453 SCSUINT32_COP *outtbulptr = NULL; /*to store unsigned int32 of outtb */
5456 int i = 0; /* return 0 by default */
5458 /* variable for param */
5460 void **outtbptr = NULL;
5466 outtbtyp = scs_imp->outtbtyp;
5467 outtbptr = scs_imp->outtbptr;
5468 funtyp = scs_imp->funtyp;
5469 inplnk = scs_imp->inplnk;
5470 inpptr = scs_imp->inpptr;
5472 /* if-then-else blk */
5473 if (funtyp[kf - 1] == -1)
5475 switch (outtbtyp[-1 + inplnk[inpptr[kf - 1] - 1]])
5478 outtbdptr = (SCSREAL_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5479 cond = (*outtbdptr <= 0.);
5483 outtbdptr = (SCSCOMPLEX_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5484 cond = (*outtbdptr <= 0.);
5488 outtbcptr = (SCSINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5489 cond = (*outtbcptr <= 0);
5493 outtbsptr = (SCSINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5494 cond = (*outtbsptr <= 0);
5498 outtblptr = (SCSINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5499 cond = (*outtblptr <= 0);
5503 outtbucptr = (SCSUINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5504 cond = (*outtbucptr <= 0);
5508 outtbusptr = (SCSUINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5509 cond = (*outtbusptr <= 0);
5513 outtbulptr = (SCSUINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5514 cond = (*outtbulptr <= 0);
5517 default : /* Add a message here */
5531 else if (funtyp[kf - 1] == -2)
5533 switch (outtbtyp[-1 + inplnk[inpptr[kf - 1] - 1]])
5536 outtbdptr = (SCSREAL_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5537 i = Max(Min((int) * outtbdptr, scs_imp->blocks[kf - 1].nevout), 1);
5541 outtbdptr = (SCSCOMPLEX_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5542 i = Max(Min((int) * outtbdptr, scs_imp->blocks[kf - 1].nevout), 1);
5546 outtbcptr = (SCSINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5547 i = Max(Min((int) * outtbcptr,
5548 scs_imp->blocks[kf - 1].nevout), 1);
5552 outtbsptr = (SCSINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5553 i = Max(Min((int) * outtbsptr, scs_imp->blocks[kf - 1].nevout), 1);
5557 outtblptr = (SCSINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5558 i = Max(Min((int) * outtblptr, scs_imp->blocks[kf - 1].nevout), 1);
5562 outtbucptr = (SCSUINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5563 i = Max(Min((int) * outtbucptr, scs_imp->blocks[kf - 1].nevout), 1);
5567 outtbusptr = (SCSUINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5568 i = Max(Min((int) * outtbusptr, scs_imp->blocks[kf - 1].nevout), 1);
5572 outtbulptr = (SCSUINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5573 i = Max(Min((int) * outtbulptr, scs_imp->blocks[kf - 1].nevout), 1);
5576 default : /* Add a message here */
5583 /*--------------------------------------------------------------------------*/
5584 /* Subroutine synchro_g_nev */
5585 static int synchro_g_nev(ScicosImport *scs_imp, double *g, int kf, int *ierr)
5587 /* synchro blocks with zcross computation */
5588 /* Copyright INRIA */
5589 SCSREAL_COP *outtbdptr = NULL; /*to store double of outtb*/
5590 SCSINT8_COP *outtbcptr = NULL; /*to store int8 of outtb*/
5591 SCSINT16_COP *outtbsptr = NULL; /*to store int16 of outtb*/
5592 SCSINT32_COP *outtblptr = NULL; /*to store int32 of outtb*/
5593 SCSUINT8_COP *outtbucptr = NULL; /*to store unsigned int8 of outtb */
5594 SCSUINT16_COP *outtbusptr = NULL; /*to store unsigned int16 of outtb */
5595 SCSUINT32_COP *outtbulptr = NULL; /*to store unsigned int32 of outtb */
5598 int i = 0; /* return 0 by default */
5601 /* variable for param */
5602 int *outtbtyp = NULL;
5603 void **outtbptr = NULL;
5610 outtbtyp = scs_imp->outtbtyp;
5611 outtbptr = scs_imp->outtbptr;
5612 funtyp = scs_imp->funtyp;
5613 inplnk = scs_imp->inplnk;
5614 inpptr = scs_imp->inpptr;
5615 zcptr = scs_imp->zcptr;
5617 /* if-then-else blk */
5618 if (funtyp[kf - 1] == -1)
5620 switch (outtbtyp[-1 + inplnk[inpptr[kf - 1] - 1]])
5623 outtbdptr = (SCSREAL_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5624 g[zcptr[kf - 1] - 1] = *outtbdptr;
5625 cond = (*outtbdptr <= 0.);
5629 outtbdptr = (SCSCOMPLEX_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5630 g[zcptr[kf - 1] - 1] = *outtbdptr;
5631 cond = (*outtbdptr <= 0.);
5635 outtbcptr = (SCSINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5636 g[zcptr[kf - 1] - 1] = (double) * outtbcptr;
5637 cond = (*outtbcptr <= 0);
5641 outtbsptr = (SCSINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5642 g[zcptr[kf - 1] - 1] = (double) * outtbsptr;
5643 cond = (*outtbsptr <= 0);
5647 outtblptr = (SCSINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5648 g[zcptr[kf - 1] - 1] = (double) * outtblptr;
5649 cond = (*outtblptr <= 0);
5653 outtbucptr = (SCSUINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5654 g[zcptr[kf - 1] - 1] = (double) * outtbucptr;
5655 cond = (*outtbucptr <= 0);
5659 outtbusptr = (SCSUINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5660 g[zcptr[kf - 1] - 1] = (double) * outtbusptr;
5661 cond = (*outtbusptr <= 0);
5665 outtbulptr = (SCSUINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5666 g[zcptr[kf - 1] - 1] = (double) * outtbulptr;
5667 cond = (*outtbulptr <= 0);
5670 default : /* Add a message here */
5684 else if (funtyp[kf - 1] == -2)
5686 switch (outtbtyp[-1 + inplnk[inpptr[kf - 1] - 1]])
5689 outtbdptr = (SCSREAL_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5690 for (jj = 0; jj < scs_imp->blocks[kf - 1].nevout - 1; jj++)
5692 g[zcptr[kf - 1] - 1 + jj] = *outtbdptr - (double)(jj + 2);
5694 i = Max(Min((int) * outtbdptr, scs_imp->blocks[kf - 1].nevout), 1);
5698 outtbdptr = (SCSCOMPLEX_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5699 for (jj = 0; jj < scs_imp->blocks[kf - 1].nevout - 1; jj++)
5701 g[zcptr[kf - 1] - 1 + jj] = *outtbdptr - (double)(jj + 2);
5703 i = Max(Min((int) * outtbdptr, scs_imp->blocks[kf - 1].nevout), 1);
5707 outtbcptr = (SCSINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5708 for (jj = 0; jj < scs_imp->blocks[kf - 1].nevout - 1; jj++)
5710 g[zcptr[kf - 1] - 1 + jj] = (double) * outtbcptr - (double)(jj + 2);
5712 i = Max(Min((int) * outtbcptr, scs_imp->blocks[kf - 1].nevout), 1);
5716 outtbsptr = (SCSINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5717 for (jj = 0; jj < scs_imp->blocks[kf - 1].nevout - 1; jj++)
5719 g[zcptr[kf - 1] - 1 + jj] = (double) * outtbsptr - (double)(jj + 2);
5721 i = Max(Min((int) * outtbsptr, scs_imp->blocks[kf - 1].nevout), 1);
5725 outtblptr = (SCSINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5726 for (jj = 0; jj < scs_imp->blocks[kf - 1].nevout - 1; jj++)
5728 g[zcptr[kf - 1] - 1 + jj] = (double) * outtblptr - (double)(jj + 2);
5730 i = Max(Min((int) * outtblptr, scs_imp->blocks[kf - 1].nevout), 1);
5734 outtbucptr = (SCSUINT8_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5735 for (jj = 0; jj < scs_imp->blocks[kf - 1].nevout - 1; jj++)
5737 g[zcptr[kf - 1] - 1 + jj] = (double) * outtbucptr - (double)(jj + 2);
5739 i = Max(Min((int) * outtbucptr, scs_imp->blocks[kf - 1].nevout), 1);
5743 outtbusptr = (SCSUINT16_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5744 for (jj = 0; jj < scs_imp->blocks[kf - 1].nevout - 1; jj++)
5746 g[zcptr[kf - 1] - 1 + jj] = (double) * outtbusptr - (double)(jj + 2);
5748 i = Max(Min((int) * outtbusptr, scs_imp->blocks[kf - 1].nevout), 1);
5752 outtbulptr = (SCSUINT32_COP *)outtbptr[-1 + inplnk[inpptr[kf - 1] - 1]];
5753 for (jj = 0; jj < scs_imp->blocks[kf - 1].nevout - 1; jj++)
5755 g[zcptr[kf - 1] - 1 + jj] = (double) * outtbulptr - (double)(jj + 2);
5757 i = Max(Min((int) * outtbulptr, scs_imp->blocks[kf - 1].nevout), 1);
5760 default : /* Add a message here */
5766 } /* synchro_g_nev */
5767 /*--------------------------------------------------------------------------*/
5769 static void FREE_blocks()
5772 for (kf = 0; kf < nblk; ++kf)