01ae0fd001a3b65cba2f81d7a5e8e4d6e8e722db
[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 int C2F(dgesv)(int*, int*, double*, int*, int*, double*, int*, int*);
7 extern double C2F(dlamch)(char*);
8 typedef void (*resfunc)(double*, double*, double*, double*, int*, double*, int*);
9
10 int 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     int neq = *nequations;
24     double SQuround = sqrt(C2F(dlamch)("P"));
25
26     tx = *tOld;
27
28     e = (double*)calloc(neq, sizeof(double));
29
30     for (i = 0 ; i < neq ; ++i)
31     {
32         del = Max(SQuround * Max(fabs(actual[i]), fabs(*h * actualP[i])), 1. / rewt[i]);
33         del *= (*h * actualP[i] >= 0) ? 1 : -1;
34         ysave   =  actual[i];
35         ypsave  =  actualP[i];
36         actual[i]  += del;
37         actualP[i] += *cj * del;
38         res(&tx, actual, actualP, e, ires, rpar, ipar);
39
40         if (ires < 0)
41         {
42             return 0;
43         }
44
45         delinv = 1. / del;
46         for (j = 0 ; j < neq ; j++)
47         {
48             wp[nrow + j] = (e[j] - savr[j]) * delinv;
49             iwp[nrow + j] = i + 1;
50             iwp[nrow + j + neq * neq] = j + 1;
51         }
52         nrow       += neq;
53         actual[i]  =  ysave;
54         actualP[i] =  ypsave;
55     }
56
57     free(e);
58
59     return (*ier);
60 }
61
62 int C2F(psol1)(int *nequations, double *tOld, double *actual, double *actualP,
63                double *savr, double *wk, double *cj, double *wght, double *wp,
64                int *iwp, double *b, double *eplin, int *ier, double *dummy1, int *dummy2)
65 {
66     int nColB = 1;
67     int info = 0;
68     int *ipiv = NULL;
69
70     ipiv = (int *) malloc(*nequations * sizeof(int));
71
72     C2F(dgesv) (nequations, &nColB, wp, nequations, ipiv, b, nequations, &info);
73
74     free(ipiv);
75     return 0;
76 }
77