f1c9c6a1a63349e49a6cda0a9b7f5b0d64b7aaa1
[scilab.git] / 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 }