11 MODULE levenberg_marquardt
14 INTEGER,
PARAMETER :: dp = selected_real_kind(15, 307)
17 PUBLIC :: dp, lmdif1, lmdif, lmder1, lmder, enorm
118 SUBROUTINE lmdif1(fcn, m, n, x, fvec, tol, epsfcn, factor, info, iwa)
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(:)
133 SUBROUTINE fcn(m, n, x, fvec, iflag)
136 INTEGER,
INTENT(IN) :: m, n
137 REAL,
INTENT(IN) :: x(:)
138 REAL,
INTENT(IN OUT) :: fvec(:)
139 INTEGER,
INTENT(IN OUT) :: iflag
143 INTEGER :: maxfev, mode, nfev, nprint
144 REAL :: ftol, gtol, xtol, fjac(m,n)
146 REAL,
PARAMETER :: zero = 0.0
153 IF (n <= 0 .OR. m < n .OR. tol < zero)
RETURN 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
168 END SUBROUTINE lmdif1
337 SUBROUTINE lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, &
338 mode, factor, nprint, info, nfev, fjac, ipvt)
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(:,:)
356 INTEGER,
INTENT(OUT) :: ipvt(:)
361 SUBROUTINE fcn(m, n, x, fvec, iflag)
364 INTEGER,
INTENT(IN) :: m, n
365 REAL,
INTENT(IN) :: x(:)
366 REAL,
INTENT(IN OUT) :: fvec(:)
367 INTEGER,
INTENT(IN OUT) :: iflag
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, &
381 epsmch = epsilon(zero)
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
393 IF (diag(j) <= zero)
GO TO 300
399 CALL fcn(m, n, x, fvec, iflag)
401 IF (iflag < 0)
GO TO 300
402 fnorm = enorm(m, fvec)
414 CALL fdjac2(fcn, m, n, x, fvec, fjac, iflag, epsfcn)
416 IF (iflag < 0)
GO TO 300
420 IF (nprint <= 0)
GO TO 40
422 IF (mod(iter-1,nprint) == 0)
CALL fcn(m, n, x, fvec, iflag)
423 IF (iflag < 0)
GO TO 300
427 40
CALL qrfac(m, n, fjac, .true., ipvt, wa1, wa2)
432 IF (iter /= 1)
GO TO 80
433 IF (mode == 2)
GO TO 60
436 IF (wa2(j) == zero) diag(j) = one
442 60 wa3(1:n) = diag(1:n)*x(1:n)
443 xnorm = enorm(n, wa3)
445 IF (delta == zero) delta = factor
449 80 wa4(1:m) = fvec(1:m)
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)
455 wa4(i) = wa4(i) + fjac(i,j)*temp
457 120 fjac(j,j) = wa1(j)
464 IF (fnorm == zero)
GO TO 170
467 IF (wa2(l) == zero) cycle
470 sum = sum + fjac(i,j)*(qtf(i)/fnorm)
472 gnorm = max(gnorm, abs(sum/wa2(l)))
477 170
IF (gnorm <= gtol) info = 4
478 IF (info /= 0)
GO TO 300
482 IF (mode == 2)
GO TO 200
484 diag(j) = max(diag(j), wa2(j))
491 200
CALL lmpar(n, fjac, ipvt, diag, qtf, delta, par, wa1, wa2)
497 wa2(j) = x(j) + wa1(j)
498 wa3(j) = diag(j)*wa1(j)
500 pnorm = enorm(n, wa3)
504 IF (iter == 1) delta = min(delta, pnorm)
509 CALL fcn(m, n, wa2, wa4, iflag)
511 IF (iflag < 0)
GO TO 300
512 fnorm1 = enorm(m, wa4)
517 IF (p1*fnorm1 < fnorm) actred = one - (fnorm1/fnorm)**2
527 wa3(i) = wa3(i) + fjac(i,j)*temp
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)
538 IF (prered /= zero) ratio = actred/prered
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)
549 IF (par /= zero .AND. ratio < p75)
GO TO 260
556 260
IF (ratio < p0001)
GO TO 290
562 wa2(j) = diag(j)*x(j)
565 xnorm = enorm(n, wa2)
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
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
588 IF (ratio < p0001)
GO TO 200
596 300
IF (iflag < 0) info = iflag
598 IF (nprint > 0)
CALL fcn(m, n, x, fvec, iflag)
724 SUBROUTINE lmder1(fcn, m, n, x, fvec, fjac, tol, info, ipvt)
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(:,:)
731 REAL,
INTENT(IN) :: tol
732 INTEGER,
INTENT(OUT) :: info
733 INTEGER,
INTENT(IN OUT) :: ipvt(:)
739 SUBROUTINE fcn(m, n, x, fvec, fjac, iflag)
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
750 INTEGER :: maxfev, mode, nfev, njev, nprint
751 REAL :: ftol, gtol, xtol
752 REAL,
PARAMETER :: factor = 100.0, zero = 0.0
758 IF ( n <= 0 .OR. m < n .OR. tol < zero )
GO TO 10
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
776 END SUBROUTINE lmder1
948 SUBROUTINE lmder(fcn, m, n, x, fvec, fjac, ftol, xtol, gtol, maxfev, &
949 mode, factor, nprint, info, nfev, njev, ipvt)
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(:,:)
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(:)
974 SUBROUTINE fcn(m, n, x, fvec, fjac, iflag)
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
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, &
996 epsmch = epsilon(zero)
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
1009 IF (diag(j) <= zero)
GO TO 300
1015 CALL fcn(m, n, x, fvec, fjac, iflag)
1017 IF (iflag < 0)
GO TO 300
1018 fnorm = enorm(m, fvec)
1030 CALL fcn(m, n, x, fvec, fjac, iflag)
1032 IF (iflag < 0)
GO TO 300
1036 IF (nprint <= 0)
GO TO 40
1038 IF (mod(iter-1,nprint) == 0)
CALL fcn(m, n, x, fvec, fjac, iflag)
1039 IF (iflag < 0)
GO TO 300
1043 40
CALL qrfac(m, n, fjac, .true., ipvt, wa1, wa2)
1048 IF (iter /= 1)
GO TO 80
1049 IF (mode == 2)
GO TO 60
1052 IF (wa2(j) == zero) diag(j) = one
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
1065 80 wa4(1:m) = fvec(1:m)
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)
1071 wa4(i) = wa4(i) + fjac(i,j)*temp
1073 120 fjac(j,j) = wa1(j)
1080 IF (fnorm == zero)
GO TO 170
1083 IF (wa2(l) == zero) cycle
1086 sum = sum + fjac(i,j)*(qtf(i)/fnorm)
1088 gnorm = max(gnorm,abs(sum/wa2(l)))
1093 170
IF (gnorm <= gtol) info = 4
1094 IF (info /= 0)
GO TO 300
1098 IF (mode == 2)
GO TO 200
1100 diag(j) = max(diag(j), wa2(j))
1107 200
CALL lmpar(n, fjac, ipvt, diag, qtf, delta, par, wa1, wa2)
1113 wa2(j) = x(j) + wa1(j)
1114 wa3(j) = diag(j)*wa1(j)
1116 pnorm = enorm(n, wa3)
1120 IF (iter == 1) delta = min(delta,pnorm)
1125 CALL fcn(m, n, wa2, wa4, fjac, iflag)
1127 IF (iflag < 0)
GO TO 300
1128 fnorm1 = enorm(m, wa4)
1133 IF (p1*fnorm1 < fnorm) actred = one - (fnorm1/fnorm)**2
1142 wa3(1:j) = wa3(1:j) + fjac(1:j,j)*temp
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)
1152 IF (prered /= zero) ratio = actred/prered
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)
1163 IF (par /= zero .AND. ratio < p75)
GO TO 260
1170 260
IF (ratio < p0001)
GO TO 290
1176 wa2(j) = diag(j)*x(j)
1178 fvec(1:m) = wa4(1:m)
1179 xnorm = enorm(n,wa2)
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
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
1202 IF (ratio < p0001)
GO TO 200
1210 300
IF (iflag < 0) info = iflag
1212 IF (nprint > 0)
CALL fcn(m, n, x, fvec, fjac, iflag)
1217 END SUBROUTINE lmder
1221 SUBROUTINE lmpar(n, r, ipvt, diag, qtb, delta, par, x, sdiag)
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(:)
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
1343 IF (r(j,j) == zero .AND. nsing == n) nsing = j - 1
1344 IF (nsing < n) wa1(j) = zero
1349 wa1(j) = wa1(j)/r(j,j)
1352 wa1(1:jm1) = wa1(1:jm1) - r(1:jm1,j)*temp
1365 wa2(1:n) = diag(1:n)*x(1:n)
1366 dxnorm = enorm(n, wa2)
1368 IF (fp <= p1*delta)
GO TO 220
1375 IF (nsing < n)
GO TO 120
1378 wa1(j) = diag(l)*(wa2(l)/dxnorm)
1381 sum = dot_product( r(1:j-1,j), wa1(1:j-1) )
1382 wa1(j) = (wa1(j) - sum)/r(j,j)
1385 parl = ((fp/delta)/temp)/temp
1390 sum = dot_product( r(1:j,j), qtb(1:j) )
1392 wa1(j) = sum/diag(l)
1394 gnorm = enorm(n,wa1)
1396 IF (paru == zero) paru = dwarf/min(delta,p1)
1403 IF (par == zero) par = gnorm/dxnorm
1411 IF (par == zero) par = max(dwarf, p001*paru)
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)
1424 IF (abs(fp) <= p1*delta .OR. parl == zero .AND. fp <= temp &
1425 .AND. temp < zero .OR. iter == 10)
GO TO 220
1431 wa1(j) = diag(l)*(wa2(l)/dxnorm)
1434 wa1(j) = wa1(j)/sdiag(j)
1437 wa1(jp1:n) = wa1(jp1:n) - r(jp1:n,j)*temp
1440 parc = ((fp/delta)/temp)/temp
1444 IF (fp > zero) parl = max(parl,par)
1445 IF (fp < zero) paru = min(paru,par)
1449 par = max(parl, par+parc)
1457 220
IF (iter == 0) par = zero
1462 END SUBROUTINE lmpar
1466 SUBROUTINE qrfac(m, n, a, pivot, ipvt, rdiag, acnorm)
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(:)
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
1560 epsmch = epsilon(zero)
1565 acnorm(j) = enorm(m,a(1:,j))
1566 rdiag(j) = acnorm(j)
1568 IF (pivot) ipvt(j) = j
1575 IF (.NOT.pivot)
GO TO 40
1581 IF (rdiag(k) > rdiag(kmax)) kmax = k
1583 IF (kmax == j)
GO TO 40
1589 rdiag(kmax) = rdiag(j)
1592 ipvt(j) = ipvt(kmax)
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
1608 sum = dot_product( a(j:m,j), a(j:m,k) )
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))
1624 END SUBROUTINE qrfac
1628 SUBROUTINE qrsolv(n, r, ipvt, diag, qtb, x, sdiag)
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(:)
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
1737 IF (diag(l) == zero) cycle
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)
1756 tan = sdiag(k)/r(k,k)
1757 cos = p5/sqrt(p25 + p25*tan**2)
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
1773 temp = cos*r(i,k) + sin*sdiag(i)
1774 sdiag(i) = -sin*r(i,k) + cos*sdiag(i)
1791 IF (sdiag(j) == zero .AND. nsing == n) nsing = j - 1
1792 IF (nsing < n) wa(j) = zero
1797 sum = dot_product( r(j+1:nsing,j), wa(j+1:nsing) )
1798 wa(j) = (wa(j) - sum)/sdiag(j)
1811 END SUBROUTINE qrsolv
1849 FUNCTION enorm(n,x)
RESULT(fn_val)
1851 INTEGER,
INTENT(IN) :: n
1852 REAL,
INTENT(IN) :: x(:)
1856 REAL :: agiant, floatn, s1, s2, s3, xabs, x1max, x3max
1857 REAL,
PARAMETER :: one = 1.0, zero = 0.0, rdwarf = 3.834e-20, &
1866 agiant = rgiant/floatn
1869 IF (xabs > rdwarf .AND. xabs < agiant)
GO TO 70
1870 IF (xabs <= rdwarf)
GO TO 30
1874 IF (xabs <= x1max)
GO TO 10
1875 s1 = one + s1*(x1max/xabs)**2
1879 10 s1 = s1 + (xabs/x1max)**2
1885 30
IF (xabs <= x3max)
GO TO 40
1886 s3 = one + s3*(x3max/xabs)**2
1890 40
IF (xabs /= zero) s3 = s3 + (xabs/x3max)**2
1896 70 s2 = s2 + xabs**2
1901 IF (s1 /= zero)
THEN 1902 fn_val = x1max*sqrt(s1 + (s2/x1max)/x1max)
1905 IF (s2 == zero)
THEN 1906 fn_val = x3max*sqrt(s3)
1908 IF (s2 >= x3max)
THEN 1909 fn_val = sqrt(s2*(one + (x3max/s2)*(x3max*s3)))
1911 fn_val = sqrt(x3max*((s2/x3max) + (x3max*s3)))
1920 SUBROUTINE fdjac2(fcn, m, n, x, fvec, fjac, iflag, epsfcn)
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(:,:)
1932 INTEGER,
INTENT(IN OUT) :: iflag
1933 REAL,
INTENT(IN) :: epsfcn
1936 SUBROUTINE fcn(m, n, x, fvec, iflag)
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
2019 REAL :: eps, epsmch, h, temp, wa(m)
2020 REAL,
PARAMETER :: zero = 0.0
2024 epsmch = epsilon(zero)
2026 eps = sqrt(max(epsfcn, epsmch))
2031 IF (h == zero) h = eps
2033 CALL fcn(m, n, x, wa, iflag)
2036 fjac(1:m,j) = (wa(1:m) - fvec(1:m))/h
2043 END SUBROUTINE fdjac2
2046 END MODULE levenberg_marquardt