1 #include <stdio.h>
2 #include <stdlib.h>
3 #include "machine.h"
4 #include "core_math.h"
5 #include "Ex-daskr.h"
6 #include "numericconstants_interface.h"
7
pjac1(resfunc res,int * ires,int * nequations,double * tOld,double * actual,double * actualP,double * rewt,double * savr,double * wk,double * h,double * cj,double * wp,int * iwp,int * ier,double * rpar,int * ipar)8 void pjac1( resfunc res, int *ires, int *nequations, double *tOld, double *actual, double *actualP,
9 double *rewt, double *savr, double *wk, double *h, double *cj, double *wp, int *iwp,
10 int *ier, double *rpar, int *ipar)
11 {
12 int i = 0;
13 int j = 0;
14 int info = 0;
15 int nrow = 0;
16 double tx = 0;
17 double del = 0;
18 double delinv = 0;
19 double ysave = 0;
20 double ypsave = 0;
21 double * e = NULL;
22
23 int neq = *nequations;
24 double SQuround = sqrt(nc_eps_machine());
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 del = (actual[i] + del) - actual[i];
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
53 /* NaN test */
54 if (ISNAN(wp[nrow + j]))
55 {
56 *ier = -1;
57 free(e);
58 return;
59 }
60 }
61 nrow += neq;
62 actual[i] = ysave;
63 actualP[i] = ypsave;
64 }
65
66 /* Proceed to LU factorization of P. */
67 C2F(dgefa) (wp, nequations, nequations, iwp, &info);
68 if (info != 0)
69 {
70 *ier = -1;
71 }
72
73 free(e);
74 }
75
psol1(int * nequations,double * tOld,double * actual,double * actualP,double * savr,double * wk,double * cj,double * wght,double * wp,int * iwp,double * b,double * eplin,int * ier,double * dummy1,int * dummy2)76 void psol1( int *nequations, double *tOld, double *actual, double *actualP,
77 double *savr, double *wk, double *cj, double *wght, double *wp,
78 int *iwp, double *b, double *eplin, int *ier, double *dummy1, int *dummy2)
79 {
80 int i = 0, job = 0;
81
82 C2F(dgesl) (wp, nequations, nequations, iwp, b, &job);
83
84 /* NaN test */
85 for (i = 0; i < *nequations; ++i)
86 {
87 if (ISNAN(b[i]))
88 {
89 /* Indicate a recoverable error, meaning that the step will be retried with the same step size,
90 but with a call to 'jacpsol' to update necessary data, unless the Jacobian data is current,
91 if (b[i] - b[i] != 0) in which case the step will be retried with a smaller step size. */
92 *ier = -1;
93 return;
94 }
95 }
96 }
97