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

   1 /*
   2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   3  * Copyright (C) Scilab Enterprises - 2013 - Paul Bignier
   4  *
   5  * This file must be used under the terms of the CeCILL.
   6  * This source file is licensed as described in the file COPYING, which
   7  * you should have received as part of this distribution.  The terms
   8  * are also available at
   9  * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
  10  *
  11  */
  12
  13 #include <stdio.h>
  14 #include <stdlib.h>
  15 #include <stdarg.h>
  16 #include <string.h>
  17 #include <machine.h>
  18 #include <float.h>   // DBL_EPSILON, to define UNIT_ROUNDOFF
  19 #include <math.h>    // fabs() and pow() functions
  20
  21 #include "ddaskr.h"
  22
  23 #define NO_FPRINTF_OUTPUT 1
  24
  25 /*
  26  * Control constant for tolerances
  27  *
  28  * Scicos only uses scalar tolerances, so we only need the scalar-scalar (SS) value for info[1].
  29  * --------------------------------
  30  */
  31 #define DDAS_SS           0
  32
  33 /* Flags for stop mode */
  34 #define DDAS_NORMAL       0
  35 #define DDAS_ONE_STEP     1
  36
  37 /*
  38  * Flags for initial values calculation.
  39  *
  40  * Indicating whether the components have been qualified (differential / algebraic).
  41  * --------------------------------
  42  */
  43 #define DDAS_YA_YDP_INIT  1
  44 #define DDAS_Y_INIT       2
  45
  46 #define UNIT_ROUNDOFF     DBL_EPSILON
  47
  48 /* =============================
  49  *
  50  *           ddaskr
  51  *
  52  * =============================
  53  *
  54  * Actual solving function, from 'ODEPACK' in 'differential_equations' module.
  55  * Since we use ddaskr's built-in jacobian function, set jacpsol type to DDasJacPsolFn.
  56  */
  57
  58 extern void C2F(ddaskr) (DDASResFn res, int *neq, realtype *t, realtype *y, realtype *yp, realtype *tout, int *info, realtype *reltol, realtype *abstol, int *istate, struct DDrWork_t *rwork, int *lrw, int *iwork, int *liw, double *dummy1, int *dummy2, DDASJacPsolFn jacpsol, DDASPsolFn psol, DDASRootFn grblk, int *ng, int *jroot);
  59
  60 /* =============================
  61  *
  62  *         DDaskrCreate
  63  *
  64  * =============================
  65  *
  66  * DDaskrCreate creates an internal memory block for a problem to be solved by DDASKR.
  67  * If successful, DDaskrCreate returns a pointer to the problem memory.
  68  * This pointer should be passed to DDaskrInit.
  69  * If an initialization error occurs,
  70  * DDaskrCreate prints an error message to standard err and returns NULL.
  71  */
  72
  73 void * DDaskrCreate (int * neq, int ng, int solverIndex)
  74 {
  75     int lIw = 0, lRw = 0, LENWP = 0, LENIWP = 0;
  76     DDaskrMem ddaskr_mem = NULL;
  77
  78     /* Allocate the problem memory space */
  79     ddaskr_mem = NULL;
  80     ddaskr_mem = (DDaskrMem) malloc(sizeof(struct DDaskrMemRec));
  81     if (ddaskr_mem == NULL)
  82     {
  83         DDASProcessError(NULL, 0, "DDASKR", "DDaskrCreate", MSG_MEM_FAIL);
  84         return (NULL);
  85     }
  86
  87     /* Zero out ddas_mem */
  88     memset(ddaskr_mem, 0, sizeof(struct DDaskrMemRec));
  89
  90     /* Set the 'rwork' and 'iwork' workspaces lengths by default
  91        LENWP and LENIWP are lentghts of segments of rwork and iwork (respectively),
  92        that will contain factored preconditioner matrix information */
  93     LENWP  = (*neq) * (*neq);
  94     LENIWP = (*neq);
  95     lRw    = 60 + (*neq) * (max(MAXORD_DEFAULT + 4, 7) + (*neq)) + 3 * ng;
  96     lIw    = 40 + 2 * (*neq);
  97
  98     /* If we are going to use the Krylov method, resize the workspaces adequately */
  99     if (solverIndex == 102)
 100     {
 101         lRw = 101 + 18 * (*neq) + 3 * ng + LENWP;
 102         lIw = 40 + (*neq) + LENIWP;
 103     }
 104
 105     /* Copy the variables into the problem memory space */
 106     ddaskr_mem->nEquations = neq;
 107     ddaskr_mem->user_data  = NULL;
 108     ddaskr_mem->iState     = 0;
 109     ddaskr_mem->info       = NULL;
 110     ddaskr_mem->rwork      = NULL;
 111     ddaskr_mem->lrw        = lRw;
 112     ddaskr_mem->iwork      = NULL;
 113     ddaskr_mem->liw        = lIw;
 114     ddaskr_mem->ehfun      = NULL;
 115     ddaskr_mem->g_fun      = NULL;
 116     ddaskr_mem->ng_fun     = ng;
 117     ddaskr_mem->jroot      = NULL;
 118     ddaskr_mem->solver     = solverIndex;
 119     ddaskr_mem->jacpsol    = NULL;
 120     ddaskr_mem->psol       = NULL;
 121     ddaskr_mem->rpar       = NULL;
 122     ddaskr_mem->ipar       = NULL;
 123
 124     return ((void *) ddaskr_mem);
 125 }
 126
 127 /* Shortcuts to problem memory space parameters */
 128 # define res        ddas_mem->res
 129 # define nEq        ddas_mem->nEquations
 130 # define user_data  ddas_mem->user_data
 131 # define yVec       ddas_mem->yVector
 132 # define ypVec      ddas_mem->yPrimeVector
 133 # define tStart     ddas_mem->tStart
 134 # define info       ddas_mem->info
 135 # define relTol     ddas_mem->relTol
 136 # define absTol     ddas_mem->absTol
 137 # define iState     ddas_mem->iState
 138 # define rwork      ddas_mem->rwork
 139 # define lrw        ddas_mem->lrw
 140 # define iwork      ddas_mem->iwork
 141 # define liw        ddas_mem->liw
 142 # define maxnhIC    ddas_mem->maxnhIC
 143 # define g_fun      ddas_mem->g_fun
 144 # define ng_fun     ddas_mem->ng_fun
 145 # define jroot      ddas_mem->jroot
 146 # define solver     ddas_mem->solver
 147 # define jacpsol    ddas_mem->jacpsol
 148 # define psol       ddas_mem->psol
 149 # define rpar       ddas_mem->rpar
 150 # define ipar       ddas_mem->ipar
 151
 152 /* =============================
 153  *
 154  *          DDaskrInit
 155  *
 156  * =============================
 157  *
 158  * DDaskrInit allocates and initializes memory for a problem.
 159  * All problem inputs are checked for errors. If any error occurs during initialization,
 160  * it is reported to the file whose file pointer is errfp and an error flag is returned.
 161  * Otherwise, it returns IDA_SUCCESS.
 162  */
 163
 164 int DDaskrInit (void * ddaskr_mem, DDASResFn Res, realtype t0, N_Vector yy0, N_Vector yp0, DDASJacPsolFn Jacpsol, DDASPsolFn Psol)
 165 {
 166     DDaskrMem ddas_mem = NULL;
 167
 168     /* Check the input arguments */
 169
 170     if (ddaskr_mem == NULL)
 171     {
 172         DDASProcessError(NULL, IDA_MEM_NULL, "DDASKR", "DDaskrInit", MSG_NO_MEM);
 173         return (IDA_MEM_NULL);
 174     }
 175     ddas_mem = (DDaskrMem) ddaskr_mem;
 176
 177     if (yy0 == NULL)
 178     {
 179         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrInit", MSG_Y0_NULL);
 180         return (IDA_ILL_INPUT);
 181     }
 182
 183     if (yp0 == NULL)
 184     {
 185         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrInit", MSG_YP0_NULL);
 186         return (IDA_ILL_INPUT);
 187     }
 188
 189     if (Res == NULL)
 190     {
 191         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrInit", MSG_RES_NULL);
 192         return (IDA_ILL_INPUT);
 193     }
 194
 195     /* Jacpsol = NULL or Psol = NULL is a problem only if the user decided to use the GMRes solver */
 196     if (solver == 102)
 197     {
 198         if (Jacpsol == NULL || Psol == NULL)
 199         {
 200             DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrInit", MSG_BAD_KRY_INPUT);
 201             return (IDA_ILL_INPUT);
 202         }
 203     }
 204
 205     /* Copy the arguments into the problem memory space */
 206     res     = Res;
 207     yVec    = NV_DATA_S(yy0);
 208     ypVec   = NV_DATA_S(yp0);
 209     tStart  = t0;
 210     jacpsol = Jacpsol;
 211     psol    = Psol;
 212
 213     /* Allocate the info[20] tab to zero, used to store parameters (zero is default value for mostof them) */
 214     info = calloc(20, sizeof(int));
 215
 216     /* info[11] = 1 => Krylov method selected
 217        info[14] = 1 => providing jacobian function (to evaluate and LU-factor the preconditioner) */
 218     if (solver == 102)
 219     {
 220         info[11] = 1;
 221         info[14] = 1;
 222     }
 223
 224     /* Allocate rwork and iwork workspaces and set them to zero.
 225        Their size is lrw and liw, respectively */
 226     rwork = (struct DDrWork_t *) calloc(lrw, sizeof(realtype));
 227     iwork = calloc(liw, sizeof(int));
 228
 229     /* Save their lengths in iwork */
 230     iwork[16] = lrw;
 231     iwork[17] = liw;
 232
 233     /* Solve the problem without invoking any special inequality constraints */
 234     info[9] = 0;
 235
 236     /* Default values for heuristic control quantities in the initial values calculation */
 237     iwork[31] = (info[11] == 0) ? 5 : 15;               // maxnit, depends on dense / Krylov method
 238     iwork[32] = (info[11] == 0) ? 6  : 2;               // maxnj, depends on dense / Krylov method
 239     iwork[33] = 5;                                                              // maxnh
 240     iwork[34] = 0;                                                              // lsoff, flag to turn linesearch on / off (0 = on)
 241     rwork->steptol = pow(UNIT_ROUNDOFF, 2. / 3);        // steptol
 242     rwork->epinit = 0.01;                                               // epinit, swing factor in the Newton iteration convergence test.
 243     maxnhIC = 5;                                                                // maxnh for IC calculation
 244
 245     return (IDA_SUCCESS);
 246 }
 247
 248 /* =============================
 249  *
 250  *         DDaskrReInit
 251  *
 252  * =============================
 253  *
 254  * DDaskrReInit reinitializes DDASKR's memory for a problem,
 255  * assuming it has already been allocated in a prior DDaskrInit call.
 256  * All problem specification inputs are checked for errors.
 257  * If any error occurs during initialization, it is reported to the file whose file pointer is errfp.
 258  * The return value is IDA_SUCCESS = 0 if no errors occurred, or a negative value otherwise.
 259  */
 260
 261 int DDaskrReInit (void * ddaskr_mem, realtype tOld, N_Vector yy0, N_Vector yp0)
 262 {
 263     DDaskrMem ddas_mem = NULL;
 264
 265     /* Check the input arguments */
 266
 267     if (ddaskr_mem == NULL)
 268     {
 269         DDASProcessError(NULL, IDA_MEM_NULL, "DDASKR", "DDaskrReInit", MSG_NO_MEM);
 270         return (IDA_MEM_NULL);
 271     }
 272     ddas_mem = (DDaskrMem) ddaskr_mem;
 273
 274     if (yy0 == NULL)
 275     {
 276         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrInit", MSG_Y0_NULL);
 277         return (IDA_ILL_INPUT);
 278     }
 279
 280     if (yp0 == NULL)
 281     {
 282         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrInit", MSG_YP0_NULL);
 283         return (IDA_ILL_INPUT);
 284     }
 285
 286     /* Reset the problem memory space variables to the arguments */
 287     *nEq    = NV_LENGTH_S(yy0);
 288     yVec    = NV_DATA_S(yy0);
 289     ypVec   = NV_DATA_S(yp0);
 290     tStart  = tOld;
 291     iState  = 0;
 292     info[0] = 0;
 293
 294     /* Tell DDaskr to get new consitent values */
 295     info[10] = 1;
 296
 297     return (IDA_SUCCESS);
 298 }
 299
 300 /* =============================
 301  *
 302  *       DDaskrSStolerances
 303  *
 304  * =============================
 305  *
 306  * This function specifies the scalar integration tolerances.
 307  * It MUST be called before the first call to DDaskr.
 308  */
 309
 310 int DDaskrSStolerances (void * ddaskr_mem, realtype reltol, realtype abstol)
 311 {
 312     DDaskrMem ddas_mem = NULL;
 313
 314     if (ddaskr_mem == NULL)
 315     {
 316         DDASProcessError(NULL, IDA_MEM_NULL, "DDaskr", "DDaskrSStolerances", MSG_NO_MEM);
 317         return (IDA_MEM_NULL);
 318     }
 319     ddas_mem = (DDaskrMem) ddaskr_mem;
 320
 321     /* Check inputs */
 322
 323     if (reltol < 0.)
 324     {
 325         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrSStolerances", MSG_BAD_RTOL);
 326         return (IDA_ILL_INPUT);
 327     }
 328
 329     if (abstol < 0.)
 330     {
 331         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrSStolerances", MSG_BAD_ATOL);
 332         return (IDA_ILL_INPUT);
 333     }
 334
 335     /* Copy tolerances into memory */
 336
 337     relTol  = reltol;
 338     absTol  = abstol;
 339     info[1] = DDAS_SS;
 340
 341     return (IDA_SUCCESS);
 342 }
 343
 344 /* =============================
 345  *
 346  *         DDaskrRootInit
 347  *
 348  * =============================
 349  *
 350  * DDaskrRootInit initializes a rootfinding problem to be solved during the integration of the ODE system.
 351  * It loads the root function pointer and the number of root functions, and allocates workspace memory.
 352  * The return value is IDA_SUCCESS = 0 if no errors occurred, or a negative value otherwise.
 353  */
 354
 355 int DDaskrRootInit (void * ddaskr_mem, int ng, DDASRootFn g)
 356 {
 357     DDaskrMem ddas_mem = NULL;
 358     int nrt = 0;
 359
 360     if (ddaskr_mem == NULL)
 361     {
 362         DDASProcessError(NULL, IDA_MEM_NULL, "DDASKR", "DDaskrRootInit", MSG_NO_MEM);
 363         return (IDA_MEM_NULL);
 364     }
 365     ddas_mem = (DDaskrMem) ddaskr_mem;
 366
 367     if (g == NULL)
 368     {
 369         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrRootInit", MSG_ROOT_FUNC_NULL);
 370         return (IDA_ILL_INPUT);
 371     }
 372
 373     g_fun  = g;
 374     nrt    = (ng < 0) ? 0 : ng;
 375     ng_fun = nrt;
 376
 377     /* Allocate jroot and set it to zero */
 378     if (ng > 0)
 379     {
 380         jroot = calloc(ng, sizeof(int));
 381     }
 382
 383     return (IDA_SUCCESS);
 384 }
 385
 386 /* =============================
 387  *
 388  *       DDaskrSetUserData
 389  *
 390  * =============================
 391  *
 392  * Sets a pointer to user_data that will be passed to the user's res function
 393  * every time a user-supplied function is called.
 394  */
 395
 396 int DDaskrSetUserData (void * ddaskr_mem, void * User_data)
 397 {
 398     DDaskrMem ddas_mem = NULL;
 399
 400     if (ddaskr_mem == NULL)
 401     {
 402         DDASProcessError(NULL, IDA_MEM_NULL, "DDASKR", "DDaskrSetUserData", MSG_NO_MEM);
 403         return(IDA_MEM_NULL);
 404     }
 405     ddas_mem = (DDaskrMem) ddaskr_mem;
 406
 407     user_data = User_data;
 408
 409     return(IDA_SUCCESS);
 410 }
 411
 412 /* =============================
 413  *
 414  *       DDaskrSetMaxStep
 415  *
 416  * =============================
 417  *
 418  * Sets the maximum step size, stocked in rwork->hmax.
 419  * Sets info[6] to 1 for rwork->hmax to be taken in consideration by ddaskr().
 420  */
 421
 422 int DDaskrSetMaxStep (void * ddaskr_mem, realtype hMax)
 423 {
 424     DDaskrMem ddas_mem = NULL;
 425
 426     if (ddaskr_mem == NULL)
 427     {
 428         DDASProcessError(NULL, IDA_MEM_NULL, "DDASKR", "DDaskrSetMaxStep", MSG_NO_MEM);
 429         return (IDA_MEM_NULL);
 430     }
 431     ddas_mem = (DDaskrMem) ddaskr_mem;
 432
 433     if (info[6] == 0)
 434     {
 435         info[6] = 1;
 436     }
 437     rwork->hmax = hMax;
 438
 439     return (IDA_SUCCESS);
 440 }
 441
 442 /* =============================
 443  *
 444  *       DDaskrSetStopTime
 445  *
 446  * =============================
 447  *
 448  * Specifies the time beyond which the integration is not to proceed, stocked in rwork->tcrit.
 449  * Sets info[3] to 1 for rwork->tcrit to be taken in consideration by ddaskr().
 450  */
 451
 452 int DDaskrSetStopTime (void * ddaskr_mem, realtype tCrit)
 453 {
 454     DDaskrMem ddas_mem = NULL;
 455
 456     if (ddaskr_mem == NULL)
 457     {
 458         DDASProcessError(NULL, IDA_MEM_NULL, "DDASKR", "DDaskrSetStopTime", MSG_NO_MEM);
 459         return (IDA_MEM_NULL);
 460     }
 461     ddas_mem = (DDaskrMem) ddaskr_mem;
 462
 463     if (info[3] == 0)
 464     {
 465         info[3] = 1;
 466     }
 467
 468     rwork->tcrit = tCrit;
 469
 470     return (IDA_SUCCESS);
 471 }
 472
 473 /* =============================
 474  *
 475  *     DDaskrSetMaxNumSteps
 476  *
 477  * =============================
 478  *
 479  * Sets the maximum number of steps in an integration interval.
 480  * Ensure that ddaskr will consider it via flag info[16], and stock it in iwork[33].
 481  */
 482
 483 int DDaskrSetMaxNumSteps (void * ddaskr_mem, int maxnh)
 484 {
 485     DDaskrMem ddas_mem = NULL;
 486
 487     if (ddaskr_mem == NULL)
 488     {
 489         DDASProcessError(NULL, IDA_MEM_NULL, "DDASKR", "DDaskrSetMaxNumSteps", MSG_NO_MEM);
 490         return (IDA_MEM_NULL);
 491     }
 492     ddas_mem = (DDaskrMem) ddaskr_mem;
 493
 494     if (maxnh <= 0)
 495     {
 496         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "IDA", "DDaskrSetMaxNumSteps", MSG_BAD_MAXNH);
 497         return (IDA_ILL_INPUT);
 498     }
 499
 500     if (info[16] == 0)
 501     {
 502         info[16] = 1;
 503     }
 504
 505     iwork[33] = maxnh;
 506
 507     return (IDA_SUCCESS);
 508 }
 509
 510 /* =============================
 511  *
 512  *     DDaskrSetMaxNumJacsIC
 513  *
 514  * =============================
 515  *
 516  * Sets the maximum number of Jacobian or preconditioner evaluations.
 517  * Ensure that ddaskr will consider it via flag info[16], and stock it in iwork[32].
 518  */
 519
 520 int DDaskrSetMaxNumJacsIC (void * ddaskr_mem, int maxnj)
 521 {
 522     DDaskrMem ddas_mem = NULL;
 523
 524     if (ddaskr_mem == NULL)
 525     {
 526         DDASProcessError(NULL, IDA_MEM_NULL, "DDASKR", "DDaskrSetMaxNumJacsIC", MSG_NO_MEM);
 527         return (IDA_MEM_NULL);
 528     }
 529     ddas_mem = (DDaskrMem) ddaskr_mem;
 530
 531     if (maxnj <= 0)
 532     {
 533         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "IDA", "DDaskrSetMaxNumJacsIC", MSG_BAD_MAXNJ);
 534         return (IDA_ILL_INPUT);
 535     }
 536
 537     if (info[16] == 0)
 538     {
 539         info[16] = 1;
 540     }
 541
 542     iwork[32] = maxnj;
 543
 544     return (IDA_SUCCESS);
 545 }
 546
 547 /* =============================
 548  *
 549  *    DDaskrSetMaxNumItersIC
 550  *
 551  * =============================
 552  *
 553  * Sets the maximum number of Newton iterations per Jacobian or preconditioner evaluation.
 554  * Ensure that ddaskr will consider it via flag info[16], and stock it in iwork[31].
 555  */
 556
 557 int DDaskrSetMaxNumItersIC (void * ddaskr_mem, int maxnit)
 558 {
 559     DDaskrMem ddas_mem = NULL;
 560
 561     if (ddaskr_mem == NULL)
 562     {
 563         DDASProcessError(NULL, IDA_MEM_NULL, "DDASKR", "DDaskrSetMaxNumItersIC", MSG_NO_MEM);
 564         return (IDA_MEM_NULL);
 565     }
 566     ddas_mem = (DDaskrMem) ddaskr_mem;
 567
 568     if (maxnit <= 0)
 569     {
 570         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "IDA", "DDaskrSetMaxNumItersIC", MSG_BAD_MAXNIT);
 571         return (IDA_ILL_INPUT);
 572     }
 573
 574     if (info[16] == 0)
 575     {
 576         info[16] = 1;
 577     }
 578
 579     iwork[31] = maxnit;
 580
 581     return (IDA_SUCCESS);
 582 }
 583
 584 /* =============================
 585  *
 586  *    DDaskrSetMaxNumStepsIC
 587  *
 588  * =============================
 589  *
 590  * Sets the maximum number of values of the artificial stepsize parameter H to be tried if info[10] = 1,
 591  * during the Initial Conditions calculation.
 592  * Ensure that ddaskr will consider it via flag info[16], and stock it in iwork[33].
 593  */
 594
 595 int DDaskrSetMaxNumStepsIC (void * ddaskr_mem, int MaxnhIC)
 596 {
 597     DDaskrMem ddas_mem = NULL;
 598
 599     if (ddaskr_mem == NULL)
 600     {
 601         DDASProcessError(NULL, IDA_MEM_NULL, "DDASKR", "DDaskrSetMaxNumStepsIC", MSG_NO_MEM);
 602         return (IDA_MEM_NULL);
 603     }
 604     ddas_mem = (DDaskrMem) ddaskr_mem;
 605
 606     if (MaxnhIC <= 0)
 607     {
 608         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "IDA", "DDaskrSetMaxNumStepsIC", MSG_BAD_MAXNH);
 609         return (IDA_ILL_INPUT);
 610     }
 611
 612     if (info[16] == 0)
 613     {
 614         info[16] = 1;
 615     }
 616
 617     maxnhIC = MaxnhIC;
 618
 619     return (IDA_SUCCESS);
 620 }
 621
 622 /* =============================
 623  *
 624  *   DDaskrSetLineSearchOffIC
 625  *
 626  * =============================
 627  *
 628  * Sets the flag to turn off the linesearch algorithm.
 629  * lsoff = 0 means linesearch is on, lsoff = 1 means it is turned off.
 630  * Ensure that ddaskr will consider it via flag info[16], and stock it in iwork[34].
 631  */
 632
 633 int DDaskrSetLineSearchOffIC (void * ddaskr_mem, int lsoff)
 634 {
 635     DDaskrMem ddas_mem = NULL;
 636
 637     if (ddaskr_mem == NULL)
 638     {
 639         DDASProcessError(NULL, IDA_MEM_NULL, "DDASKR", "DDaskrSetLineSearchOffIC", MSG_NO_MEM);
 640         return (IDA_MEM_NULL);
 641     }
 642     ddas_mem = (DDaskrMem) ddaskr_mem;
 643
 644     if (info[16] == 0)
 645     {
 646         info[16] = 1;
 647     }
 648
 649     if (lsoff)
 650     {
 651         iwork[34] = 1;
 652     }
 653
 654     return (IDA_SUCCESS);
 655 }
 656
 657 /* =============================
 658  *
 659  *          DDaskrSetID
 660  *
 661  * =============================
 662  *
 663  * Specifies which components are differential and which ones are algrebraic, in order to get consistent initial values.
 664  * They are stocked in xproperty[neq] in the form:
 665  *  - xproperty[i] = 1 => differential component
 666  *  - xproperty[i] = 0 => algebraic component
 667  */
 668
 669 int DDaskrSetId (void * ddaskr_mem, N_Vector xproperty)
 670 {
 671     DDaskrMem ddas_mem = NULL;
 672     realtype * temp = NULL;
 673     int i = 0, LID = 0;
 674
 675     if (ddaskr_mem == NULL)
 676     {
 677         DDASProcessError(NULL, IDA_MEM_NULL, "DDASKR", "DDaskrSetID", MSG_NO_MEM);
 678         return (IDA_MEM_NULL);
 679     }
 680     ddas_mem = (DDaskrMem) ddaskr_mem;
 681
 682     /* Copy xproperty data */
 683     temp = NULL;
 684     temp = NV_DATA_S(xproperty);
 685
 686     /* Inform ddaskr that we are going to define components
 687        This is also a flag to compute consistent initial values */
 688     if (info[10] == 0)
 689     {
 690         info[10] = 1;
 691     }
 692
 693     /* Since we don't use the non-negative constraints in scicos, LID = 40 (otherwise, LID = 40+neq) */
 694     LID = (info[9] == 0) ? 40 : 40 + *nEq;
 695
 696     /* Stock xproperty in a segment of iwork
 697        temp[i] = 0 or 1, but that segment needs to be -1 or 1 (respectively) */
 698     for (i = 0; i < *nEq; ++i)
 699     {
 700         iwork[i + LID] = (temp[i] == 0) ? -1 : 1;
 701     }
 702
 703     return (IDA_SUCCESS);
 704 }
 705
 706 /* =============================
 707  *
 708  *          DDaskrSolve
 709  *
 710  * =============================
 711  *
 712  * This routine is the main driver of DDASKR.
 713  *
 714  * It integrates and looks for roots over a time interval defined by the user.
 715  *
 716  * The first time that DDaskr is called for a successfully initialized problem,
 717  * it computes a tentative initial step size h.
 718  *
 719  * DDaskr supports five modes, specified by itask: DDAS_NORMAL, DDAS_ONE_STEP, DDAS_MESH_STEP, DDAS_NORMAL_TSTOP, and DDAS_ONE_STEP_TSTOP.
 720  *
 721  * In the DDAS_NORMAL mode, the solver steps until it reaches or passes tout and then interpolates to obtain y(tOut).
 722  * In the DDAS_ONE_STEP mode, it takes one internal step and returns.
 723  * DDAS_MESH_STEP means stop at the first internal mesh point at or beyond t = tout and return.
 724  * DDAS_NORMAL_TSTOP and DDAS_ONE_STEP_TSTOP are similar to DDAS_NORMAL and DDAS_ONE_STEP, respectively,
 725  * but the integration never proceeds past tcrit (which must have been defined through a call to DDaskrSetStopTime).
 726  *
 727  * It returns IDA_ROOT_RETURN if a root was detected, IDA_SUCCESS if the integration went correctly,
 728  * or a corresponding error flag.
 729  */
 730
 731 int DDaskrSolve (void * ddaskr_mem, realtype tOut, realtype * tOld, N_Vector yOut, N_Vector ypOut, int itask)
 732 {
 733     DDaskrMem ddas_mem = NULL;
 734
 735     /* Check the input arguments */
 736
 737     if (ddaskr_mem == NULL)
 738     {
 739         DDASProcessError(NULL, IDA_MEM_NULL, "DDASKR", "DDaskrSolve", MSG_NO_MEM);
 740         return (IDA_MEM_NULL);
 741     }
 742     ddas_mem = (DDaskrMem) ddaskr_mem;
 743
 744     if (yOut == NULL)
 745     {
 746         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrSolve", MSG_YRET_NULL);
 747         return (IDA_ILL_INPUT);
 748     }
 749
 750     if (ypOut == NULL)
 751     {
 752         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrSolve", MSG_YPRET_NULL);
 753         return (IDA_ILL_INPUT);
 754     }
 755
 756     if (tOld == NULL)
 757     {
 758         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrSolve", MSG_TRET_NULL);
 759         return(IDA_ILL_INPUT);
 760     }
 761
 762     if ((itask != DDAS_NORMAL) && (itask != DDAS_ONE_STEP))
 763     {
 764         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrSolve", MSG_BAD_ITASK);
 765         return (IDA_ILL_INPUT);
 766     }
 767
 768     /* Retrieve nEq if it has changed, use a copy of the solution vector and stock the simulation times */
 769     *nEq   = NV_LENGTH_S(yOut);
 770     yVec   = NV_DATA_S(yOut);
 771     ypVec  = NV_DATA_S(ypOut);
 772     tStart = *tOld;
 773
 774     /* Save the task mode in info[2] */
 775     info[2] = itask;
 776
 777     /* Launch the simulation with the memory space parameters.
 778        ddaskr() will update yVec, iState, rwork, iwork and jroot */
 779     C2F(ddaskr) (res, nEq, &tStart, yVec, ypVec, &tOut, info, &relTol, &absTol, &iState, rwork, &lrw, iwork, &liw, rpar, ipar, jacpsol, psol, g_fun, &ng_fun, jroot);
 780
 781     /* Increment the start time */
 782     *tOld  = tStart;
 783
 784     /* For continuation calls, avoiding recomputation of consistent values (if info[10] used to be 1) */
 785     info[10] = 0;
 786
 787     /* ddaskr() stocked the completion status in iState; return accordingly */
 788     switch (iState)
 789     {
 790         case 5:
 791             return (IDA_ROOT_RETURN);
 792         case 6:
 793             return (IDA_ZERO_DETACH_RETURN);
 794         case -1:
 795             DDASProcessError(ddas_mem, IDA_TOO_MUCH_WORK, "DDASKR", "DDaskrSolve", MSG_MAX_STEPS, tStart);
 796             return (IDA_TOO_MUCH_WORK);
 797         case -2:
 798             DDASProcessError(ddas_mem, IDA_TOO_MUCH_ACC, "DDASKR", "DDaskrSolve", MSG_TOO_MUCH_ACC, tStart);
 799             return (IDA_TOO_MUCH_ACC);
 800         case -3:
 801             DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrSolve", MSG_BAD_ATOL, tStart);
 802             return (IDA_ILL_INPUT);
 803         case -6:
 804             DDASProcessError(ddas_mem, IDA_ERR_FAIL, "DDASKR", "DDaskrSolve", MSG_ERR_FAILS, tStart);
 805             return (IDA_ERR_FAIL);
 806         case -7:
 807             DDASProcessError(ddas_mem, IDA_CONV_FAIL, "DDASKR", "DDaskrSolve", MSG_CONV_FAILS, tStart);
 808             return (IDA_CONV_FAIL);
 809         case -8:
 810             DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrSolve", MSG_SINGULAR);
 811             return (IDA_ILL_INPUT);
 812         case -9:
 813             DDASProcessError(ddas_mem, IDA_CONV_FAIL, "DDASKR", "DDaskrSolve", MSG_CONV_FAILS, tStart);
 814             return (IDA_CONV_FAIL);
 815         case -10:
 816             DDASProcessError(ddas_mem, IDA_CONV_FAIL, "DDASKR", "DDaskrSolve", MSG_CONV_FAILS, tStart);
 817             return (IDA_CONV_FAIL);
 818         case -11:
 819             DDASProcessError(ddas_mem, IDA_RES_FAIL, "DDASKR", "DDaskrSolve", MSG_RES_NONRECOV, tStart);
 820             return (IDA_RES_FAIL);
 821         case -12:
 822             DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrSolve", MSG_IC_FAIL_CONSTR);
 823             return (IDA_ILL_INPUT);
 824         case -33:
 825             DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrSolve", MSG_BAD_INPUT);
 826             return (IDA_ILL_INPUT);
 827         default:
 828             return (IDA_SUCCESS);
 829     }
 830 }
 831
 832 /* =============================
 833  *
 834  *         DDaskrCalcIC
 835  *
 836  * =============================
 837  *
 838  * Computing consistent initial values for the problem.
 839  * This is done by launching the solver with a "stop after consistent initial values computation" flag
 840  * (info[10] and info[13]).
 841  */
 842
 843 int DDaskrCalcIC (void * ddaskr_mem, int icopt, realtype tout1)
 844 {
 845     DDaskrMem ddas_mem = NULL;
 846     double tdist = 0, troundoff = 0, maxnhTemp = 0;
 847
 848     if (ddaskr_mem == NULL)
 849     {
 850         DDASProcessError(NULL, IDA_MEM_NULL, "DDASKR", "DDaskrCalcIC", MSG_NO_MEM);
 851         return (IDA_MEM_NULL);
 852     }
 853     ddas_mem = (DDaskrMem) ddaskr_mem;
 854
 855     /* icopt is a flag to determine whether the DDaskrSetID has been called,
 856        in which case it is now known which components are differential / algebraic.
 857        Here, we only use DDAS_YA_YDP_INIT (DDaskr has been called) */
 858     if (icopt != DDAS_YA_YDP_INIT && icopt != DDAS_Y_INIT)
 859     {
 860         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrCalcIC", MSG_IC_BAD_ICOPT);
 861         return (IDA_ILL_INPUT);
 862     }
 863
 864     /* Checking for valid tout1, not too close to tStart */
 865     tdist = fabs(tout1 - tStart);
 866     troundoff = 2 * UNIT_ROUNDOFF * (fabs(tStart) + fabs(tout1));
 867     if (tdist < troundoff)
 868     {
 869         DDASProcessError(ddas_mem, IDA_ILL_INPUT, "DDASKR", "DDaskrCalcIC", MSG_IC_TOO_CLOSE);
 870         return (IDA_ILL_INPUT);
 871     }
 872
 873     /* info[10] = icopt => Flag on which initial values to compute (differential, algebraic or both)
 874        info[13] = 1 => Do not proceed to integration after calculation */
 875     info[10] = icopt;
 876     if (info[13] == 0)
 877     {
 878         info[13] = 1;
 879     }
 880
 881     /* Giving maxnh a specific value for IC calculation, switching back right after */
 882     if (info[16] == 1)
 883     {
 884         maxnhTemp = iwork[33];
 885         iwork[33] = maxnhIC;
 886     }
 887
 888     C2F(ddaskr) (res, nEq, &tStart, yVec, ypVec, &tout1, info, &relTol, &absTol, &iState, rwork, &lrw, iwork, &liw, rpar, ipar, jacpsol, psol, g_fun, &ng_fun, jroot);
 889
 890     if (info[16] == 1)
 891     {
 892         iwork[33] = maxnhTemp;
 893     }
 894
 895     /* The continuation of the program will not need initial values computation again, unless ReInit */
 896     info[10] = 0;
 897     info[13] = 0;
 898
 899     switch (iState)
 900     {
 901         case 4:
 902             return (IDA_SUCCESS);
 903         default:
 904             DDASProcessError(ddas_mem, IDA_CONV_FAIL, "DDASKR", "DDaskrCalcIC", MSG_IC_CONV_FAILED);
 905             return (IDA_CONV_FAIL);
 906     }
 907 }
 908
 909 /* =============================
 910  *
 911  *     DDaskrGetConsistentIC
 912  *
 913  * =============================
 914  *
 915  * Following on DDasCalcIC, copying yy0 and yp0 (computed consistent values) into the memory space.
 916  */
 917
 918 int DDaskrGetConsistentIC (void * ddaskr_mem, N_Vector yy0, N_Vector yp0)
 919 {
 920     DDaskrMem ddas_mem = NULL;
 921
 922     if (ddaskr_mem == NULL)
 923     {
 924         DDASProcessError(NULL, IDA_MEM_NULL, "DDASKR", "DDaskrGetConsistentIC", MSG_NO_MEM);
 925         return (IDA_MEM_NULL);
 926     }
 927     ddas_mem = (DDaskrMem) ddaskr_mem;
 928
 929     if (yy0 != NULL)
 930     {
 931         NV_DATA_S(yy0) = yVec;
 932     }
 933     if (yp0 != NULL)
 934     {
 935         NV_DATA_S(yp0) = ypVec;
 936     }
 937
 938     return (IDA_SUCCESS);
 939 }
 940
 941 /* =============================
 942  *
 943  *       DDaskrGetRootInfo
 944  *
 945  * =============================
 946  *
 947  * Updates rootsfound[] to the computed roots stocked in jroot[].
 948  */
 949
 950 int DDaskrGetRootInfo (void * ddaskr_mem, int * rootsfound)
 951 {
 952     DDaskrMem ddas_mem = NULL;
 953
 954     if (ddaskr_mem == NULL)
 955     {
 956         DDASProcessError(NULL, IDA_MEM_NULL, "DDASKR", "DDaskrGetRootInfo", MSG_NO_MEM);
 957         return (IDA_MEM_NULL);
 958     }
 959     ddas_mem = (DDaskrMem) ddaskr_mem;
 960
 961     /* Copy jroot to rootsfound */
 962     memcpy(rootsfound, jroot, ng_fun * sizeof(int));
 963
 964     return (IDA_SUCCESS);
 965 }
 966
 967 /* =============================
 968  *
 969  *            DDASFree
 970  *
 971  * =============================
 972  *
 973  * This routine frees the problem memory allocated by DDaskrInit.
 974  */
 975
 976 void DDaskrFree (void ** ddaskr_mem)
 977 {
 978     DDaskrMem ddas_mem = NULL;
 979
 980     if (*ddaskr_mem == NULL)
 981     {
 982         return;
 983     }
 984     ddas_mem = (DDaskrMem) (*ddaskr_mem);
 985
 986     /* Free the inner vectors */
 987     DDASFreeVectors (ddas_mem);
 988
 989     free (*ddaskr_mem);
 990     *ddaskr_mem = NULL;
 991 }
 992
 993 /* =============================
 994  *
 995  *         DDASFreeVectors
 996  *
 997  * =============================
 998  *
 999  * Frees the problem memory space vectors.
1000  */
1001
1002 void DDASFreeVectors (DDaskrMem ddas_mem)
1003 {
1004     /* info, rwork, iwork and jroot have been allocated; free them */
1005     free (info);
1006     free (rwork);
1007     free (iwork);
1008     free (jroot);
1009 }
1010
1011 #define ehfun    ddas_mem->ehfun
1012
1013 /* =============================
1014  *
1015  *     DDaskrSetErrHandlerFn
1016  *
1017  * =============================
1018  *
1019  * Specifies the error handler function.
1020  */
1021
1022 int DDaskrSetErrHandlerFn (void * ddaskr_mem, DDASErrHandlerFn ehFun, void * eh_data)
1023 {
1024     DDaskrMem ddas_mem = NULL;
1025
1026     if (ddaskr_mem == NULL)
1027     {
1028         DDASProcessError(NULL, IDA_MEM_NULL, "DDASKR", "DDaskrSetErrHandlerFn", MSG_NO_MEM);
1029         return(IDA_MEM_NULL);
1030     }
1031
1032     ddas_mem = (DDaskrMem) ddaskr_mem;
1033
1034     ehfun = ehFun;
1035
1036     return(IDA_SUCCESS);
1037 }
1038
1039 /* =============================
1040  *
1041  *         DDASProcessError
1042  *
1043  * =============================
1044  *
1045  * Error handling function.
1046  */
1047
1048 void DDASProcessError (DDaskrMem ddas_mem, int error_code, const char *module, const char *fname, const char *msgfmt, ...)
1049 {
1050     va_list ap;
1051     char msg[256];
1052
1053     /* Initialize the argument pointer variable
1054        (msgfmt is the last required argument to DDASProcessError) */
1055     va_start(ap, msgfmt);
1056
1057     if (ddas_mem == NULL)      /* We write to stderr */
1058     {
1059 #ifndef NO_FPRINTF_OUTPUT
1060         fprintf(stderr, "\n[%s ERROR]  %s\n  ", module, fname);
1061         fprintf(stderr, msgfmt);
1062         fprintf(stderr, "\n\n");
1063 #endif
1064     }
1065     else                     /* We can call ehfun */
1066     {
1067         /* Compose the message */
1068         vsprintf(msg, msgfmt, ap);
1069
1070         /* Call ehfun */
1071         ehfun(error_code, module, fname, msg, NULL);
1072     }
1073
1074     /* Finalize argument processing */
1075     va_end(ap);
1076 }