MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
levenberg_marquardt.f90
1 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2 ! MINPACK routines which are used by both LMDIF & LMDER
3 ! 25 October 2001:
4 ! Changed INTENT of iflag in several places to IN OUT.
5 ! Changed INTENT of fvec to IN OUT in user routine FCN.
6 ! Removed arguments diag and qtv from LMDIF & LMDER.
7 ! Replaced several DO loops with array operations.
8 ! amiller @ bigpond.net.au
9 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
10 
11 MODULE levenberg_marquardt
12 
13  IMPLICIT NONE
14  INTEGER, PARAMETER :: dp = selected_real_kind(15, 307)
15 
16  PRIVATE
17  PUBLIC :: dp, lmdif1, lmdif, lmder1, lmder, enorm
18 
19 CONTAINS
20 
21 
22  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
23  ! subroutine lmdif1
24  !
25  ! The purpose of lmdif1 is to minimize the sum of the squares of m nonlinear
26  ! functions in n variables by a modification of the Levenberg-Marquardt
27  ! algorithm. This is done by using the more general least-squares
28  ! solver lmdif. The user must provide a subroutine which calculates the
29  ! functions. The jacobian is then calculated by a forward-difference
30  ! approximation.
31  !
32  ! the subroutine statement is
33  !
34  ! subroutine lmdif1(fcn, m, n, x, fvec, tol, epsfcn, info, iwa)
35  !
36  ! where
37  !
38  ! fcn is the name of the user-supplied subroutine which calculates
39  ! the functions. fcn must be declared in an external statement in the
40  ! user calling program, and should be written as follows.
41 
42  ! subroutine fcn(m, n, x, fvec, iflag)
43  ! integer m, n, iflag
44  ! REAL x(n), fvec(m)
45  ! ----------
46  ! calculate the functions at x and return this vector in fvec.
47  ! ----------
48  ! return
49  ! end
50  !
51  ! the value of iflag should not be changed by fcn unless
52  ! the user wants to terminate execution of lmdif1.
53  ! In this case set iflag to a negative integer.
54  !
55  ! m is a positive integer input variable set to the number of functions.
56  !
57  ! n is a positive integer input variable set to the number of variables.
58  ! n must not exceed m.
59  !
60  ! x is an array of length n. On input x must contain an initial estimate
61  ! of the solution vector. On output x contains the final estimate of
62  ! the solution vector.
63  !
64  ! fvec is an output array of length m which contains
65  ! the functions evaluated at the output x.
66  !
67  ! tol is a nonnegative input variable. Termination occurs when the
68  ! algorithm estimates either that the relative error in the sum of
69  ! squares is at most tol or that the relative error between x and the
70  ! solution is at most tol.
71  !
72  ! info is an integer output variable. If the user has terminated execution,
73  ! info is set to the (negative) value of iflag. See description of fcn.
74  ! Otherwise, info is set as follows.
75  !
76  ! info = 0 improper input parameters.
77  !
78  ! info = 1 algorithm estimates that the relative error
79  ! in the sum of squares is at most tol.
80  !
81  ! info = 2 algorithm estimates that the relative error
82  ! between x and the solution is at most tol.
83  !
84  ! info = 3 conditions for info = 1 and info = 2 both hold.
85  !
86  ! info = 4 fvec is orthogonal to the columns of the
87  ! jacobian to machine precision.
88  !
89  ! info = 5 number of calls to fcn has reached or exceeded 200*(n+1).
90  !
91  ! info = 6 tol is too small. no further reduction in
92  ! the sum of squares is possible.
93  !
94  ! info = 7 tol is too small. No further improvement in
95  ! the approximate solution x is possible.
96  !
97  ! iwa is an integer work array of length n.
98  !
99  ! wa is a work array of length lwa.
100  !
101  ! lwa is a positive integer input variable not less than m*n+5*n+m.
102  !
103  ! subprograms called
104  !
105  ! user-supplied ...... fcn
106  !
107  ! minpack-supplied ... lmdif
108  !
109  ! argonne national laboratory. minpack project. march 1980.
110  ! burton s. garbow, kenneth e. hillstrom, jorge j. more
111  !
112  ! Code converted using TO_F90 by Alan Miller
113  ! Date: 1999-12-11 Time: 00:51:44
114  ! N.B. Arguments WA & LWA have been removed.
115  !
116  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
117 
118  SUBROUTINE lmdif1(fcn, m, n, x, fvec, tol, epsfcn, factor, info, iwa)
119 
120  INTEGER, INTENT(IN) :: m
121  INTEGER, INTENT(IN) :: n
122  REAL, INTENT(IN OUT) :: x(:)
123  REAL, INTENT(OUT) :: fvec(:)
124  REAL, INTENT(IN) :: tol
125  REAL, INTENT(IN OUT) :: epsfcn
126  REAL, INTENT(IN) :: factor
127  INTEGER, INTENT(OUT) :: info
128  INTEGER, INTENT(OUT) :: iwa(:)
129 
130  ! EXTERNAL fcn
131 
132  INTERFACE
133  SUBROUTINE fcn(m, n, x, fvec, iflag)
134  IMPLICIT NONE
135  ! INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 307)
136  INTEGER, INTENT(IN) :: m, n
137  REAL, INTENT(IN) :: x(:)
138  REAL, INTENT(IN OUT) :: fvec(:)
139  INTEGER, INTENT(IN OUT) :: iflag
140  END SUBROUTINE fcn
141  END INTERFACE
142 
143  INTEGER :: maxfev, mode, nfev, nprint
144  REAL :: ftol, gtol, xtol, fjac(m,n)
145  !REAL, PARAMETER :: factor = 100._dp, zero = 0.0_dp
146  REAL, PARAMETER :: zero = 0.0
147 
148 
149  info = 0
150 
151  ! check the input parameters for errors.
152 
153  IF (n <= 0 .OR. m < n .OR. tol < zero) RETURN
154 
155  ! call lmdif.
156 
157  maxfev = 200*(n + 1)
158  ftol = tol
159  xtol = tol
160  gtol = zero
161  !epsfcn = zero
162  mode = 1
163  nprint = 0
164  CALL lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, &
165  mode, factor, nprint, info, nfev, fjac, iwa)
166  IF (info == 8) info = 4
167 
168  END SUBROUTINE lmdif1
169 
170 
171 
172  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
173  ! subroutine lmdif
174 
175  ! The purpose of lmdif is to minimize the sum of the squares of m nonlinear
176  ! functions in n variables by a modification of the Levenberg-Marquardt
177  ! algorithm. The user must provide a subroutine which calculates the
178  ! functions. The jacobian is then calculated by a forward-difference
179  ! approximation.
180 
181  ! the subroutine statement is
182 
183  ! subroutine lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn,
184  ! diag, mode, factor, nprint, info, nfev, fjac,
185  ! ldfjac, ipvt, qtf, wa1, wa2, wa3, wa4)
186 
187  ! N.B. 7 of these arguments have been removed in this version.
188 
189  ! where
190 
191  ! fcn is the name of the user-supplied subroutine which calculates the
192  ! functions. fcn must be declared in an external statement in the user
193  ! calling program, and should be written as follows.
194 
195  ! subroutine fcn(m, n, x, fvec, iflag)
196  ! integer m, n, iflag
197  ! REAL x(:), fvec(m)
198  ! ----------
199  ! calculate the functions at x and return this vector in fvec.
200  ! ----------
201  ! return
202  ! end
203 
204  ! the value of iflag should not be changed by fcn unless
205  ! the user wants to terminate execution of lmdif.
206  ! in this case set iflag to a negative integer.
207 
208  ! m is a positive integer input variable set to the number of functions.
209 
210  ! n is a positive integer input variable set to the number of variables.
211  ! n must not exceed m.
212 
213  ! x is an array of length n. On input x must contain an initial estimate
214  ! of the solution vector. On output x contains the final estimate of the
215  ! solution vector.
216 
217  ! fvec is an output array of length m which contains
218  ! the functions evaluated at the output x.
219 
220  ! ftol is a nonnegative input variable. Termination occurs when both the
221  ! actual and predicted relative reductions in the sum of squares are at
222  ! most ftol. Therefore, ftol measures the relative error desired
223  ! in the sum of squares.
224 
225  ! xtol is a nonnegative input variable. Termination occurs when the
226  ! relative error between two consecutive iterates is at most xtol.
227  ! Therefore, xtol measures the relative error desired in the approximate
228  ! solution.
229 
230  ! gtol is a nonnegative input variable. Termination occurs when the cosine
231  ! of the angle between fvec and any column of the jacobian is at most
232  ! gtol in absolute value. Therefore, gtol measures the orthogonality
233  ! desired between the function vector and the columns of the jacobian.
234 
235  ! maxfev is a positive integer input variable. Termination occurs when the
236  ! number of calls to fcn is at least maxfev by the end of an iteration.
237 
238  ! epsfcn is an input variable used in determining a suitable step length
239  ! for the forward-difference approximation. This approximation assumes
240  ! that the relative errors in the functions are of the order of epsfcn.
241  ! If epsfcn is less than the machine precision, it is assumed that the
242  ! relative errors in the functions are of the order of the machine
243  ! precision.
244 
245  ! diag is an array of length n. If mode = 1 (see below), diag is
246  ! internally set. If mode = 2, diag must contain positive entries that
247  ! serve as multiplicative scale factors for the variables.
248 
249  ! mode is an integer input variable. If mode = 1, the variables will be
250  ! scaled internally. If mode = 2, the scaling is specified by the input
251  ! diag. other values of mode are equivalent to mode = 1.
252 
253  ! factor is a positive input variable used in determining the initial step
254  ! bound. This bound is set to the product of factor and the euclidean
255  ! norm of diag*x if nonzero, or else to factor itself. In most cases
256  ! factor should lie in the interval (.1,100.). 100. is a generally
257  ! recommended value.
258 
259  ! nprint is an integer input variable that enables controlled printing of
260  ! iterates if it is positive. In this case, fcn is called with iflag = 0
261  ! at the beginning of the first iteration and every nprint iterations
262  ! thereafter and immediately prior to return, with x and fvec available
263  ! for printing. If nprint is not positive, no special calls
264  ! of fcn with iflag = 0 are made.
265 
266  ! info is an integer output variable. If the user has terminated
267  ! execution, info is set to the (negative) value of iflag.
268  ! See description of fcn. Otherwise, info is set as follows.
269 
270  ! info = 0 improper input parameters.
271 
272  ! info = 1 both actual and predicted relative reductions
273  ! in the sum of squares are at most ftol.
274 
275  ! info = 2 relative error between two consecutive iterates <= xtol.
276 
277  ! info = 3 conditions for info = 1 and info = 2 both hold.
278 
279  ! info = 4 the cosine of the angle between fvec and any column of
280  ! the Jacobian is at most gtol in absolute value.
281 
282  ! info = 5 number of calls to fcn has reached or exceeded maxfev.
283 
284  ! info = 6 ftol is too small. no further reduction in
285  ! the sum of squares is possible.
286 
287  ! info = 7 xtol is too small. no further improvement in
288  ! the approximate solution x is possible.
289 
290  ! info = 8 gtol is too small. fvec is orthogonal to the
291  ! columns of the jacobian to machine precision.
292 
293  ! nfev is an integer output variable set to the number of calls to fcn.
294 
295  ! fjac is an output m by n array. the upper n by n submatrix
296  ! of fjac contains an upper triangular matrix r with
297  ! diagonal elements of nonincreasing magnitude such that
298 
299  ! t t t
300  ! p *(jac *jac)*p = r *r,
301 
302  ! where p is a permutation matrix and jac is the final calculated
303  ! Jacobian. Column j of p is column ipvt(j) (see below) of the
304  ! identity matrix. the lower trapezoidal part of fjac contains
305  ! information generated during the computation of r.
306 
307  ! ldfjac is a positive integer input variable not less than m
308  ! which specifies the leading dimension of the array fjac.
309 
310  ! ipvt is an integer output array of length n. ipvt defines a permutation
311  ! matrix p such that jac*p = q*r, where jac is the final calculated
312  ! jacobian, q is orthogonal (not stored), and r is upper triangular
313  ! with diagonal elements of nonincreasing magnitude.
314  ! Column j of p is column ipvt(j) of the identity matrix.
315 
316  ! qtf is an output array of length n which contains
317  ! the first n elements of the vector (q transpose)*fvec.
318 
319  ! wa1, wa2, and wa3 are work arrays of length n.
320 
321  ! wa4 is a work array of length m.
322 
323  ! subprograms called
324 
325  ! user-supplied ...... fcn
326 
327  ! minpack-supplied ... dpmpar,enorm,fdjac2,lmpar,qrfac
328 
329  ! fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod
330  !
331  ! argonne national laboratory. minpack project. march 1980.
332  ! burton s. garbow, kenneth e. hillstrom, jorge j. more
333  !
334  ! N.B. Arguments LDFJAC, DIAG, QTF, WA1, WA2, WA3 & WA4 have been removed.
335  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
336 
337  SUBROUTINE lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, &
338  mode, factor, nprint, info, nfev, fjac, ipvt)
339 
340 
341  INTEGER, INTENT(IN) :: m
342  INTEGER, INTENT(IN) :: n
343  REAL, INTENT(IN OUT) :: x(:)
344  REAL, INTENT(OUT) :: fvec(:)
345  REAL, INTENT(IN) :: ftol
346  REAL, INTENT(IN) :: xtol
347  REAL, INTENT(IN OUT) :: gtol
348  INTEGER, INTENT(IN OUT) :: maxfev
349  REAL, INTENT(IN OUT) :: epsfcn
350  INTEGER, INTENT(IN) :: mode
351  REAL, INTENT(IN) :: factor
352  INTEGER, INTENT(IN) :: nprint
353  INTEGER, INTENT(OUT) :: info
354  INTEGER, INTENT(OUT) :: nfev
355  REAL, INTENT(OUT) :: fjac(:,:) ! fjac(ldfjac,n)
356  INTEGER, INTENT(OUT) :: ipvt(:)
357 
358  ! EXTERNAL fcn
359 
360  INTERFACE
361  SUBROUTINE fcn(m, n, x, fvec, iflag)
362  IMPLICIT NONE
363  ! INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 307)
364  INTEGER, INTENT(IN) :: m, n
365  REAL, INTENT(IN) :: x(:)
366  REAL, INTENT(IN OUT) :: fvec(:)
367  INTEGER, INTENT(IN OUT) :: iflag
368  END SUBROUTINE fcn
369  END INTERFACE
370 
371  INTEGER :: i, iflag, iter, j, l
372  REAL :: actred, delta, dirder, epsmch, fnorm, fnorm1, gnorm, &
373  par, pnorm, prered, ratio, sum, temp, temp1, temp2, xnorm
374  REAL :: diag(n), qtf(n), wa1(n), wa2(n), wa3(n), wa4(m)
375  REAL, PARAMETER :: one = 1.0, p1 = 0.1, p5 = 0.5, &
376  p25 = 0.25, p75 = 0.75, p0001 = 0.0001, &
377  zero = 0.0
378 
379  ! epsmch is the machine precision.
380 
381  epsmch = epsilon(zero)
382 
383  info = 0
384  iflag = 0
385  nfev = 0
386 
387  ! check the input parameters for errors.
388 
389  IF (n <= 0 .OR. m < n .OR. ftol < zero .OR. xtol < zero .OR. gtol < zero &
390  .OR. maxfev <= 0 .OR. factor <= zero) GO TO 300
391  IF (mode /= 2) GO TO 20
392  DO j = 1, n
393  IF (diag(j) <= zero) GO TO 300
394  END DO
395 
396  ! evaluate the function at the starting point and calculate its norm.
397 
398 20 iflag = 1
399  CALL fcn(m, n, x, fvec, iflag)
400  nfev = 1
401  IF (iflag < 0) GO TO 300
402  fnorm = enorm(m, fvec)
403 
404  ! initialize levenberg-marquardt parameter and iteration counter.
405 
406  par = zero
407  iter = 1
408 
409  ! beginning of the outer loop.
410 
411  ! calculate the jacobian matrix.
412 
413 30 iflag = 2
414  CALL fdjac2(fcn, m, n, x, fvec, fjac, iflag, epsfcn)
415  nfev = nfev + n
416  IF (iflag < 0) GO TO 300
417 
418  ! If requested, call fcn to enable printing of iterates.
419 
420  IF (nprint <= 0) GO TO 40
421  iflag = 0
422  IF (mod(iter-1,nprint) == 0) CALL fcn(m, n, x, fvec, iflag)
423  IF (iflag < 0) GO TO 300
424 
425  ! Compute the qr factorization of the jacobian.
426 
427 40 CALL qrfac(m, n, fjac, .true., ipvt, wa1, wa2)
428 
429  ! On the first iteration and if mode is 1, scale according
430  ! to the norms of the columns of the initial jacobian.
431 
432  IF (iter /= 1) GO TO 80
433  IF (mode == 2) GO TO 60
434  DO j = 1, n
435  diag(j) = wa2(j)
436  IF (wa2(j) == zero) diag(j) = one
437  END DO
438 
439  ! On the first iteration, calculate the norm of the scaled x
440  ! and initialize the step bound delta.
441 
442 60 wa3(1:n) = diag(1:n)*x(1:n)
443  xnorm = enorm(n, wa3)
444  delta = factor*xnorm
445  IF (delta == zero) delta = factor
446 
447  ! Form (q transpose)*fvec and store the first n components in qtf.
448 
449 80 wa4(1:m) = fvec(1:m)
450  DO j = 1, n
451  IF (fjac(j,j) == zero) GO TO 120
452  sum = dot_product( fjac(j:m,j), wa4(j:m) )
453  temp = -sum/fjac(j,j)
454  DO i = j, m
455  wa4(i) = wa4(i) + fjac(i,j)*temp
456  END DO
457 120 fjac(j,j) = wa1(j)
458  qtf(j) = wa4(j)
459  END DO
460 
461  ! compute the norm of the scaled gradient.
462 
463  gnorm = zero
464  IF (fnorm == zero) GO TO 170
465  DO j = 1, n
466  l = ipvt(j)
467  IF (wa2(l) == zero) cycle
468  sum = zero
469  DO i = 1, j
470  sum = sum + fjac(i,j)*(qtf(i)/fnorm)
471  END DO
472  gnorm = max(gnorm, abs(sum/wa2(l)))
473  END DO
474 
475  ! test for convergence of the gradient norm.
476 
477 170 IF (gnorm <= gtol) info = 4
478  IF (info /= 0) GO TO 300
479 
480  ! rescale if necessary.
481 
482  IF (mode == 2) GO TO 200
483  DO j = 1, n
484  diag(j) = max(diag(j), wa2(j))
485  END DO
486 
487  ! beginning of the inner loop.
488 
489  ! determine the Levenberg-Marquardt parameter.
490 
491 200 CALL lmpar(n, fjac, ipvt, diag, qtf, delta, par, wa1, wa2)
492 
493  ! store the direction p and x + p. calculate the norm of p.
494 
495  DO j = 1, n
496  wa1(j) = -wa1(j)
497  wa2(j) = x(j) + wa1(j)
498  wa3(j) = diag(j)*wa1(j)
499  END DO
500  pnorm = enorm(n, wa3)
501 
502  ! on the first iteration, adjust the initial step bound.
503 
504  IF (iter == 1) delta = min(delta, pnorm)
505 
506  ! evaluate the function at x + p and calculate its norm.
507 
508  iflag = 1
509  CALL fcn(m, n, wa2, wa4, iflag)
510  nfev = nfev + 1
511  IF (iflag < 0) GO TO 300
512  fnorm1 = enorm(m, wa4)
513 
514  ! compute the scaled actual reduction.
515 
516  actred = -one
517  IF (p1*fnorm1 < fnorm) actred = one - (fnorm1/fnorm)**2
518 
519  ! Compute the scaled predicted reduction and
520  ! the scaled directional derivative.
521 
522  DO j = 1, n
523  wa3(j) = zero
524  l = ipvt(j)
525  temp = wa1(l)
526  DO i = 1, j
527  wa3(i) = wa3(i) + fjac(i,j)*temp
528  END DO
529  END DO
530  temp1 = enorm(n,wa3)/fnorm
531  temp2 = (sqrt(par)*pnorm)/fnorm
532  prered = temp1**2 + temp2**2/p5
533  dirder = -(temp1**2 + temp2**2)
534 
535  ! compute the ratio of the actual to the predicted reduction.
536 
537  ratio = zero
538  IF (prered /= zero) ratio = actred/prered
539 
540  ! update the step bound.
541 
542  IF (ratio <= p25) THEN
543  IF (actred >= zero) temp = p5
544  IF (actred < zero) temp = p5*dirder/(dirder + p5*actred)
545  IF (p1*fnorm1 >= fnorm .OR. temp < p1) temp = p1
546  delta = temp*min(delta,pnorm/p1)
547  par = par/temp
548  ELSE
549  IF (par /= zero .AND. ratio < p75) GO TO 260
550  delta = pnorm/p5
551  par = p5*par
552  END IF
553 
554  ! test for successful iteration.
555 
556 260 IF (ratio < p0001) GO TO 290
557 
558  ! successful iteration. update x, fvec, and their norms.
559 
560  DO j = 1, n
561  x(j) = wa2(j)
562  wa2(j) = diag(j)*x(j)
563  END DO
564  fvec(1:m) = wa4(1:m)
565  xnorm = enorm(n, wa2)
566  fnorm = fnorm1
567  iter = iter + 1
568 
569  ! tests for convergence.
570 
571 290 IF (abs(actred) <= ftol .AND. prered <= ftol .AND. p5*ratio <= one) info = 1
572  IF (delta <= xtol*xnorm) info = 2
573  IF (abs(actred) <= ftol .AND. prered <= ftol &
574  .AND. p5*ratio <= one .AND. info == 2) info = 3
575  IF (info /= 0) GO TO 300
576 
577  ! tests for termination and stringent tolerances.
578 
579  IF (nfev >= maxfev) info = 5
580  IF (abs(actred) <= epsmch .AND. prered <= epsmch &
581  .AND. p5*ratio <= one) info = 6
582  IF (delta <= epsmch*xnorm) info = 7
583  IF (gnorm <= epsmch) info = 8
584  IF (info /= 0) GO TO 300
585 
586  ! end of the inner loop. repeat if iteration unsuccessful.
587 
588  IF (ratio < p0001) GO TO 200
589 
590  ! end of the outer loop.
591 
592  GO TO 30
593 
594  ! termination, either normal or user imposed.
595 
596 300 IF (iflag < 0) info = iflag
597  iflag = 0
598  IF (nprint > 0) CALL fcn(m, n, x, fvec, iflag)
599  RETURN
600 
601  ! last card of subroutine lmdif.
602 
603  END SUBROUTINE lmdif
604 
605 
606 
607  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
608  ! subroutine lmder1
609 
610  ! The purpose of lmder1 is to minimize the sum of the squares of
611  ! m nonlinear functions in n variables by a modification of the
612  ! levenberg-marquardt algorithm. This is done by using the more
613  ! general least-squares solver lmder. The user must provide a
614  ! subroutine which calculates the functions and the jacobian.
615 
616  ! the subroutine statement is
617 
618  ! subroutine lmder1(fcn, m, n, x, fvec, fjac, tol, info, ipvt)
619 
620  ! where
621 
622  ! fcn is the name of the user-supplied subroutine which
623  ! calculates the functions and the jacobian. fcn must
624  ! be declared in an interface statement in the user
625  ! calling program, and should be written as follows.
626 
627  ! subroutine fcn(m, n, x, fvec, fjac, iflag)
628  ! integer :: m, n, ldfjac, iflag
629  ! REAL :: x(:), fvec(:), fjac(:,:)
630  ! ----------
631  ! if iflag = 1 calculate the functions at x and
632  ! return this vector in fvec. do not alter fjac.
633  ! if iflag = 2 calculate the jacobian at x and
634  ! return this matrix in fjac. do not alter fvec.
635  ! ----------
636  ! return
637  ! end
638 
639  ! the value of iflag should not be changed by fcn unless
640  ! the user wants to terminate execution of lmder1.
641  ! in this case set iflag to a negative integer.
642 
643  ! m is a positive integer input variable set to the number of functions.
644 
645  ! n is a positive integer input variable set to the number
646  ! of variables. n must not exceed m.
647 
648  ! x is an array of length n. on input x must contain
649  ! an initial estimate of the solution vector. on output x
650  ! contains the final estimate of the solution vector.
651 
652  ! fvec is an output array of length m which contains
653  ! the functions evaluated at the output x.
654 
655  ! fjac is an output m by n array. the upper n by n submatrix
656  ! of fjac contains an upper triangular matrix r with
657  ! diagonal elements of nonincreasing magnitude such that
658 
659  ! t t t
660  ! p *(jac *jac)*p = r *r,
661 
662  ! where p is a permutation matrix and jac is the final calculated
663  ! Jacobian. Column j of p is column ipvt(j) (see below) of the
664  ! identity matrix. The lower trapezoidal part of fjac contains
665  ! information generated during the computation of r.
666 
667  ! ldfjac is a positive integer input variable not less than m
668  ! which specifies the leading dimension of the array fjac.
669 
670  ! tol is a nonnegative input variable. termination occurs
671  ! when the algorithm estimates either that the relative
672  ! error in the sum of squares is at most tol or that
673  ! the relative error between x and the solution is at most tol.
674 
675  ! info is an integer output variable. If the user has terminated
676  ! execution, info is set to the (negative) value of iflag.
677  ! See description of fcn. Otherwise, info is set as follows.
678 
679  ! info = 0 improper input parameters.
680 
681  ! info = 1 algorithm estimates that the relative error
682  ! in the sum of squares is at most tol.
683 
684  ! info = 2 algorithm estimates that the relative error
685  ! between x and the solution is at most tol.
686 
687  ! info = 3 conditions for info = 1 and info = 2 both hold.
688 
689  ! info = 4 fvec is orthogonal to the columns of the
690  ! jacobian to machine precision.
691 
692  ! info = 5 number of calls to fcn with iflag = 1 has reached 100*(n+1).
693 
694  ! info = 6 tol is too small. No further reduction in
695  ! the sum of squares is possible.
696 
697  ! info = 7 tol is too small. No further improvement in
698  ! the approximate solution x is possible.
699 
700  ! ipvt is an integer output array of length n. ipvt
701  ! defines a permutation matrix p such that jac*p = q*r,
702  ! where jac is the final calculated jacobian, q is
703  ! orthogonal (not stored), and r is upper triangular
704  ! with diagonal elements of nonincreasing magnitude.
705  ! column j of p is column ipvt(j) of the identity matrix.
706 
707  ! wa is a work array of length lwa.
708 
709  ! lwa is a positive integer input variable not less than 5*n+m.
710 
711  ! subprograms called
712 
713  ! user-supplied ...... fcn
714 
715  ! minpack-supplied ... lmder
716 
717  ! argonne national laboratory. minpack project. march 1980.
718  ! burton s. garbow, kenneth e. hillstrom, jorge j. more
719  !
720  ! N.B. Arguments LDFJAC, WA & LWA have been removed.
721  !
722  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
723 
724  SUBROUTINE lmder1(fcn, m, n, x, fvec, fjac, tol, info, ipvt)
725 
726  INTEGER, INTENT(IN) :: m
727  INTEGER, INTENT(IN) :: n
728  REAL, INTENT(IN OUT) :: x(:)
729  REAL, INTENT(OUT) :: fvec(:)
730  REAL, INTENT(IN OUT) :: fjac(:,:) ! fjac(ldfjac,n)
731  REAL, INTENT(IN) :: tol
732  INTEGER, INTENT(OUT) :: info
733  INTEGER, INTENT(IN OUT) :: ipvt(:)
734 
735 
736  ! EXTERNAL fcn
737 
738  INTERFACE
739  SUBROUTINE fcn(m, n, x, fvec, fjac, iflag)
740  IMPLICIT NONE
741  ! INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 307)
742  INTEGER, INTENT(IN) :: m, n
743  REAL, INTENT(IN) :: x(:)
744  REAL, INTENT(IN OUT) :: fvec(:)
745  REAL, INTENT(OUT) :: fjac(:,:)
746  INTEGER, INTENT(IN OUT) :: iflag
747  END SUBROUTINE fcn
748  END INTERFACE
749 
750  INTEGER :: maxfev, mode, nfev, njev, nprint
751  REAL :: ftol, gtol, xtol
752  REAL, PARAMETER :: factor = 100.0, zero = 0.0
753 
754  info = 0
755 
756  ! check the input parameters for errors.
757 
758  IF ( n <= 0 .OR. m < n .OR. tol < zero ) GO TO 10
759 
760  ! call lmder.
761 
762  maxfev = 100*(n + 1)
763  ftol = tol
764  xtol = tol
765  gtol = zero
766  mode = 1
767  nprint = 0
768  CALL lmder(fcn, m, n, x, fvec, fjac, ftol, xtol, gtol, maxfev, &
769  mode, factor, nprint, info, nfev, njev, ipvt)
770  IF (info == 8) info = 4
771 
772 10 RETURN
773 
774  ! last card of subroutine lmder1.
775 
776  END SUBROUTINE lmder1
777 
778 
779 
780  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
781  ! subroutine lmder
782 
783  ! the purpose of lmder is to minimize the sum of the squares of
784  ! m nonlinear functions in n variables by a modification of
785  ! the levenberg-marquardt algorithm. the user must provide a
786  ! subroutine which calculates the functions and the jacobian.
787 
788  ! the subroutine statement is
789 
790  ! subroutine lmder(fcn,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol,
791  ! maxfev,diag,mode,factor,nprint,info,nfev,
792  ! njev,ipvt,qtf,wa1,wa2,wa3,wa4)
793 
794  ! where
795 
796  ! fcn is the name of the user-supplied subroutine which
797  ! calculates the functions and the jacobian. fcn must
798  ! be declared in an external statement in the user
799  ! calling program, and should be written as follows.
800 
801  ! subroutine fcn(m,n,x,fvec,fjac,ldfjac,iflag)
802  ! integer m,n,ldfjac,iflag
803  ! REAL x(:),fvec(m),fjac(ldfjac,n)
804  ! ----------
805  ! if iflag = 1 calculate the functions at x and
806  ! return this vector in fvec. do not alter fjac.
807  ! if iflag = 2 calculate the jacobian at x and
808  ! return this matrix in fjac. Do not alter fvec.
809  ! ----------
810  ! return
811  ! end
812 
813  ! the value of iflag should not be changed by fcn unless
814  ! the user wants to terminate execution of lmder.
815  ! in this case set iflag to a negative integer.
816 
817  ! m is a positive integer input variable set to the number
818  ! of functions.
819 
820  ! n is a positive integer input variable set to the number
821  ! of variables. n must not exceed m.
822 
823  ! x is an array of length n. on input x must contain
824  ! an initial estimate of the solution vector. on output x
825  ! contains the final estimate of the solution vector.
826 
827  ! fvec is an output array of length m which contains
828  ! the functions evaluated at the output x.
829 
830  ! fjac is an output m by n array. the upper n by n submatrix
831  ! of fjac contains an upper triangular matrix r with
832  ! diagonal elements of nonincreasing magnitude such that
833 
834  ! t t t
835  ! p *(jac *jac)*p = r *r
836 
837  ! where p is a permutation matrix and jac is the final calculated
838  ! jacobian. Column j of p is column ipvt(j) (see below) of the
839  ! identity matrix. The lower trapezoidal part of fjac contains
840  ! information generated during the computation of r.
841 
842  ! ldfjac is a positive integer input variable not less than m
843  ! which specifies the leading dimension of the array fjac.
844 
845  ! ftol is a nonnegative input variable. Termination occurs when both
846  ! the actual and predicted relative reductions in the sum of squares
847  ! are at most ftol. Therefore, ftol measures the relative error
848  ! desired in the sum of squares.
849 
850  ! xtol is a nonnegative input variable. termination
851  ! occurs when the relative error between two consecutive
852  ! iterates is at most xtol. therefore, xtol measures the
853  ! relative error desired in the approximate solution.
854 
855  ! gtol is a nonnegative input variable. Termination occurs when the
856  ! cosine of the angle between fvec and any column of the jacobian is
857  ! at most gtol in absolute value. Therefore, gtol measures the
858  ! orthogonality desired between the function vector and the columns
859  ! of the jacobian.
860 
861  ! maxfev is a positive integer input variable. Termination occurs when
862  ! the number of calls to fcn with iflag = 1 has reached maxfev.
863 
864  ! diag is an array of length n. If mode = 1 (see below), diag is
865  ! internally set. If mode = 2, diag must contain positive entries
866  ! that serve as multiplicative scale factors for the variables.
867 
868  ! mode is an integer input variable. if mode = 1, the
869  ! variables will be scaled internally. if mode = 2,
870  ! the scaling is specified by the input diag. other
871  ! values of mode are equivalent to mode = 1.
872 
873  ! factor is a positive input variable used in determining the
874  ! initial step bound. this bound is set to the product of
875  ! factor and the euclidean norm of diag*x if nonzero, or else
876  ! to factor itself. in most cases factor should lie in the
877  ! interval (.1,100.).100. is a generally recommended value.
878 
879  ! nprint is an integer input variable that enables controlled printing
880  ! of iterates if it is positive. In this case, fcn is called with
881  ! iflag = 0 at the beginning of the first iteration and every nprint
882  ! iterations thereafter and immediately prior to return, with x, fvec,
883  ! and fjac available for printing. fvec and fjac should not be
884  ! altered. If nprint is not positive, no special calls of fcn with
885  ! iflag = 0 are made.
886 
887  ! info is an integer output variable. If the user has terminated
888  ! execution, info is set to the (negative) value of iflag.
889  ! See description of fcn. Otherwise, info is set as follows.
890 
891  ! info = 0 improper input parameters.
892 
893  ! info = 1 both actual and predicted relative reductions
894  ! in the sum of squares are at most ftol.
895 
896  ! info = 2 relative error between two consecutive iterates
897  ! is at most xtol.
898 
899  ! info = 3 conditions for info = 1 and info = 2 both hold.
900 
901  ! info = 4 the cosine of the angle between fvec and any column of
902  ! the jacobian is at most gtol in absolute value.
903 
904  ! info = 5 number of calls to fcn with iflag = 1 has reached maxfev.
905 
906  ! info = 6 ftol is too small. No further reduction in
907  ! the sum of squares is possible.
908 
909  ! info = 7 xtol is too small. No further improvement in
910  ! the approximate solution x is possible.
911 
912  ! info = 8 gtol is too small. fvec is orthogonal to the
913  ! columns of the jacobian to machine precision.
914 
915  ! nfev is an integer output variable set to the number of
916  ! calls to fcn with iflag = 1.
917 
918  ! njev is an integer output variable set to the number of
919  ! calls to fcn with iflag = 2.
920 
921  ! ipvt is an integer output array of length n. ipvt
922  ! defines a permutation matrix p such that jac*p = q*r,
923  ! where jac is the final calculated jacobian, q is
924  ! orthogonal (not stored), and r is upper triangular
925  ! with diagonal elements of nonincreasing magnitude.
926  ! column j of p is column ipvt(j) of the identity matrix.
927 
928  ! qtf is an output array of length n which contains
929  ! the first n elements of the vector (q transpose)*fvec.
930 
931  ! wa1, wa2, and wa3 are work arrays of length n.
932 
933  ! wa4 is a work array of length m.
934 
935  ! subprograms called
936 
937  ! user-supplied ...... fcn
938 
939  ! minpack-supplied ... dpmpar,enorm,lmpar,qrfac
940 
941  ! fortran-supplied ... ABS,MAX,MIN,SQRT,mod
942 
943  ! argonne national laboratory. minpack project. march 1980.
944  ! burton s. garbow, kenneth e. hillstrom, jorge j. more
945 
946  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
947 
948  SUBROUTINE lmder(fcn, m, n, x, fvec, fjac, ftol, xtol, gtol, maxfev, &
949  mode, factor, nprint, info, nfev, njev, ipvt)
950 
951  ! Code converted using TO_F90 by Alan Miller
952  ! Date: 1999-12-09 Time: 12:45:50
953 
954  ! N.B. Arguments LDFJAC, DIAG, QTF, WA1, WA2, WA3 & WA4 have been removed.
955 
956  INTEGER, INTENT(IN) :: m
957  INTEGER, INTENT(IN) :: n
958  REAL, INTENT(IN OUT) :: x(:)
959  REAL, INTENT(OUT) :: fvec(m)
960  REAL, INTENT(OUT) :: fjac(:,:) ! fjac(ldfjac,n)
961  REAL, INTENT(IN) :: ftol
962  REAL, INTENT(IN) :: xtol
963  REAL, INTENT(IN OUT) :: gtol
964  INTEGER, INTENT(IN OUT) :: maxfev
965  INTEGER, INTENT(IN) :: mode
966  REAL, INTENT(IN) :: factor
967  INTEGER, INTENT(IN) :: nprint
968  INTEGER, INTENT(OUT) :: info
969  INTEGER, INTENT(OUT) :: nfev
970  INTEGER, INTENT(OUT) :: njev
971  INTEGER, INTENT(OUT) :: ipvt(:)
972 
973  INTERFACE
974  SUBROUTINE fcn(m, n, x, fvec, fjac, iflag)
975  IMPLICIT NONE
976  ! INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 307)
977  INTEGER, INTENT(IN) :: m, n
978  REAL, INTENT(IN) :: x(:)
979  REAL, INTENT(IN OUT) :: fvec(:)
980  REAL, INTENT(OUT) :: fjac(:,:)
981  INTEGER, INTENT(IN OUT) :: iflag
982  END SUBROUTINE fcn
983  END INTERFACE
984 
985 
986  INTEGER :: i, iflag, iter, j, l
987  REAL :: actred, delta, dirder, epsmch, fnorm, fnorm1, gnorm, &
988  par, pnorm, prered, ratio, sum, temp, temp1, temp2, xnorm
989  REAL :: diag(n), qtf(n), wa1(n), wa2(n), wa3(n), wa4(m)
990  REAL, PARAMETER :: one = 1.0, p1 = 0.1, p5 = 0.5, &
991  p25 = 0.25, p75 = 0.75, p0001 = 0.0001, &
992  zero = 0.0
993 
994  ! epsmch is the machine precision.
995 
996  epsmch = epsilon(zero)
997 
998  info = 0
999  iflag = 0
1000  nfev = 0
1001  njev = 0
1002 
1003  ! check the input parameters for errors.
1004 
1005  IF (n <= 0 .OR. m < n .OR. ftol < zero .OR. xtol < zero .OR. gtol < zero &
1006  .OR. maxfev <= 0 .OR. factor <= zero) GO TO 300
1007  IF (mode /= 2) GO TO 20
1008  DO j = 1, n
1009  IF (diag(j) <= zero) GO TO 300
1010  END DO
1011 
1012  ! evaluate the function at the starting point and calculate its norm.
1013 
1014 20 iflag = 1
1015  CALL fcn(m, n, x, fvec, fjac, iflag)
1016  nfev = 1
1017  IF (iflag < 0) GO TO 300
1018  fnorm = enorm(m, fvec)
1019 
1020  ! initialize levenberg-marquardt parameter and iteration counter.
1021 
1022  par = zero
1023  iter = 1
1024 
1025  ! beginning of the outer loop.
1026 
1027  ! calculate the jacobian matrix.
1028 
1029 30 iflag = 2
1030  CALL fcn(m, n, x, fvec, fjac, iflag)
1031  njev = njev + 1
1032  IF (iflag < 0) GO TO 300
1033 
1034  ! if requested, call fcn to enable printing of iterates.
1035 
1036  IF (nprint <= 0) GO TO 40
1037  iflag = 0
1038  IF (mod(iter-1,nprint) == 0) CALL fcn(m, n, x, fvec, fjac, iflag)
1039  IF (iflag < 0) GO TO 300
1040 
1041  ! compute the qr factorization of the jacobian.
1042 
1043 40 CALL qrfac(m, n, fjac, .true., ipvt, wa1, wa2)
1044 
1045  ! on the first iteration and if mode is 1, scale according
1046  ! to the norms of the columns of the initial jacobian.
1047 
1048  IF (iter /= 1) GO TO 80
1049  IF (mode == 2) GO TO 60
1050  DO j = 1, n
1051  diag(j) = wa2(j)
1052  IF (wa2(j) == zero) diag(j) = one
1053  END DO
1054 
1055  ! on the first iteration, calculate the norm of the scaled x
1056  ! and initialize the step bound delta.
1057 
1058 60 wa3(1:n) = diag(1:n)*x(1:n)
1059  xnorm = enorm(n,wa3)
1060  delta = factor*xnorm
1061  IF (delta == zero) delta = factor
1062 
1063  ! form (q transpose)*fvec and store the first n components in qtf.
1064 
1065 80 wa4(1:m) = fvec(1:m)
1066  DO j = 1, n
1067  IF (fjac(j,j) == zero) GO TO 120
1068  sum = dot_product( fjac(j:m,j), wa4(j:m) )
1069  temp = -sum/fjac(j,j)
1070  DO i = j, m
1071  wa4(i) = wa4(i) + fjac(i,j)*temp
1072  END DO
1073 120 fjac(j,j) = wa1(j)
1074  qtf(j) = wa4(j)
1075  END DO
1076 
1077  ! compute the norm of the scaled gradient.
1078 
1079  gnorm = zero
1080  IF (fnorm == zero) GO TO 170
1081  DO j = 1, n
1082  l = ipvt(j)
1083  IF (wa2(l) == zero) cycle
1084  sum = zero
1085  DO i = 1, j
1086  sum = sum + fjac(i,j)*(qtf(i)/fnorm)
1087  END DO
1088  gnorm = max(gnorm,abs(sum/wa2(l)))
1089  END DO
1090 
1091  ! test for convergence of the gradient norm.
1092 
1093 170 IF (gnorm <= gtol) info = 4
1094  IF (info /= 0) GO TO 300
1095 
1096  ! rescale if necessary.
1097 
1098  IF (mode == 2) GO TO 200
1099  DO j = 1, n
1100  diag(j) = max(diag(j), wa2(j))
1101  END DO
1102 
1103  ! beginning of the inner loop.
1104 
1105  ! determine the levenberg-marquardt parameter.
1106 
1107 200 CALL lmpar(n, fjac, ipvt, diag, qtf, delta, par, wa1, wa2)
1108 
1109  ! store the direction p and x + p. calculate the norm of p.
1110 
1111  DO j = 1, n
1112  wa1(j) = -wa1(j)
1113  wa2(j) = x(j) + wa1(j)
1114  wa3(j) = diag(j)*wa1(j)
1115  END DO
1116  pnorm = enorm(n, wa3)
1117 
1118  ! on the first iteration, adjust the initial step bound.
1119 
1120  IF (iter == 1) delta = min(delta,pnorm)
1121 
1122  ! evaluate the function at x + p and calculate its norm.
1123 
1124  iflag = 1
1125  CALL fcn(m, n, wa2, wa4, fjac, iflag)
1126  nfev = nfev + 1
1127  IF (iflag < 0) GO TO 300
1128  fnorm1 = enorm(m, wa4)
1129 
1130  ! compute the scaled actual reduction.
1131 
1132  actred = -one
1133  IF (p1*fnorm1 < fnorm) actred = one - (fnorm1/fnorm)**2
1134 
1135  ! compute the scaled predicted reduction and
1136  ! the scaled directional derivative.
1137 
1138  DO j = 1, n
1139  wa3(j) = zero
1140  l = ipvt(j)
1141  temp = wa1(l)
1142  wa3(1:j) = wa3(1:j) + fjac(1:j,j)*temp
1143  END DO
1144  temp1 = enorm(n,wa3)/fnorm
1145  temp2 = (sqrt(par)*pnorm)/fnorm
1146  prered = temp1**2 + temp2**2/p5
1147  dirder = -(temp1**2 + temp2**2)
1148 
1149  ! compute the ratio of the actual to the predicted reduction.
1150 
1151  ratio = zero
1152  IF (prered /= zero) ratio = actred/prered
1153 
1154  ! update the step bound.
1155 
1156  IF (ratio <= p25) THEN
1157  IF (actred >= zero) temp = p5
1158  IF (actred < zero) temp = p5*dirder/(dirder + p5*actred)
1159  IF (p1*fnorm1 >= fnorm .OR. temp < p1) temp = p1
1160  delta = temp*min(delta, pnorm/p1)
1161  par = par/temp
1162  ELSE
1163  IF (par /= zero .AND. ratio < p75) GO TO 260
1164  delta = pnorm/p5
1165  par = p5*par
1166  END IF
1167 
1168  ! test for successful iteration.
1169 
1170 260 IF (ratio < p0001) GO TO 290
1171 
1172  ! successful iteration. update x, fvec, and their norms.
1173 
1174  DO j = 1, n
1175  x(j) = wa2(j)
1176  wa2(j) = diag(j)*x(j)
1177  END DO
1178  fvec(1:m) = wa4(1:m)
1179  xnorm = enorm(n,wa2)
1180  fnorm = fnorm1
1181  iter = iter + 1
1182 
1183  ! tests for convergence.
1184 
1185 290 IF (abs(actred) <= ftol .AND. prered <= ftol .AND. p5*ratio <= one) info = 1
1186  IF (delta <= xtol*xnorm) info = 2
1187  IF (abs(actred) <= ftol .AND. prered <= ftol &
1188  .AND. p5*ratio <= one .AND. info == 2) info = 3
1189  IF (info /= 0) GO TO 300
1190 
1191  ! tests for termination and stringent tolerances.
1192 
1193  IF (nfev >= maxfev) info = 5
1194  IF (abs(actred) <= epsmch .AND. prered <= epsmch &
1195  .AND. p5*ratio <= one) info = 6
1196  IF (delta <= epsmch*xnorm) info = 7
1197  IF (gnorm <= epsmch) info = 8
1198  IF (info /= 0) GO TO 300
1199 
1200  ! end of the inner loop. repeat if iteration unsuccessful.
1201 
1202  IF (ratio < p0001) GO TO 200
1203 
1204  ! end of the outer loop.
1205 
1206  GO TO 30
1207 
1208  ! termination, either normal or user imposed.
1209 
1210 300 IF (iflag < 0) info = iflag
1211  iflag = 0
1212  IF (nprint > 0) CALL fcn(m, n, x, fvec, fjac, iflag)
1213  RETURN
1214 
1215  ! last card of subroutine lmder.
1216 
1217  END SUBROUTINE lmder
1218 
1219 
1220 
1221  SUBROUTINE lmpar(n, r, ipvt, diag, qtb, delta, par, x, sdiag)
1222 
1223  ! Code converted using TO_F90 by Alan Miller
1224  ! Date: 1999-12-09 Time: 12:46:12
1225 
1226  ! N.B. Arguments LDR, WA1 & WA2 have been removed.
1227 
1228  INTEGER, INTENT(IN) :: n
1229  REAL, INTENT(IN OUT) :: r(:,:)
1230  INTEGER, INTENT(IN) :: ipvt(:)
1231  REAL, INTENT(IN) :: diag(:)
1232  REAL, INTENT(IN) :: qtb(:)
1233  REAL, INTENT(IN) :: delta
1234  REAL, INTENT(OUT) :: par
1235  REAL, INTENT(OUT) :: x(:)
1236  REAL, INTENT(OUT) :: sdiag(:)
1237 
1238  ! **********
1239 
1240  ! subroutine lmpar
1241 
1242  ! Given an m by n matrix a, an n by n nonsingular diagonal matrix d,
1243  ! an m-vector b, and a positive number delta, the problem is to determine a
1244  ! value for the parameter par such that if x solves the system
1245 
1246  ! a*x = b , sqrt(par)*d*x = 0 ,
1247 
1248  ! in the least squares sense, and dxnorm is the euclidean
1249  ! norm of d*x, then either par is zero and
1250 
1251  ! (dxnorm-delta) <= 0.1*delta ,
1252 
1253  ! or par is positive and
1254 
1255  ! abs(dxnorm-delta) <= 0.1*delta .
1256 
1257  ! This subroutine completes the solution of the problem if it is provided
1258  ! with the necessary information from the r factorization, with column
1259  ! qpivoting, of a. That is, if a*p = q*r, where p is a permutation matrix,
1260  ! q has orthogonal columns, and r is an upper triangular matrix with diagonal
1261  ! elements of nonincreasing magnitude, then lmpar expects the full upper
1262  ! triangle of r, the permutation matrix p, and the first n components of
1263  ! (q transpose)*b.
1264  ! On output lmpar also provides an upper triangular matrix s such that
1265 
1266  ! t t t
1267  ! p *(a *a + par*d*d)*p = s *s .
1268 
1269  ! s is employed within lmpar and may be of separate interest.
1270 
1271  ! Only a few iterations are generally needed for convergence of the algorithm.
1272  ! If, however, the limit of 10 iterations is reached, then the output par
1273  ! will contain the best value obtained so far.
1274 
1275  ! the subroutine statement is
1276 
1277  ! subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag, wa1,wa2)
1278 
1279  ! where
1280 
1281  ! n is a positive integer input variable set to the order of r.
1282 
1283  ! r is an n by n array. on input the full upper triangle
1284  ! must contain the full upper triangle of the matrix r.
1285  ! On output the full upper triangle is unaltered, and the
1286  ! strict lower triangle contains the strict upper triangle
1287  ! (transposed) of the upper triangular matrix s.
1288 
1289  ! ldr is a positive integer input variable not less than n
1290  ! which specifies the leading dimension of the array r.
1291 
1292  ! ipvt is an integer input array of length n which defines the
1293  ! permutation matrix p such that a*p = q*r. column j of p
1294  ! is column ipvt(j) of the identity matrix.
1295 
1296  ! diag is an input array of length n which must contain the
1297  ! diagonal elements of the matrix d.
1298 
1299  ! qtb is an input array of length n which must contain the first
1300  ! n elements of the vector (q transpose)*b.
1301 
1302  ! delta is a positive input variable which specifies an upper
1303  ! bound on the euclidean norm of d*x.
1304 
1305  ! par is a nonnegative variable. on input par contains an
1306  ! initial estimate of the levenberg-marquardt parameter.
1307  ! on output par contains the final estimate.
1308 
1309  ! x is an output array of length n which contains the least
1310  ! squares solution of the system a*x = b, sqrt(par)*d*x = 0,
1311  ! for the output par.
1312 
1313  ! sdiag is an output array of length n which contains the
1314  ! diagonal elements of the upper triangular matrix s.
1315 
1316  ! wa1 and wa2 are work arrays of length n.
1317 
1318  ! subprograms called
1319 
1320  ! minpack-supplied ... dpmpar,enorm,qrsolv
1321 
1322  ! fortran-supplied ... ABS,MAX,MIN,SQRT
1323 
1324  ! argonne national laboratory. minpack project. march 1980.
1325  ! burton s. garbow, kenneth e. hillstrom, jorge j. more
1326 
1327  ! **********
1328  INTEGER :: iter, j, jm1, jp1, k, l, nsing
1329  REAL :: dxnorm, dwarf, fp, gnorm, parc, parl, paru, sum, temp
1330  REAL :: wa1(n), wa2(n)
1331  REAL, PARAMETER :: p1 = 0.1, p001 = 0.001, zero = 0.0
1332 
1333  ! dwarf is the smallest positive magnitude.
1334 
1335  dwarf = tiny(zero)
1336 
1337  ! compute and store in x the gauss-newton direction. if the
1338  ! jacobian is rank-deficient, obtain a least squares solution.
1339 
1340  nsing = n
1341  DO j = 1, n
1342  wa1(j) = qtb(j)
1343  IF (r(j,j) == zero .AND. nsing == n) nsing = j - 1
1344  IF (nsing < n) wa1(j) = zero
1345  END DO
1346 
1347  DO k = 1, nsing
1348  j = nsing - k + 1
1349  wa1(j) = wa1(j)/r(j,j)
1350  temp = wa1(j)
1351  jm1 = j - 1
1352  wa1(1:jm1) = wa1(1:jm1) - r(1:jm1,j)*temp
1353  END DO
1354 
1355  DO j = 1, n
1356  l = ipvt(j)
1357  x(l) = wa1(j)
1358  END DO
1359 
1360  ! initialize the iteration counter.
1361  ! evaluate the function at the origin, and test
1362  ! for acceptance of the gauss-newton direction.
1363 
1364  iter = 0
1365  wa2(1:n) = diag(1:n)*x(1:n)
1366  dxnorm = enorm(n, wa2)
1367  fp = dxnorm - delta
1368  IF (fp <= p1*delta) GO TO 220
1369 
1370  ! if the jacobian is not rank deficient, the newton
1371  ! step provides a lower bound, parl, for the zero of
1372  ! the function. Otherwise set this bound to zero.
1373 
1374  parl = zero
1375  IF (nsing < n) GO TO 120
1376  DO j = 1, n
1377  l = ipvt(j)
1378  wa1(j) = diag(l)*(wa2(l)/dxnorm)
1379  END DO
1380  DO j = 1, n
1381  sum = dot_product( r(1:j-1,j), wa1(1:j-1) )
1382  wa1(j) = (wa1(j) - sum)/r(j,j)
1383  END DO
1384  temp = enorm(n,wa1)
1385  parl = ((fp/delta)/temp)/temp
1386 
1387  ! calculate an upper bound, paru, for the zero of the function.
1388 
1389 120 DO j = 1, n
1390  sum = dot_product( r(1:j,j), qtb(1:j) )
1391  l = ipvt(j)
1392  wa1(j) = sum/diag(l)
1393  END DO
1394  gnorm = enorm(n,wa1)
1395  paru = gnorm/delta
1396  IF (paru == zero) paru = dwarf/min(delta,p1)
1397 
1398  ! if the input par lies outside of the interval (parl,paru),
1399  ! set par to the closer endpoint.
1400 
1401  par = max(par,parl)
1402  par = min(par,paru)
1403  IF (par == zero) par = gnorm/dxnorm
1404 
1405  ! beginning of an iteration.
1406 
1407 150 iter = iter + 1
1408 
1409  ! evaluate the function at the current value of par.
1410 
1411  IF (par == zero) par = max(dwarf, p001*paru)
1412  temp = sqrt(par)
1413  wa1(1:n) = temp*diag(1:n)
1414  CALL qrsolv(n, r, ipvt, wa1, qtb, x, sdiag)
1415  wa2(1:n) = diag(1:n)*x(1:n)
1416  dxnorm = enorm(n, wa2)
1417  temp = fp
1418  fp = dxnorm - delta
1419 
1420  ! if the function is small enough, accept the current value
1421  ! of par. also test for the exceptional cases where parl
1422  ! is zero or the number of iterations has reached 10.
1423 
1424  IF (abs(fp) <= p1*delta .OR. parl == zero .AND. fp <= temp &
1425  .AND. temp < zero .OR. iter == 10) GO TO 220
1426 
1427  ! compute the newton correction.
1428 
1429  DO j = 1, n
1430  l = ipvt(j)
1431  wa1(j) = diag(l)*(wa2(l)/dxnorm)
1432  END DO
1433  DO j = 1, n
1434  wa1(j) = wa1(j)/sdiag(j)
1435  temp = wa1(j)
1436  jp1 = j + 1
1437  wa1(jp1:n) = wa1(jp1:n) - r(jp1:n,j)*temp
1438  END DO
1439  temp = enorm(n,wa1)
1440  parc = ((fp/delta)/temp)/temp
1441 
1442  ! depending on the sign of the function, update parl or paru.
1443 
1444  IF (fp > zero) parl = max(parl,par)
1445  IF (fp < zero) paru = min(paru,par)
1446 
1447  ! compute an improved estimate for par.
1448 
1449  par = max(parl, par+parc)
1450 
1451  ! end of an iteration.
1452 
1453  GO TO 150
1454 
1455  ! termination.
1456 
1457 220 IF (iter == 0) par = zero
1458  RETURN
1459 
1460  ! last card of subroutine lmpar.
1461 
1462  END SUBROUTINE lmpar
1463 
1464 
1465 
1466  SUBROUTINE qrfac(m, n, a, pivot, ipvt, rdiag, acnorm)
1467 
1468  ! Code converted using TO_F90 by Alan Miller
1469  ! Date: 1999-12-09 Time: 12:46:17
1470 
1471  ! N.B. Arguments LDA, LIPVT & WA have been removed.
1472 
1473  INTEGER, INTENT(IN) :: m
1474  INTEGER, INTENT(IN) :: n
1475  REAL, INTENT(IN OUT) :: a(:,:)
1476  LOGICAL, INTENT(IN) :: pivot
1477  INTEGER, INTENT(OUT) :: ipvt(:)
1478  REAL, INTENT(OUT) :: rdiag(:)
1479  REAL, INTENT(OUT) :: acnorm(:)
1480 
1481  ! **********
1482 
1483  ! subroutine qrfac
1484 
1485  ! This subroutine uses Householder transformations with column pivoting
1486  ! (optional) to compute a qr factorization of the m by n matrix a.
1487  ! That is, qrfac determines an orthogonal matrix q, a permutation matrix p,
1488  ! and an upper trapezoidal matrix r with diagonal elements of nonincreasing
1489  ! magnitude, such that a*p = q*r. The householder transformation for
1490  ! column k, k = 1,2,...,min(m,n), is of the form
1491 
1492  ! t
1493  ! i - (1/u(k))*u*u
1494 
1495  ! where u has zeros in the first k-1 positions. The form of this
1496  ! transformation and the method of pivoting first appeared in the
1497  ! corresponding linpack subroutine.
1498 
1499  ! the subroutine statement is
1500 
1501  ! subroutine qrfac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm, wa)
1502 
1503  ! N.B. 3 of these arguments have been omitted in this version.
1504 
1505  ! where
1506 
1507  ! m is a positive integer input variable set to the number of rows of a.
1508 
1509  ! n is a positive integer input variable set to the number of columns of a.
1510 
1511  ! a is an m by n array. On input a contains the matrix for
1512  ! which the qr factorization is to be computed. On output
1513  ! the strict upper trapezoidal part of a contains the strict
1514  ! upper trapezoidal part of r, and the lower trapezoidal
1515  ! part of a contains a factored form of q (the non-trivial
1516  ! elements of the u vectors described above).
1517 
1518  ! lda is a positive integer input variable not less than m
1519  ! which specifies the leading dimension of the array a.
1520 
1521  ! pivot is a logical input variable. If pivot is set true,
1522  ! then column pivoting is enforced. If pivot is set false,
1523  ! then no column pivoting is done.
1524 
1525  ! ipvt is an integer output array of length lipvt. ipvt
1526  ! defines the permutation matrix p such that a*p = q*r.
1527  ! Column j of p is column ipvt(j) of the identity matrix.
1528  ! If pivot is false, ipvt is not referenced.
1529 
1530  ! lipvt is a positive integer input variable. If pivot is false,
1531  ! then lipvt may be as small as 1. If pivot is true, then
1532  ! lipvt must be at least n.
1533 
1534  ! rdiag is an output array of length n which contains the
1535  ! diagonal elements of r.
1536 
1537  ! acnorm is an output array of length n which contains the norms of the
1538  ! corresponding columns of the input matrix a.
1539  ! If this information is not needed, then acnorm can coincide with rdiag.
1540 
1541  ! wa is a work array of length n. If pivot is false, then wa
1542  ! can coincide with rdiag.
1543 
1544  ! subprograms called
1545 
1546  ! minpack-supplied ... dpmpar,enorm
1547 
1548  ! fortran-supplied ... MAX,SQRT,MIN
1549 
1550  ! argonne national laboratory. minpack project. march 1980.
1551  ! burton s. garbow, kenneth e. hillstrom, jorge j. more
1552 
1553  ! **********
1554  INTEGER :: i, j, jp1, k, kmax, minmn
1555  REAL :: ajnorm, epsmch, sum, temp, wa(n)
1556  REAL, PARAMETER :: one = 1.0, p05 = 0.05, zero = 0.0
1557 
1558  ! epsmch is the machine precision.
1559 
1560  epsmch = epsilon(zero)
1561 
1562  ! compute the initial column norms and initialize several arrays.
1563 
1564  DO j = 1, n
1565  acnorm(j) = enorm(m,a(1:,j))
1566  rdiag(j) = acnorm(j)
1567  wa(j) = rdiag(j)
1568  IF (pivot) ipvt(j) = j
1569  END DO
1570 
1571  ! Reduce a to r with Householder transformations.
1572 
1573  minmn = min(m,n)
1574  DO j = 1, minmn
1575  IF (.NOT.pivot) GO TO 40
1576 
1577  ! Bring the column of largest norm into the pivot position.
1578 
1579  kmax = j
1580  DO k = j, n
1581  IF (rdiag(k) > rdiag(kmax)) kmax = k
1582  END DO
1583  IF (kmax == j) GO TO 40
1584  DO i = 1, m
1585  temp = a(i,j)
1586  a(i,j) = a(i,kmax)
1587  a(i,kmax) = temp
1588  END DO
1589  rdiag(kmax) = rdiag(j)
1590  wa(kmax) = wa(j)
1591  k = ipvt(j)
1592  ipvt(j) = ipvt(kmax)
1593  ipvt(kmax) = k
1594 
1595  ! Compute the Householder transformation to reduce the
1596  ! j-th column of a to a multiple of the j-th unit vector.
1597 
1598 40 ajnorm = enorm(m-j+1, a(j:,j))
1599  IF (ajnorm == zero) cycle
1600  IF (a(j,j) < zero) ajnorm = -ajnorm
1601  a(j:m,j) = a(j:m,j)/ajnorm
1602  a(j,j) = a(j,j) + one
1603 
1604  ! Apply the transformation to the remaining columns and update the norms.
1605 
1606  jp1 = j + 1
1607  DO k = jp1, n
1608  sum = dot_product( a(j:m,j), a(j:m,k) )
1609  temp = sum/a(j,j)
1610  a(j:m,k) = a(j:m,k) - temp*a(j:m,j)
1611  IF (.NOT.pivot .OR. rdiag(k) == zero) cycle
1612  temp = a(j,k)/rdiag(k)
1613  rdiag(k) = rdiag(k)*sqrt(max(zero, one-temp**2))
1614  IF (p05*(rdiag(k)/wa(k))**2 > epsmch) cycle
1615  rdiag(k) = enorm(m-j, a(jp1:,k))
1616  wa(k) = rdiag(k)
1617  END DO
1618  rdiag(j) = -ajnorm
1619  END DO
1620  RETURN
1621 
1622  ! last card of subroutine qrfac.
1623 
1624  END SUBROUTINE qrfac
1625 
1626 
1627 
1628  SUBROUTINE qrsolv(n, r, ipvt, diag, qtb, x, sdiag)
1629 
1630  ! N.B. Arguments LDR & WA have been removed.
1631 
1632  INTEGER, INTENT(IN) :: n
1633  REAL, INTENT(IN OUT) :: r(:,:)
1634  INTEGER, INTENT(IN) :: ipvt(:)
1635  REAL, INTENT(IN) :: diag(:)
1636  REAL, INTENT(IN) :: qtb(:)
1637  REAL, INTENT(OUT) :: x(:)
1638  REAL, INTENT(OUT) :: sdiag(:)
1639 
1640  ! **********
1641 
1642  ! subroutine qrsolv
1643 
1644  ! Given an m by n matrix a, an n by n diagonal matrix d, and an m-vector b,
1645  ! the problem is to determine an x which solves the system
1646 
1647  ! a*x = b , d*x = 0 ,
1648 
1649  ! in the least squares sense.
1650 
1651  ! This subroutine completes the solution of the problem if it is provided
1652  ! with the necessary information from the qr factorization, with column
1653  ! pivoting, of a. That is, if a*p = q*r, where p is a permutation matrix,
1654  ! q has orthogonal columns, and r is an upper triangular matrix with diagonal
1655  ! elements of nonincreasing magnitude, then qrsolv expects the full upper
1656  ! triangle of r, the permutation matrix p, and the first n components of
1657  ! (q transpose)*b. The system a*x = b, d*x = 0, is then equivalent to
1658 
1659  ! t t
1660  ! r*z = q *b , p *d*p*z = 0 ,
1661 
1662  ! where x = p*z. if this system does not have full rank,
1663  ! then a least squares solution is obtained. On output qrsolv
1664  ! also provides an upper triangular matrix s such that
1665 
1666  ! t t t
1667  ! p *(a *a + d*d)*p = s *s .
1668 
1669  ! s is computed within qrsolv and may be of separate interest.
1670 
1671  ! the subroutine statement is
1672 
1673  ! subroutine qrsolv(n, r, ldr, ipvt, diag, qtb, x, sdiag, wa)
1674 
1675  ! N.B. Arguments LDR and WA have been removed in this version.
1676 
1677  ! where
1678 
1679  ! n is a positive integer input variable set to the order of r.
1680 
1681  ! r is an n by n array. On input the full upper triangle must contain
1682  ! the full upper triangle of the matrix r.
1683  ! On output the full upper triangle is unaltered, and the strict lower
1684  ! triangle contains the strict upper triangle (transposed) of the
1685  ! upper triangular matrix s.
1686 
1687  ! ldr is a positive integer input variable not less than n
1688  ! which specifies the leading dimension of the array r.
1689 
1690  ! ipvt is an integer input array of length n which defines the
1691  ! permutation matrix p such that a*p = q*r. Column j of p
1692  ! is column ipvt(j) of the identity matrix.
1693 
1694  ! diag is an input array of length n which must contain the
1695  ! diagonal elements of the matrix d.
1696 
1697  ! qtb is an input array of length n which must contain the first
1698  ! n elements of the vector (q transpose)*b.
1699 
1700  ! x is an output array of length n which contains the least
1701  ! squares solution of the system a*x = b, d*x = 0.
1702 
1703  ! sdiag is an output array of length n which contains the
1704  ! diagonal elements of the upper triangular matrix s.
1705 
1706  ! wa is a work array of length n.
1707 
1708  ! subprograms called
1709 
1710  ! fortran-supplied ... ABS,SQRT
1711 
1712  ! argonne national laboratory. minpack project. march 1980.
1713  ! burton s. garbow, kenneth e. hillstrom, jorge j. more
1714 
1715  ! **********
1716  INTEGER :: i, j, k, kp1, l, nsing
1717  REAL :: cos, cotan, qtbpj, sin, sum, tan, temp, wa(n)
1718  REAL, PARAMETER :: p5 = 0.5, p25 = 0.25, zero = 0.0
1719 
1720  ! Copy r and (q transpose)*b to preserve input and initialize s.
1721  ! In particular, save the diagonal elements of r in x.
1722 
1723  DO j = 1, n
1724  r(j:n,j) = r(j,j:n)
1725  x(j) = r(j,j)
1726  wa(j) = qtb(j)
1727  END DO
1728 
1729  ! Eliminate the diagonal matrix d using a givens rotation.
1730 
1731  DO j = 1, n
1732 
1733  ! Prepare the row of d to be eliminated, locating the
1734  ! diagonal element using p from the qr factorization.
1735 
1736  l = ipvt(j)
1737  IF (diag(l) == zero) cycle
1738  sdiag(j:n) = zero
1739  sdiag(j) = diag(l)
1740 
1741  ! The transformations to eliminate the row of d modify only a single
1742  ! element of (q transpose)*b beyond the first n, which is initially zero.
1743 
1744  qtbpj = zero
1745  DO k = j, n
1746 
1747  ! Determine a givens rotation which eliminates the
1748  ! appropriate element in the current row of d.
1749 
1750  IF (sdiag(k) == zero) cycle
1751  IF (abs(r(k,k)) < abs(sdiag(k))) THEN
1752  cotan = r(k,k)/sdiag(k)
1753  sin = p5/sqrt(p25 + p25*cotan**2)
1754  cos = sin*cotan
1755  ELSE
1756  tan = sdiag(k)/r(k,k)
1757  cos = p5/sqrt(p25 + p25*tan**2)
1758  sin = cos*tan
1759  END IF
1760 
1761  ! Compute the modified diagonal element of r and
1762  ! the modified element of ((q transpose)*b,0).
1763 
1764  r(k,k) = cos*r(k,k) + sin*sdiag(k)
1765  temp = cos*wa(k) + sin*qtbpj
1766  qtbpj = -sin*wa(k) + cos*qtbpj
1767  wa(k) = temp
1768 
1769  ! Accumulate the tranformation in the row of s.
1770 
1771  kp1 = k + 1
1772  DO i = kp1, n
1773  temp = cos*r(i,k) + sin*sdiag(i)
1774  sdiag(i) = -sin*r(i,k) + cos*sdiag(i)
1775  r(i,k) = temp
1776  END DO
1777  END DO
1778 
1779  ! Store the diagonal element of s and restore
1780  ! the corresponding diagonal element of r.
1781 
1782  sdiag(j) = r(j,j)
1783  r(j,j) = x(j)
1784  END DO
1785 
1786  ! Solve the triangular system for z. If the system is singular,
1787  ! then obtain a least squares solution.
1788 
1789  nsing = n
1790  DO j = 1, n
1791  IF (sdiag(j) == zero .AND. nsing == n) nsing = j - 1
1792  IF (nsing < n) wa(j) = zero
1793  END DO
1794 
1795  DO k = 1, nsing
1796  j = nsing - k + 1
1797  sum = dot_product( r(j+1:nsing,j), wa(j+1:nsing) )
1798  wa(j) = (wa(j) - sum)/sdiag(j)
1799  END DO
1800 
1801  ! Permute the components of z back to components of x.
1802 
1803  DO j = 1, n
1804  l = ipvt(j)
1805  x(l) = wa(j)
1806  END DO
1807  RETURN
1808 
1809  ! last card of subroutine qrsolv.
1810 
1811  END SUBROUTINE qrsolv
1812 
1813 
1814 
1815  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1816  ! function enorm
1817 
1818  ! given an n-vector x, this function calculates the euclidean norm of x.
1819 
1820  ! the euclidean norm is computed by accumulating the sum of squares in
1821  ! three different sums. The sums of squares for the small and large
1822  ! components are scaled so that no overflows occur. Non-destructive
1823  ! underflows are permitted. Underflows and overflows do not occur in the
1824  ! computation of the unscaled sum of squares for the intermediate
1825  ! components. The definitions of small, intermediate and large components
1826  ! depend on two constants, rdwarf and rgiant. The main restrictions on
1827  ! these constants are that rdwarf**2 not underflow and rgiant**2 not
1828  ! overflow. The constants given here are suitable for every known computer.
1829 
1830  ! the function statement is
1831 
1832  ! REAL (dp) function enorm(n,x)
1833 
1834  ! where
1835 
1836  ! n is a positive integer input variable.
1837 
1838  ! x is an input array of length n.
1839 
1840  ! subprograms called
1841 
1842  ! fortran-supplied ... ABS,SQRT
1843 
1844  ! argonne national laboratory. minpack project. march 1980.
1845  ! burton s. garbow, kenneth e. hillstrom, jorge j. more
1846  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1847 
1848 
1849  FUNCTION enorm(n,x) RESULT(fn_val)
1850 
1851  INTEGER, INTENT(IN) :: n
1852  REAL, INTENT(IN) :: x(:)
1853  REAL :: fn_val
1854 
1855  INTEGER :: i
1856  REAL :: agiant, floatn, s1, s2, s3, xabs, x1max, x3max
1857  REAL, PARAMETER :: one = 1.0, zero = 0.0, rdwarf = 3.834e-20, &
1858  rgiant = 1.304e+19
1859 
1860  s1 = zero
1861  s2 = zero
1862  s3 = zero
1863  x1max = zero
1864  x3max = zero
1865  floatn = n
1866  agiant = rgiant/floatn
1867  DO i = 1, n
1868  xabs = abs(x(i))
1869  IF (xabs > rdwarf .AND. xabs < agiant) GO TO 70
1870  IF (xabs <= rdwarf) GO TO 30
1871 
1872  ! sum for large components.
1873 
1874  IF (xabs <= x1max) GO TO 10
1875  s1 = one + s1*(x1max/xabs)**2
1876  x1max = xabs
1877  GO TO 20
1878 
1879 10 s1 = s1 + (xabs/x1max)**2
1880 
1881 20 GO TO 60
1882 
1883  ! sum for small components.
1884 
1885 30 IF (xabs <= x3max) GO TO 40
1886  s3 = one + s3*(x3max/xabs)**2
1887  x3max = xabs
1888  GO TO 60
1889 
1890 40 IF (xabs /= zero) s3 = s3 + (xabs/x3max)**2
1891 
1892 60 cycle
1893 
1894  ! sum for intermediate components.
1895 
1896 70 s2 = s2 + xabs**2
1897  END DO
1898 
1899  ! calculation of norm.
1900 
1901  IF (s1 /= zero) THEN
1902  fn_val = x1max*sqrt(s1 + (s2/x1max)/x1max)
1903  ELSE
1904 
1905  IF (s2 == zero) THEN
1906  fn_val = x3max*sqrt(s3)
1907  ELSE
1908  IF (s2 >= x3max) THEN
1909  fn_val = sqrt(s2*(one + (x3max/s2)*(x3max*s3)))
1910  ELSE
1911  fn_val = sqrt(x3max*((s2/x3max) + (x3max*s3)))
1912  END IF
1913  END IF
1914  END IF
1915 
1916  END FUNCTION enorm
1917 
1918 
1919 
1920  SUBROUTINE fdjac2(fcn, m, n, x, fvec, fjac, iflag, epsfcn)
1921 
1922  ! Code converted using TO_F90 by Alan Miller
1923  ! Date: 1999-12-09 Time: 12:45:44
1924 
1925  ! N.B. Arguments LDFJAC & WA have been removed.
1926 
1927  INTEGER, INTENT(IN) :: m
1928  INTEGER, INTENT(IN) :: n
1929  REAL, INTENT(IN OUT) :: x(n)
1930  REAL, INTENT(IN) :: fvec(m)
1931  REAL, INTENT(OUT) :: fjac(:,:) ! fjac(ldfjac,n)
1932  INTEGER, INTENT(IN OUT) :: iflag
1933  REAL, INTENT(IN) :: epsfcn
1934 
1935  INTERFACE
1936  SUBROUTINE fcn(m, n, x, fvec, iflag)
1937  IMPLICIT NONE
1938  INTEGER, PARAMETER :: dp = selected_real_kind(15, 307)
1939  INTEGER, INTENT(IN) :: m, n
1940  REAL, INTENT(IN) :: x(:)
1941  REAL, INTENT(IN OUT) :: fvec(:)
1942  INTEGER, INTENT(IN OUT) :: iflag
1943  END SUBROUTINE fcn
1944  END INTERFACE
1945 
1946  ! **********
1947 
1948  ! subroutine fdjac2
1949 
1950  ! this subroutine computes a forward-difference approximation
1951  ! to the m by n jacobian matrix associated with a specified
1952  ! problem of m functions in n variables.
1953 
1954  ! the subroutine statement is
1955 
1956  ! subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
1957 
1958  ! where
1959 
1960  ! fcn is the name of the user-supplied subroutine which calculates the
1961  ! functions. fcn must be declared in an external statement in the user
1962  ! calling program, and should be written as follows.
1963 
1964  ! subroutine fcn(m,n,x,fvec,iflag)
1965  ! integer m,n,iflag
1966  ! REAL (dp) x(n),fvec(m)
1967  ! ----------
1968  ! calculate the functions at x and
1969  ! return this vector in fvec.
1970  ! ----------
1971  ! return
1972  ! end
1973 
1974  ! the value of iflag should not be changed by fcn unless
1975  ! the user wants to terminate execution of fdjac2.
1976  ! in this case set iflag to a negative integer.
1977 
1978  ! m is a positive integer input variable set to the number of functions.
1979 
1980  ! n is a positive integer input variable set to the number of variables.
1981  ! n must not exceed m.
1982 
1983  ! x is an input array of length n.
1984 
1985  ! fvec is an input array of length m which must contain the
1986  ! functions evaluated at x.
1987 
1988  ! fjac is an output m by n array which contains the
1989  ! approximation to the jacobian matrix evaluated at x.
1990 
1991  ! ldfjac is a positive integer input variable not less than m
1992  ! which specifies the leading dimension of the array fjac.
1993 
1994  ! iflag is an integer variable which can be used to terminate
1995  ! the execution of fdjac2. see description of fcn.
1996 
1997  ! epsfcn is an input variable used in determining a suitable step length
1998  ! for the forward-difference approximation. This approximation assumes
1999  ! that the relative errors in the functions are of the order of epsfcn.
2000  ! If epsfcn is less than the machine precision, it is assumed that the
2001  ! relative errors in the functions are of the order of the machine
2002  ! precision.
2003 
2004  ! wa is a work array of length m.
2005 
2006  ! subprograms called
2007 
2008  ! user-supplied ...... fcn
2009 
2010  ! minpack-supplied ... dpmpar
2011 
2012  ! fortran-supplied ... ABS,MAX,SQRT
2013 
2014  ! argonne national laboratory. minpack project. march 1980.
2015  ! burton s. garbow, kenneth e. hillstrom, jorge j. more
2016 
2017  ! **********
2018  INTEGER :: j
2019  REAL :: eps, epsmch, h, temp, wa(m)
2020  REAL, PARAMETER :: zero = 0.0
2021 
2022  ! epsmch is the machine precision.
2023 
2024  epsmch = epsilon(zero)
2025 
2026  eps = sqrt(max(epsfcn, epsmch))
2027 
2028  DO j = 1, n
2029  temp = x(j)
2030  h = eps*abs(temp)
2031  IF (h == zero) h = eps
2032  x(j) = temp + h
2033  CALL fcn(m, n, x, wa, iflag)
2034  IF (iflag < 0) EXIT
2035  x(j) = temp
2036  fjac(1:m,j) = (wa(1:m) - fvec(1:m))/h
2037  END DO
2038 
2039  RETURN
2040 
2041  ! last card of subroutine fdjac2.
2042 
2043  END SUBROUTINE fdjac2
2044 
2045 
2046 END MODULE levenberg_marquardt