1 #include <math.h>
2 #include "cminpack.h"
3 #include "ssq.h"
4 #define real __cminpack_real__
5
ssqjac(int m,int n,const real * x,real * fjac,int ldfjac,int nprob)6 void ssqjac(int m, int n, const real *x, real *fjac, int ldfjac, int nprob)
7 {
8 /* Initialized data */
9
10 const real v[11] = { 4.,2.,1.,.5,.25,.167,.125,.1,.0833,.0714, .0625 };
11
12 /* System generated locals */
13 real d__1;
14
15 /* Local variables */
16 int i, j, k;
17 real s2, dx, ti;
18 int mm1, nm1;
19 real div, tpi, tmp1, tmp2, tmp3, tmp4, prod, temp;
20
21 /* ********** */
22
23 /* subroutine ssqjac */
24
25 /* this subroutine defines the jacobian matrices of eighteen */
26 /* nonlinear least squares problems. the problem dimensions are */
27 /* as described in the prologue comments of ssqfcn. */
28
29 /* the subroutine statement is */
30
31 /* subroutine ssqjac(m,n,x,fjac,ldfjac,nprob) */
32
33 /* where */
34
35 /* m and n are positive int input variables. n must not */
36 /* exceed m. */
37
38 /* x is an input array of length n. */
39
40 /* fjac is an m by n output array which contains the jacobian */
41 /* matrix of the nprob function evaluated at x. */
42
43 /* ldfjac is a positive int input variable not less than m */
44 /* which specifies the leading dimension of the array fjac. */
45
46 /* nprob is a positive int variable which defines the */
47 /* number of the problem. nprob must not exceed 18. */
48
49 /* subprograms called */
50
51 /* fortran-supplied ... datan,dcos,dexp,dsin,dsqrt */
52
53 /* argonne national laboratory. minpack project. march 1980. */
54 /* burton s. garbow, kenneth e. hillstrom, jorge j. more */
55
56 /* ********** */
57 /* Parameter adjustments */
58 --x;
59 fjac -= 1 + ldfjac;
60
61 /* Function Body */
62
63 /* jacobian routine selector. */
64
65 switch (nprob) {
66
67 /* linear function - full rank. */
68
69 case 1:
70 temp = 2. / (real) m;
71 for (j = 1; j <= n; ++j) {
72 for (i = 1; i <= m; ++i) {
73 fjac[i + j * ldfjac] = -temp;
74 }
75 fjac[j + j * ldfjac] += 1.;
76 }
77 break;
78
79 /* linear function - rank 1. */
80
81 case 2:
82 for (j = 1; j <= n; ++j) {
83 for (i = 1; i <= m; ++i) {
84 fjac[i + j * ldfjac] = (real) i * (real) j;
85 }
86 }
87 break;
88
89 /* linear function - rank 1 with zero columns and rows. */
90
91 case 3:
92 for (j = 1; j <= n; ++j) {
93 for (i = 1; i <= m; ++i) {
94 fjac[i + j * ldfjac] = 0.;
95 }
96 }
97 #if 0
98 if (nm1 < 2) {
99 goto L120;
100 }
101 for (j = 2; j <= nm1; ++j) {
102 i__2 = mm1;
103 for (i__ = 2; i__ <= mm1; ++i__) {
104 i__3 = i__ - 1;
105 fjac[i__ + j * fjac_dim1] = (doublereal) i__3 * (doublereal) j;
106 /* L100: */
107 }
108 /* L110: */
109 }
110 L120:
111 #else
112 nm1 = n - 1;
113 mm1 = m - 1;
114 if (nm1 >= 2) {
115 for (j = 2; j <= nm1; ++j) {
116 for (i = 2; i <= mm1; ++i) {
117 fjac[i + j * ldfjac] = (real) (i - 1) * (real) j;
118 }
119 }
120 }
121 #endif
122 break;
123
124 /* rosenbrock function. */
125
126 case 4:
127 fjac[1 + 1 * ldfjac] = -20. * x[1];
128 fjac[1 + 2 * ldfjac] = 10.;
129 fjac[2 + 1 * ldfjac] = -1.;
130 fjac[2 + 2 * ldfjac] = 0.;
131 break;
132
133 /* helical valley function. */
134
135 case 5:
136 tpi = 8. * atan(1.);
137 temp = x[1] * x[1] + x[2] * x[2];
138 tmp1 = tpi * temp;
139 tmp2 = sqrt(temp);
140 fjac[1 + 1 * ldfjac] = 100. * x[2] / tmp1;
141 fjac[1 + 2 * ldfjac] = -100. * x[1] / tmp1;
142 fjac[1 + 3 * ldfjac] = 10.;
143 fjac[2 + 1 * ldfjac] = 10. * x[1] / tmp2;
144 fjac[2 + 2 * ldfjac] = 10. * x[2] / tmp2;
145 fjac[2 + 3 * ldfjac] = 0.;
146 fjac[3 + 1 * ldfjac] = 0.;
147 fjac[3 + 2 * ldfjac] = 0.;
148 fjac[3 + 3 * ldfjac] = 1.;
149 break;
150
151 /* powell singular function. */
152
153 case 6:
154 for (j = 1; j <= 4; ++j) {
155 for (i = 1; i <= 4; ++i) {
156 fjac[i + j * ldfjac] = 0.;
157 }
158 }
159 fjac[1 + 1 * ldfjac] = 1.;
160 fjac[1 + 2 * ldfjac] = 10.;
161 fjac[2 + 3 * ldfjac] = sqrt(5.);
162 fjac[2 + 4 * ldfjac] = -fjac[2 + 3 * ldfjac];
163 fjac[3 + 2 * ldfjac] = 2. * (x[2] - 2. * x[3]);
164 fjac[3 + 3 * ldfjac] = -2. * fjac[3 + 2 * ldfjac];
165 fjac[4 + 1 * ldfjac] = 2. * sqrt(10.) * (x[1] - x[4]);
166 fjac[4 + 4 * ldfjac] = -fjac[4 + 1 * ldfjac];
167 break;
168
169 /* freudenstein and roth function. */
170
171 case 7:
172 fjac[1 + 1 * ldfjac] = 1.;
173 fjac[1 + 2 * ldfjac] = x[2] * (10. - 3. * x[2]) - 2.;
174 fjac[2 + 1 * ldfjac] = 1.;
175 fjac[2 + 2 * ldfjac] = x[2] * (2. + 3. * x[2]) - 14.;
176 break;
177
178 /* bard function. */
179
180 case 8:
181 for (i = 1; i <= 15; ++i) {
182 tmp1 = (real) i;
183 tmp2 = (real) (16 - i);
184 tmp3 = tmp1;
185 if (i > 8) {
186 tmp3 = tmp2;
187 }
188 /* Computing 2nd power */
189 d__1 = x[2] * tmp2 + x[3] * tmp3;
190 tmp4 = d__1 * d__1;
191 fjac[i + 1 * ldfjac] = -1.;
192 fjac[i + 2 * ldfjac] = tmp1 * tmp2 / tmp4;
193 fjac[i + 3 * ldfjac] = tmp1 * tmp3 / tmp4;
194 }
195 break;
196
197 /* kowalik and osborne function. */
198
199 case 9:
200 for (i = 1; i <= 11; ++i) {
201 tmp1 = v[i - 1] * (v[i - 1] + x[2]);
202 tmp2 = v[i - 1] * (v[i - 1] + x[3]) + x[4];
203 fjac[i + 1 * ldfjac] = -tmp1 / tmp2;
204 fjac[i + 2 * ldfjac] = -v[i - 1] * x[1] / tmp2;
205 fjac[i + 3 * ldfjac] = fjac[i + 1 * ldfjac] * fjac[i + 2 * ldfjac];
206 fjac[i + 4 * ldfjac] = fjac[i + 3 * ldfjac] / v[i - 1];
207 }
208 break;
209
210 /* meyer function. */
211
212 case 10:
213 for (i = 1; i <= 16; ++i) {
214 temp = 5. * (real) i + 45. + x[3];
215 tmp1 = x[2] / temp;
216 tmp2 = exp(tmp1);
217 fjac[i + 1 * ldfjac] = tmp2;
218 fjac[i + 2 * ldfjac] = x[1] * tmp2 / temp;
219 fjac[i + 3 * ldfjac] = -tmp1 * fjac[i + 2 * ldfjac];
220 }
221 break;
222
223 /* watson function. */
224
225 case 11:
226 for (i = 1; i <= 29; ++i) {
227 div = (real) i / 29.;
228 s2 = 0.;
229 dx = 1.;
230 for (j = 1; j <= n; ++j) {
231 s2 += dx * x[j];
232 dx = div * dx;
233 }
234 temp = 2. * div * s2;
235 dx = 1. / div;
236 for (j = 1; j <= n; ++j) {
237 fjac[i + j * ldfjac] = dx * ((real) (j - 1) - temp);
238 dx = div * dx;
239 }
240 }
241 for (j = 1; j <= n; ++j) {
242 for (i = 30; i <= 31; ++i) {
243 fjac[i + j * ldfjac] = 0.;
244 }
245 }
246 fjac[30 + 1 * ldfjac] = 1.;
247 fjac[31 + 1 * ldfjac] = -2. * x[1];
248 fjac[31 + 2 * ldfjac] = 1.;
249 break;
250
251 /* box 3-dimensional function. */
252
253 case 12:
254 for (i = 1; i <= m; ++i) {
255 temp = (real) i;
256 tmp1 = temp / 10.;
257 fjac[i + 1 * ldfjac] = -tmp1 * exp(-tmp1 * x[1]);
258 fjac[i + 2 * ldfjac] = tmp1 * exp(-tmp1 * x[2]);
259 fjac[i + 3 * ldfjac] = exp(-temp) - exp(-tmp1);
260 }
261 break;
262
263 /* jennrich and sampson function. */
264
265 case 13:
266 for (i = 1; i <= m; ++i) {
267 temp = (real) i;
268 fjac[i + 1 * ldfjac] = -temp * exp(temp * x[1]);
269 fjac[i + 2 * ldfjac] = -temp * exp(temp * x[2]);
270 }
271 break;
272
273 /* brown and dennis function. */
274
275 case 14:
276 for (i = 1; i <= m; ++i) {
277 temp = (real) i / 5.;
278 ti = sin(temp);
279 tmp1 = x[1] + temp * x[2] - exp(temp);
280 tmp2 = x[3] + ti * x[4] - cos(temp);
281 fjac[i + 1 * ldfjac] = 2. * tmp1;
282 fjac[i + 2 * ldfjac] = temp * fjac[i + 1 * ldfjac];
283 fjac[i + 3 * ldfjac] = 2. * tmp2;
284 fjac[i + 4 * ldfjac] = ti * fjac[i + 3 * ldfjac];
285 }
286 break;
287
288 /* chebyquad function. */
289
290 case 15:
291 dx = 1. / (real) (n);
292 for (j = 1; j <= n; ++j) {
293 tmp1 = 1.;
294 tmp2 = 2. * x[j] - 1.;
295 temp = 2. * tmp2;
296 tmp3 = 0.;
297 tmp4 = 2.;
298 for (i = 1; i <= m; ++i) {
299 fjac[i + j * ldfjac] = dx * tmp4;
300 ti = 4. * tmp2 + temp * tmp4 - tmp3;
301 tmp3 = tmp4;
302 tmp4 = ti;
303 ti = temp * tmp2 - tmp1;
304 tmp1 = tmp2;
305 tmp2 = ti;
306 }
307 }
308 break;
309
310 /* brown almost-linear function. */
311
312 case 16:
313 prod = 1.;
314 for (j = 1; j <= n; ++j) {
315 prod = x[j] * prod;
316 for (i = 1; i <= n; ++i) {
317 fjac[i + j * ldfjac] = 1.;
318 }
319 fjac[j + j * ldfjac] = 2.;
320 }
321 for (j = 1; j <= n; ++j) {
322 temp = x[j];
323 if (temp == 0.) {
324 temp = 1.;
325 prod = 1.;
326 for (k = 1; k <= n; ++k) {
327 if (k != j) {
328 prod = x[k] * prod;
329 }
330 }
331 }
332 fjac[n + j * ldfjac] = prod / temp;
333 }
334 break;
335
336 /* osborne 1 function. */
337
338 case 17:
339 for (i = 1; i <= 33; ++i) {
340 temp = 10. * (real) (i - 1);
341 tmp1 = exp(-x[4] * temp);
342 tmp2 = exp(-x[5] * temp);
343 fjac[i + 1 * ldfjac] = -1.;
344 fjac[i + 2 * ldfjac] = -tmp1;
345 fjac[i + 3 * ldfjac] = -tmp2;
346 fjac[i + 4 * ldfjac] = temp * x[2] * tmp1;
347 fjac[i + 5 * ldfjac] = temp * x[3] * tmp2;
348 }
349 break;
350
351 /* osborne 2 function. */
352
353 case 18:
354 for (i = 1; i <= 65; ++i) {
355 temp = (real) (i - 1) / 10.;
356 tmp1 = exp(-x[5] * temp);
357 /* Computing 2nd power */
358 d__1 = temp - x[9];
359 tmp2 = exp(-x[6] * (d__1 * d__1));
360 /* Computing 2nd power */
361 d__1 = temp - x[10];
362 tmp3 = exp(-x[7] * (d__1 * d__1));
363 /* Computing 2nd power */
364 d__1 = temp - x[11];
365 tmp4 = exp(-x[8] * (d__1 * d__1));
366 fjac[i + 1 * ldfjac] = -tmp1;
367 fjac[i + 2 * ldfjac] = -tmp2;
368 fjac[i + 3 * ldfjac] = -tmp3;
369 fjac[i + 4 * ldfjac] = -tmp4;
370 fjac[i + 5 * ldfjac] = temp * x[1] * tmp1;
371 /* Computing 2nd power */
372 d__1 = temp - x[9];
373 fjac[i + 6 * ldfjac] = x[2] * (d__1 * d__1) * tmp2;
374 /* Computing 2nd power */
375 d__1 = temp - x[10];
376 fjac[i + 7 * ldfjac] = x[3] * (d__1 * d__1) * tmp3;
377 /* Computing 2nd power */
378 d__1 = temp - x[11];
379 fjac[i + 8 * ldfjac] = x[4] * (d__1 * d__1) * tmp4;
380 fjac[i + 9 * ldfjac] = -2. * x[2] * x[6] * (temp - x[9]) * tmp2;
381 fjac[i + 10 * ldfjac] = -2. * x[3] * x[7] * (temp - x[10]) * tmp3;
382 fjac[i + 11 * ldfjac] = -2. * x[4] * x[8] * (temp - x[11]) * tmp4;
383 }
384 break;
385 }
386
387 /* last card of subroutine ssqjac. */
388
389 } /* ssqjac_ */
390
391