Differential_equations: daskr subroutines error control
[scilab.git] / scilab / modules / differential_equations / sci_gateway / c / Ex-daskr.c
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include "machine.h"
4 #include "core_math.h"
5
6 extern void C2F(dgesv)(int*, int*, double*, int*, int*, double*, int*, int*);
7 extern double C2F(dlamch)(const char*);
8 typedef void (*resfunc)(double*, double*, double*, double*, int*, double*, int*);
9
10 void C2F(pjac1)(resfunc res, int *ires, int *nequations, double *tOld, double *actual, double *actualP,
11                 double *rewt, double *savr, double *wk, double *h, double *cj, double *wp, int *iwp,
12                 int *ier, double *rpar, int *ipar)
13 {
14     int i = 0;
15     int j = 0;
16     int nrow = 0;
17     double tx = 0;
18     double del = 0;
19     double delinv = 0;
20     double ysave = 0;
21     double ypsave = 0;
22     double* e = NULL;
23
24     int neq = *nequations;
25     double SQuround = sqrt(C2F(dlamch)("P"));
26
27     tx = *tOld;
28
29     e = (double*) calloc(neq, sizeof(double));
30
31     for (i = 0 ; i < neq ; ++i)
32     {
33         del = Max(SQuround * Max(fabs(actual[i]), fabs(*h * actualP[i])), 1. / rewt[i]);
34         del *= (*h * actualP[i] >= 0) ? 1 : -1;
35         ysave   =  actual[i];
36         ypsave  =  actualP[i];
37         actual[i]  += del;
38         actualP[i] += *cj * del;
39         res(&tx, actual, actualP, e, ires, rpar, ipar);
40
41         if (*ires < 0)
42         {
43             *ier = -1;
44             free(e);
45             return;
46         }
47
48         delinv = 1. / del;
49         for (j = 0 ; j < neq ; j++)
50         {
51             wp[nrow + j] = (e[j] - savr[j]) * delinv;
52             /* NaN test */
53             if (ISNAN(wp[nrow + j]))
54             {
55                 *ier = -1;
56                 free (e);
57                 return;
58             }
59             iwp[nrow + j] = i + 1;
60             iwp[nrow + j + neq * neq] = j + 1;
61         }
62         nrow       += neq;
63         actual[i]  =  ysave;
64         actualP[i] =  ypsave;
65     }
66
67     free(e);
68 }
69
70 void C2F(psol1)(int *nequations, double *tOld, double *actual, double *actualP,
71                 double *savr, double *wk, double *cj, double *wght, double *wp,
72                 int *iwp, double *b, double *eplin, int *ier, double *dummy1, int *dummy2)
73 {
74     int i = 0;
75     int nColB = 1;
76     int info = 0;
77     int *ipiv = NULL;
78
79     ipiv = (int *) malloc(*nequations * sizeof(int));
80
81     C2F(dgesv) (nequations, &nColB, wp, nequations, ipiv, b, nequations, &info);
82
83     if (info == 0)
84     {
85         /* NaN test */
86         for (i = 0; i < *nequations; ++i)
87         {
88             if (ISNAN(b[i]))
89             {
90                 *ier = -1;
91                 free (ipiv);
92                 return;
93             }
94         }
95     }
96     else
97     {
98         *ier = -1;
99     }
100
101     free (ipiv);
102 }