scilab/modules/scicos/src/c/scicos.c

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