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