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