-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathspartan_sfdjac.c
129 lines (103 loc) · 3.92 KB
/
spartan_sfdjac.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
#include "spartan.h"
/*
J = spartan_sfdjac(userfun,x,J,coloring) Computes a finite-difference
approximation of a sparse Jacobian
Input: userfun A user-defined function of the form
f = userfun(x)
which returns f(x), a vector of length m
x A vector of length n, at which to evaluate J(x)
coloring A coloring of the columns of the Jacobian such
that coloring(j) = k implies that the jth column
of J is colored k.
A coloring is a partition of the columns of
the Jacobian into distinct groups (or colors) such
that no column of the same color has a nonzero entry
in the same row.
All columns of the same color can be estimated via
a single finite-difference. Thus, the number of
of function evaluations required is equal to the number
of colors plus one.
numcolors The number of colors in the coloring
Input/Output:
J
A m x n matrix that contains the sparsity pattern
of the Jacobian. J should be compressed-column format,
storage for the nonzero values should be allocated.
The nonzero values of J are altered.
*/
spartan_sparse *spartan_sfdjac
( spartan_function usrfun,
const double *x,
spartan_sparse *J,
spartan_index *coloring,
spartan_index numcolors
)
{
spartan_index n, j, k, color;
spartan_index *Ji, *Jp;
double *w, *f, *xp, *fp;
double *Jx;
/* For now fix epsilon as the sqrt of
machine precion. The user should be
allowed to pass in a vector of epsilons */
double epsilon = sqrt(DBL_EPSILON);
int err;
n = J->n; Ji = J->i; Jp = J->p; Jx = J->x;
/* allocate work arrays */
w = spartan_malloc(3*n,sizeof(double));
if(!w) return (NULL);
f = w;
xp = w+n;
fp = w+2*n;
/* compute f = F(x) */
err = usrfun(f,x,n);
if (err){ spartan_free(w); return (NULL); };
for (color=0; color<numcolors; color++){
/* Construct xp */
/* xp[j] = x[j] + epsilon if column j is colored color */
/* otherwise xp[j] = x[j] */
for(j=0; j<n; j++){
if (coloring[j] == color)
xp[j] = x[j] + epsilon;
else
xp[j] = x[j];
}
/* Compute fp = userfun(xp) */
err = usrfun(fp,xp,n);
if (err){ spartan_free(w); return (NULL); };
/* fp ~= f + J*(epsilon*p) */
/* fp <- (fp - f)/epsilon */
for(j=0; j<n; j++)
fp[j] = (fp[j] - f[j])/epsilon;
/* place the nonzeros in fp in J */
for(j=0; j<n; j++){
if (coloring[j] == color){
for(k = Jp[j]; k< Jp[j+1]; k++) /* Loop over nonzeros in col j */
Jx[k] = fp[Ji[k]]; /* Set J(i,j) = fp(i) */
}
}
}
return (spartan_done(J,w,NULL,1));
}
/*Global variables to store the nonzero pattern of the Jacobian
and a column coloring */
spartan_sparse *spartan_Jpattern;
spartan_index *spartan_colcoloring;
spartan_index spartan_numcolors;
spartan_function spartan_fdusr;
/* prototype for return to compute Jacobian */
int spartan_fdjacobian(spartan_sparse **J, const double *x, spartan_index n);
spartan_jacobian spartan_initsfd(spartan_function usrfun, spartan_sparse *Jpattern,
spartan_index *colcoloring, spartan_index numcolors){
spartan_fdusr = usrfun;
spartan_Jpattern = Jpattern;
spartan_colcoloring = colcoloring;
spartan_numcolors = numcolors;
return (spartan_fdjacobian);
}
int spartan_fdjacobian(spartan_sparse **J, const double *x, spartan_index n){
spartan_index j,k;
*J = spartan_sfdjac(spartan_fdusr,x,spartan_Jpattern,spartan_colcoloring, spartan_numcolors);
if (!*J) return(-1);
return (0);
}