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