Line data Source code
1 0 : subroutine pda_lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn,
2 0 : * diag,mode,factor,nprint,info,nfev,fjac,ldfjac,
3 0 : * ipvt,qtf,wa1,wa2,wa3,wa4)
4 : integer m,n,maxfev,mode,nprint,info,nfev,ldfjac
5 : integer ipvt(n)
6 : double precision ftol,xtol,gtol,epsfcn,factor
7 : double precision x(n),fvec(m),diag(n),fjac(ldfjac,n),qtf(n),
8 : * wa1(n),wa2(n),wa3(n),wa4(m)
9 : external fcn
10 : c **********
11 : c
12 : c subroutine pda_lmdif
13 : c
14 : c the purpose of pda_lmdif is to minimize the sum of the squares of
15 : c m nonlinear functions in n variables by a modification of
16 : c the levenberg-marquardt algorithm. the user must provide a
17 : c subroutine which calculates the functions. the jacobian is
18 : c then calculated by a forward-difference approximation.
19 : c
20 : c the subroutine statement is
21 : c
22 : c subroutine pda_lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn,
23 : c diag,mode,factor,nprint,info,nfev,fjac,
24 : c ldfjac,ipvt,qtf,wa1,wa2,wa3,wa4)
25 : c
26 : c where
27 : c
28 : c fcn is the name of the user-supplied subroutine which
29 : c calculates the functions. fcn must be declared
30 : c in an external statement in the user calling
31 : c program, and should be written as follows.
32 : c
33 : c subroutine fcn(m,n,x,fvec,iflag)
34 : c integer m,n,iflag
35 : c double precision x(n),fvec(m)
36 : c ----------
37 : c calculate the functions at x and
38 : c return this vector in fvec.
39 : c ----------
40 : c return
41 : c end
42 : c
43 : c the value of iflag should not be changed by fcn unless
44 : c the user wants to terminate execution of pda_lmdif.
45 : c in this case set iflag to a negative integer.
46 : c
47 : c m is a positive integer input variable set to the number
48 : c of functions.
49 : c
50 : c n is a positive integer input variable set to the number
51 : c of variables. n must not exceed m.
52 : c
53 : c x is an array of length n. on input x must contain
54 : c an initial estimate of the solution vector. on output x
55 : c contains the final estimate of the solution vector.
56 : c
57 : c fvec is an output array of length m which contains
58 : c the functions evaluated at the output x.
59 : c
60 : c ftol is a nonnegative input variable. termination
61 : c occurs when both the actual and predicted relative
62 : c reductions in the sum of squares are at most ftol.
63 : c therefore, ftol measures the relative error desired
64 : c in the sum of squares.
65 : c
66 : c xtol is a nonnegative input variable. termination
67 : c occurs when the relative error between two consecutive
68 : c iterates is at most xtol. therefore, xtol measures the
69 : c relative error desired in the approximate solution.
70 : c
71 : c gtol is a nonnegative input variable. termination
72 : c occurs when the cosine of the angle between fvec and
73 : c any column of the jacobian is at most gtol in absolute
74 : c value. therefore, gtol measures the orthogonality
75 : c desired between the function vector and the columns
76 : c of the jacobian.
77 : c
78 : c maxfev is a positive integer input variable. termination
79 : c occurs when the number of calls to fcn is at least
80 : c maxfev by the end of an iteration.
81 : c
82 : c epsfcn is an input variable used in determining a suitable
83 : c step length for the forward-difference approximation. this
84 : c approximation assumes that the relative errors in the
85 : c functions are of the order of epsfcn. if epsfcn is less
86 : c than the machine precision, it is assumed that the relative
87 : c errors in the functions are of the order of the machine
88 : c precision.
89 : c
90 : c diag is an array of length n. if mode = 1 (see
91 : c below), diag is internally set. if mode = 2, diag
92 : c must contain positive entries that serve as
93 : c multiplicative scale factors for the variables.
94 : c
95 : c mode is an integer input variable. if mode = 1, the
96 : c variables will be scaled internally. if mode = 2,
97 : c the scaling is specified by the input diag. other
98 : c values of mode are equivalent to mode = 1.
99 : c
100 : c factor is a positive input variable used in determining the
101 : c initial step bound. this bound is set to the product of
102 : c factor and the euclidean norm of diag*x if nonzero, or else
103 : c to factor itself. in most cases factor should lie in the
104 : c interval (.1,100.). 100. is a generally recommended value.
105 : c
106 : c nprint is an integer input variable that enables controlled
107 : c printing of iterates if it is positive. in this case,
108 : c fcn is called with iflag = 0 at the beginning of the first
109 : c iteration and every nprint iterations thereafter and
110 : c immediately prior to return, with x and fvec available
111 : c for printing. if nprint is not positive, no special calls
112 : c of fcn with iflag = 0 are made.
113 : c
114 : c info is an integer output variable. if the user has
115 : c terminated execution, info is set to the (negative)
116 : c value of iflag. see description of fcn. otherwise,
117 : c info is set as follows.
118 : c
119 : c info = 0 improper input parameters.
120 : c
121 : c info = 1 both actual and predicted relative reductions
122 : c in the sum of squares are at most ftol.
123 : c
124 : c info = 2 relative error between two consecutive iterates
125 : c is at most xtol.
126 : c
127 : c info = 3 conditions for info = 1 and info = 2 both hold.
128 : c
129 : c info = 4 the cosine of the angle between fvec and any
130 : c column of the jacobian is at most gtol in
131 : c absolute value.
132 : c
133 : c info = 5 number of calls to fcn has reached or
134 : c exceeded maxfev.
135 : c
136 : c info = 6 ftol is too small. no further reduction in
137 : c the sum of squares is possible.
138 : c
139 : c info = 7 xtol is too small. no further improvement in
140 : c the approximate solution x is possible.
141 : c
142 : c info = 8 gtol is too small. fvec is orthogonal to the
143 : c columns of the jacobian to machine precision.
144 : c
145 : c nfev is an integer output variable set to the number of
146 : c calls to fcn.
147 : c
148 : c fjac is an output m by n array. the upper n by n submatrix
149 : c of fjac contains an upper triangular matrix r with
150 : c diagonal elements of nonincreasing magnitude such that
151 : c
152 : c t t t
153 : c p *(jac *jac)*p = r *r,
154 : c
155 : c where p is a permutation matrix and jac is the final
156 : c calculated jacobian. column j of p is column ipvt(j)
157 : c (see below) of the identity matrix. the lower trapezoidal
158 : c part of fjac contains information generated during
159 : c the computation of r.
160 : c
161 : c ldfjac is a positive integer input variable not less than m
162 : c which specifies the leading dimension of the array fjac.
163 : c
164 : c ipvt is an integer output array of length n. ipvt
165 : c defines a permutation matrix p such that jac*p = q*r,
166 : c where jac is the final calculated jacobian, q is
167 : c orthogonal (not stored), and r is upper triangular
168 : c with diagonal elements of nonincreasing magnitude.
169 : c column j of p is column ipvt(j) of the identity matrix.
170 : c
171 : c qtf is an output array of length n which contains
172 : c the first n elements of the vector (q transpose)*fvec.
173 : c
174 : c wa1, wa2, and wa3 are work arrays of length n.
175 : c
176 : c wa4 is a work array of length m.
177 : c
178 : c subprograms called
179 : c
180 : c user-supplied ...... fcn
181 : c
182 : c minpack-supplied ... pda_dpmpar,pda_enorm,pda_fdjac2,pda_lmpar,pda_qrfac
183 : c
184 : c fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod
185 : c
186 : c argonne national laboratory. minpack project. march 1980.
187 : c burton s. garbow, kenneth e. hillstrom, jorge j. more
188 : c
189 : c **********
190 : integer i,iflag,iter,j,l
191 : double precision actred,delta,dirder,epsmch,fnorm,fnorm1,gnorm,
192 : * one,par,pnorm,prered,p1,p5,p25,p75,p0001,ratio,
193 : * sum,temp,temp1,temp2,xnorm,zero
194 : double precision pda_dpmpar,pda_enorm
195 : data one,p1,p5,p25,p75,p0001,zero
196 : * /1.0d0,1.0d-1,5.0d-1,2.5d-1,7.5d-1,1.0d-4,0.0d0/
197 : c
198 : c epsmch is the machine precision.
199 : c
200 0 : epsmch = pda_dpmpar(1)
201 : c
202 0 : info = 0
203 0 : iflag = 0
204 0 : nfev = 0
205 : c
206 : c check the input parameters for errors.
207 : c
208 : if (n .le. 0 .or. m .lt. n .or. ldfjac .lt. m
209 : * .or. ftol .lt. zero .or. xtol .lt. zero .or. gtol .lt. zero
210 0 : * .or. maxfev .le. 0 .or. factor .le. zero) go to 300
211 0 : if (mode .ne. 2) go to 20
212 0 : do 10 j = 1, n
213 0 : if (diag(j) .le. zero) go to 300
214 0 : 10 continue
215 : 20 continue
216 : c
217 : c evaluate the function at the starting point
218 : c and calculate its norm.
219 : c
220 0 : iflag = 1
221 0 : call fcn(m,n,x,fvec,iflag)
222 0 : nfev = 1
223 0 : if (iflag .lt. 0) go to 300
224 0 : fnorm = pda_enorm(m,fvec)
225 : c
226 : c initialize levenberg-marquardt parameter and iteration counter.
227 : c
228 0 : par = zero
229 0 : iter = 1
230 : c
231 : c beginning of the outer loop.
232 : c
233 : 30 continue
234 : c
235 : c calculate the jacobian matrix.
236 : c
237 0 : iflag = 2
238 0 : call pda_fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa4)
239 0 : nfev = nfev + n
240 0 : if (iflag .lt. 0) go to 300
241 : c
242 : c if requested, call fcn to enable printing of iterates.
243 : c
244 0 : if (nprint .le. 0) go to 40
245 0 : iflag = 0
246 0 : if (mod(iter-1,nprint) .eq. 0) call fcn(m,n,x,fvec,iflag)
247 0 : if (iflag .lt. 0) go to 300
248 : 40 continue
249 : c
250 : c compute the qr factorization of the jacobian.
251 : c
252 0 : call pda_qrfac(m,n,fjac,ldfjac,.true.,ipvt,n,wa1,wa2,wa3)
253 : c
254 : c on the first iteration and if mode is 1, scale according
255 : c to the norms of the columns of the initial jacobian.
256 : c
257 0 : if (iter .ne. 1) go to 80
258 0 : if (mode .eq. 2) go to 60
259 0 : do 50 j = 1, n
260 0 : diag(j) = wa2(j)
261 0 : if (wa2(j) .eq. zero) diag(j) = one
262 0 : 50 continue
263 : 60 continue
264 : c
265 : c on the first iteration, calculate the norm of the scaled x
266 : c and initialize the step bound delta.
267 : c
268 0 : do 70 j = 1, n
269 0 : wa3(j) = diag(j)*x(j)
270 0 : 70 continue
271 0 : xnorm = pda_enorm(n,wa3)
272 0 : delta = factor*xnorm
273 0 : if (delta .eq. zero) delta = factor
274 : 80 continue
275 : c
276 : c form (q transpose)*fvec and store the first n components in
277 : c qtf.
278 : c
279 0 : do 90 i = 1, m
280 0 : wa4(i) = fvec(i)
281 0 : 90 continue
282 0 : do 130 j = 1, n
283 0 : if (fjac(j,j) .eq. zero) go to 120
284 0 : sum = zero
285 0 : do 100 i = j, m
286 0 : sum = sum + fjac(i,j)*wa4(i)
287 0 : 100 continue
288 0 : temp = -sum/fjac(j,j)
289 0 : do 110 i = j, m
290 0 : wa4(i) = wa4(i) + fjac(i,j)*temp
291 0 : 110 continue
292 : 120 continue
293 0 : fjac(j,j) = wa1(j)
294 0 : qtf(j) = wa4(j)
295 0 : 130 continue
296 : c
297 : c compute the norm of the scaled gradient.
298 : c
299 0 : gnorm = zero
300 0 : if (fnorm .eq. zero) go to 170
301 0 : do 160 j = 1, n
302 0 : l = ipvt(j)
303 0 : if (wa2(l) .eq. zero) go to 150
304 0 : sum = zero
305 0 : do 140 i = 1, j
306 0 : sum = sum + fjac(i,j)*(qtf(i)/fnorm)
307 0 : 140 continue
308 0 : gnorm = dmax1(gnorm,dabs(sum/wa2(l)))
309 : 150 continue
310 0 : 160 continue
311 : 170 continue
312 : c
313 : c test for convergence of the gradient norm.
314 : c
315 0 : if (gnorm .le. gtol) info = 4
316 0 : if (info .ne. 0) go to 300
317 : c
318 : c rescale if necessary.
319 : c
320 0 : if (mode .eq. 2) go to 190
321 0 : do 180 j = 1, n
322 0 : diag(j) = dmax1(diag(j),wa2(j))
323 0 : 180 continue
324 : 190 continue
325 : c
326 : c beginning of the inner loop.
327 : c
328 : 200 continue
329 : c
330 : c determine the levenberg-marquardt parameter.
331 : c
332 : call pda_lmpar(n,fjac,ldfjac,ipvt,diag,qtf,delta,par,
333 0 : * wa1,wa2,wa3,wa4)
334 : c
335 : c store the direction p and x + p. calculate the norm of p.
336 : c
337 0 : do 210 j = 1, n
338 0 : wa1(j) = -wa1(j)
339 0 : wa2(j) = x(j) + wa1(j)
340 0 : wa3(j) = diag(j)*wa1(j)
341 0 : 210 continue
342 0 : pnorm = pda_enorm(n,wa3)
343 : c
344 : c on the first iteration, adjust the initial step bound.
345 : c
346 0 : if (iter .eq. 1) delta = dmin1(delta,pnorm)
347 : c
348 : c evaluate the function at x + p and calculate its norm.
349 : c
350 0 : iflag = 1
351 0 : call fcn(m,n,wa2,wa4,iflag)
352 0 : nfev = nfev + 1
353 0 : if (iflag .lt. 0) go to 300
354 0 : fnorm1 = pda_enorm(m,wa4)
355 : c
356 : c compute the scaled actual reduction.
357 : c
358 0 : actred = -one
359 0 : if (p1*fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
360 : c
361 : c compute the scaled predicted reduction and
362 : c the scaled directional derivative.
363 : c
364 0 : do 230 j = 1, n
365 0 : wa3(j) = zero
366 0 : l = ipvt(j)
367 0 : temp = wa1(l)
368 0 : do 220 i = 1, j
369 0 : wa3(i) = wa3(i) + fjac(i,j)*temp
370 0 : 220 continue
371 0 : 230 continue
372 0 : temp1 = pda_enorm(n,wa3)/fnorm
373 0 : temp2 = (dsqrt(par)*pnorm)/fnorm
374 0 : prered = temp1**2 + temp2**2/p5
375 0 : dirder = -(temp1**2 + temp2**2)
376 : c
377 : c compute the ratio of the actual to the predicted
378 : c reduction.
379 : c
380 0 : ratio = zero
381 0 : if (prered .ne. zero) ratio = actred/prered
382 : c
383 : c update the step bound.
384 : c
385 0 : if (ratio .gt. p25) go to 240
386 0 : if (actred .ge. zero) temp = p5
387 0 : if (actred .lt. zero)
388 0 : * temp = p5*dirder/(dirder + p5*actred)
389 0 : if (p1*fnorm1 .ge. fnorm .or. temp .lt. p1) temp = p1
390 0 : delta = temp*dmin1(delta,pnorm/p1)
391 0 : par = par/temp
392 0 : go to 260
393 : 240 continue
394 0 : if (par .ne. zero .and. ratio .lt. p75) go to 250
395 0 : delta = pnorm/p5
396 0 : par = p5*par
397 : 250 continue
398 : 260 continue
399 : c
400 : c test for successful iteration.
401 : c
402 0 : if (ratio .lt. p0001) go to 290
403 : c
404 : c successful iteration. update x, fvec, and their norms.
405 : c
406 0 : do 270 j = 1, n
407 0 : x(j) = wa2(j)
408 0 : wa2(j) = diag(j)*x(j)
409 0 : 270 continue
410 0 : do 280 i = 1, m
411 0 : fvec(i) = wa4(i)
412 0 : 280 continue
413 0 : xnorm = pda_enorm(n,wa2)
414 0 : fnorm = fnorm1
415 0 : iter = iter + 1
416 : 290 continue
417 : c
418 : c tests for convergence.
419 : c
420 : if (dabs(actred) .le. ftol .and. prered .le. ftol
421 0 : * .and. p5*ratio .le. one) info = 1
422 0 : if (delta .le. xtol*xnorm) info = 2
423 : if (dabs(actred) .le. ftol .and. prered .le. ftol
424 0 : * .and. p5*ratio .le. one .and. info .eq. 2) info = 3
425 0 : if (info .ne. 0) go to 300
426 : c
427 : c tests for termination and stringent tolerances.
428 : c
429 0 : if (nfev .ge. maxfev) info = 5
430 : if (dabs(actred) .le. epsmch .and. prered .le. epsmch
431 0 : * .and. p5*ratio .le. one) info = 6
432 0 : if (delta .le. epsmch*xnorm) info = 7
433 0 : if (gnorm .le. epsmch) info = 8
434 0 : if (info .ne. 0) go to 300
435 : c
436 : c end of the inner loop. repeat if iteration unsuccessful.
437 : c
438 0 : if (ratio .lt. p0001) go to 200
439 : c
440 : c end of the outer loop.
441 : c
442 0 : go to 30
443 : 300 continue
444 : c
445 : c termination, either normal or user imposed.
446 : c
447 0 : if (iflag .lt. 0) info = iflag
448 0 : iflag = 0
449 0 : if (nprint .gt. 0) call fcn(m,n,x,fvec,iflag)
450 0 : return
451 : c
452 : c last card of subroutine pda_lmdif.
453 : c
454 : end
|
