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