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