7 INTEGER,
PARAMETER :: dp = selected_real_kind(14, 60)
15 SUBROUTINE hbrd(fcn, n, x, fvec, epsfcn, tol, info, diag)
20 INTEGER,
INTENT(IN) :: n
21 REAL (dp),
INTENT(IN OUT) :: x(n)
22 REAL (dp),
INTENT(IN OUT) :: fvec(n)
23 REAL (dp),
INTENT(IN) :: epsfcn
24 REAL (dp),
INTENT(IN) :: tol
25 INTEGER,
INTENT(OUT) :: info
26 REAL (dp),
INTENT(OUT) :: diag(n)
30 SUBROUTINE fcn(N, X, FVEC, IFLAG)
32 INTEGER,
PARAMETER :: dp = selected_real_kind(14, 60)
33 INTEGER,
INTENT(IN) :: n
34 REAL (dp),
INTENT(IN) :: x(n)
35 REAL (dp),
INTENT(OUT) :: fvec(n)
36 INTEGER,
INTENT(IN OUT) :: iflag
102 INTEGER :: maxfev, ml, mode, mu, nfev, nprint
104 REAL (dp),
PARAMETER :: factor = 1.0_dp, zero = 0.0_dp
111 IF (n <= 0 .OR. epsfcn < zero .OR. tol < zero)
GO TO 20
121 CALL hybrd(fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, &
122 factor, nprint, info, nfev)
123 IF (info == 5) info = 4
132 SUBROUTINE hybrd(fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, &
133 factor, nprint, info, nfev)
135 INTEGER,
INTENT(IN) :: n
136 REAL (dp),
INTENT(IN OUT) :: x(n)
137 REAL (dp),
INTENT(IN OUT) :: fvec(n)
138 REAL (dp),
INTENT(IN) :: xtol
139 INTEGER,
INTENT(IN OUT) :: maxfev
140 INTEGER,
INTENT(IN OUT) :: ml
141 INTEGER,
INTENT(IN) :: mu
142 REAL (dp),
INTENT(IN) :: epsfcn
143 REAL (dp),
INTENT(OUT) :: diag(n)
144 INTEGER,
INTENT(IN) :: mode
145 REAL (dp),
INTENT(IN) :: factor
146 INTEGER,
INTENT(IN OUT) :: nprint
147 INTEGER,
INTENT(OUT) :: info
148 INTEGER,
INTENT(OUT) :: nfev
152 SUBROUTINE fcn(N, X, FVEC, IFLAG)
154 INTEGER,
PARAMETER :: dp = selected_real_kind(14, 60)
155 INTEGER,
INTENT(IN) :: n
156 REAL (dp),
INTENT(IN) :: x(n)
157 REAL (dp),
INTENT(OUT) :: fvec(n)
158 INTEGER,
INTENT(IN OUT) :: iflag
307 INTEGER :: i, iflag, iter, j, jm1, l, lr, msum, ncfail, ncsuc, nslow1, &
310 LOGICAL :: jeval, sing
311 REAL (dp) :: actred, delta, epsmch, fnorm, fnorm1, pnorm, prered, &
312 ratio, sum, temp, xnorm
313 REAL (dp),
PARAMETER :: one = 1.0_dp, p1 = 0.1_dp, p5 = 0.5_dp, &
314 p001 = 0.001_dp, p0001 = 0.0001_dp, zero = 0.0_dp
317 REAL (dp) :: fjac(n,n), r(n*(n+1)/2), qtf(n), wa1(n), wa2(n), &
322 epsmch = epsilon(1.0_dp)
331 IF (n > 0 .AND. xtol >= zero .AND. maxfev > 0 .AND. ml >= 0 .AND. mu >= &
332 0 .AND. factor > zero )
THEN 340 CALL fcn(n, x, fvec, iflag)
343 fnorm = enorm(n, fvec)
347 msum = min(ml+mu+1,n)
364 CALL fdjac1(fcn, n, x, fvec, fjac, n, iflag, ml, mu, epsfcn, wa1, wa2)
370 CALL qrfac(n, n, fjac, n, .false., iwa, 1, wa1, wa2, wa3)
379 IF (wa2(j) == zero) diag(j) = one
386 wa3(1:n) = diag(1:n) * x(1:n)
387 xnorm = enorm(n, wa3)
388 delta = factor * xnorm
389 IF (delta == zero) delta = factor
396 IF (fjac(j,j) /= zero)
THEN 399 sum = sum + fjac(i,j) * qtf(i)
401 temp = -sum / fjac(j,j)
403 qtf(i) = qtf(i) + fjac(i,j) * temp
421 IF (wa1(j) == zero) sing = .true.
426 CALL qform(n, n, fjac, n, wa1)
432 diag(j) = max(diag(j), wa2(j))
440 120
IF (nprint > 0)
THEN 442 IF (mod(iter-1, nprint) == 0)
CALL fcn(n, x, fvec, iflag)
443 IF (iflag < 0)
GO TO 190
448 CALL dogleg(n, r, lr, diag, qtf, delta, wa1, wa2, wa3)
454 wa2(j) = x(j) + wa1(j)
455 wa3(j) = diag(j) * wa1(j)
457 pnorm = enorm(n, wa3)
461 IF (iter == 1) delta = min(delta, pnorm)
466 CALL fcn(n, wa2, wa4, iflag)
469 fnorm1 = enorm(n, wa4)
474 IF (fnorm1 < fnorm) actred = one - (fnorm1/fnorm) ** 2
482 sum = sum + r(l) * wa1(j)
485 wa3(i) = qtf(i) + sum
489 IF (temp < fnorm) prered = one - (temp/fnorm) ** 2
494 IF (prered > zero) ratio = actred / prered
505 IF (ratio >= p5 .OR. ncsuc > 1) delta = max(delta,pnorm/p5)
506 IF (abs(ratio-one) <= p1) delta = pnorm / p5
511 IF (ratio >= p0001)
THEN 517 wa2(j) = diag(j) * x(j)
520 xnorm = enorm(n, wa2)
528 IF (actred >= p001) nslow1 = 0
529 IF (jeval) nslow2 = nslow2 + 1
530 IF (actred >= p1) nslow2 = 0
534 IF (delta <= xtol*xnorm .OR. fnorm == zero) info = 1
539 IF (nfev >= maxfev) info = 2
540 IF (p1*max(p1*delta, pnorm) <= epsmch*xnorm) info = 3
541 IF (nslow2 == 5) info = 4
542 IF (nslow1 == 10) info = 5
548 IF (ncfail /= 2)
THEN 556 sum = sum + fjac(i,j) * wa4(i)
558 wa2(j) = (sum-wa3(j)) / pnorm
559 wa1(j) = diag(j) * ((diag(j)*wa1(j))/pnorm)
560 IF (ratio >= p0001) qtf(j) = sum
565 CALL r1updt(n, n, r, lr, wa1, wa2, wa3, sing)
566 CALL r1mpyq(n, n, fjac, n, wa2, wa3)
567 CALL r1mpyq(1, n, qtf, 1, wa2, wa3)
587 190
IF (iflag < 0) info = iflag
589 IF (nprint > 0)
CALL fcn(n, x, fvec, iflag)
598 SUBROUTINE dogleg(n, r, lr, diag, qtb, delta, x, wa1, wa2)
600 INTEGER,
INTENT(IN) :: n
601 INTEGER,
INTENT(IN) :: lr
602 REAL (dp),
INTENT(IN) :: r(lr)
603 REAL (dp),
INTENT(IN) :: diag(n)
604 REAL (dp),
INTENT(IN) :: qtb(n)
605 REAL (dp),
INTENT(IN) :: delta
606 REAL (dp),
INTENT(IN OUT) :: x(n)
607 REAL (dp),
INTENT(OUT) :: wa1(n)
608 REAL (dp),
INTENT(OUT) :: wa2(n)
668 INTEGER :: i, j, jj, jp1, k, l
669 REAL (dp) :: alpha, bnorm, epsmch, gnorm, qnorm, sgnorm, sum, temp
673 epsmch = epsilon(1.0_dp)
677 jj = (n*(n+1)) / 2 + 1
686 sum = sum + r(l) * x(i)
691 IF (temp == 0.0)
THEN 694 temp = max(temp,abs(r(l)))
698 IF (temp == 0.0) temp = epsmch
700 x(j) = (qtb(j)-sum) / temp
707 wa2(j) = diag(j) * x(j)
709 qnorm = enorm(n, wa2)
710 IF (qnorm > delta)
THEN 719 wa1(i) = wa1(i) + r(l) * temp
722 wa1(j) = wa1(j) / diag(j)
728 gnorm = enorm(n, wa1)
730 alpha = delta / qnorm
731 IF (gnorm /= 0.0)
THEN 737 wa1(j) = (wa1(j)/gnorm) / diag(j)
743 sum = sum + r(l) * wa1(i)
749 sgnorm = (gnorm/temp) / temp
754 IF (sgnorm < delta)
THEN 760 bnorm = enorm(n, qtb)
761 temp = (bnorm/gnorm) * (bnorm/qnorm) * (sgnorm/delta)
762 temp = temp - (delta/qnorm) * (sgnorm/delta) ** 2 + sqrt(( &
763 temp-(delta/qnorm))**2+(1.0-(delta/qnorm)**2)*(1.0-( sgnorm/delta)**2))
764 alpha = ((delta/qnorm)*(1.0-(sgnorm/delta)**2)) / temp
771 temp = (1.0-alpha) * min(sgnorm,delta)
773 x(j) = temp * wa1(j) + alpha * x(j)
777 END SUBROUTINE dogleg
780 SUBROUTINE fdjac1(fcn, n, x, fvec, fjac, ldfjac, iflag, ml, mu, epsfcn, &
783 INTEGER,
INTENT(IN) :: n
784 REAL (dp),
INTENT(IN OUT) :: x(n)
785 REAL (dp),
INTENT(IN) :: fvec(n)
786 INTEGER,
INTENT(IN) :: ldfjac
787 REAL (dp),
INTENT(OUT) :: fjac(ldfjac,n)
788 INTEGER,
INTENT(IN OUT) :: iflag
789 INTEGER,
INTENT(IN) :: ml
790 INTEGER,
INTENT(IN) :: mu
791 REAL (dp),
INTENT(IN) :: epsfcn
792 REAL (dp),
INTENT(IN OUT) :: wa1(n)
793 REAL (dp),
INTENT(OUT) :: wa2(n)
797 SUBROUTINE fcn(N, X, FVEC, IFLAG)
799 INTEGER,
PARAMETER :: dp = selected_real_kind(14, 60)
800 INTEGER,
INTENT(IN) :: n
801 REAL (dp),
INTENT(IN) :: x(n)
802 REAL (dp),
INTENT(OUT) :: fvec(n)
803 INTEGER,
INTENT(IN OUT) :: iflag
889 INTEGER :: i, j, k, msum
890 REAL (dp) :: eps, epsmch, h, temp
891 REAL (dp),
PARAMETER :: zero = 0.0_dp
895 epsmch = epsilon(1.0_dp)
897 eps = sqrt(max(epsfcn, epsmch))
906 IF (h == zero) h = eps
908 CALL fcn(n, x, wa1, iflag)
912 fjac(i,j) = (wa1(i)-fvec(i)) / h
922 h = eps * abs(wa2(j))
923 IF (h == zero) h = eps
926 CALL fcn(n, x, wa1, iflag)
930 h = eps * abs(wa2(j))
931 IF (h == zero) h = eps
934 IF (i >= j-mu .AND. i <= j+ml) fjac(i,j) = (wa1(i)-fvec(i)) / h
943 END SUBROUTINE fdjac1
947 SUBROUTINE qform(m, n, q, ldq, wa)
949 INTEGER,
INTENT(IN) :: m
950 INTEGER,
INTENT(IN) :: n
951 INTEGER,
INTENT(IN) :: ldq
952 REAL (dp),
INTENT(OUT) :: q(ldq,m)
953 REAL (dp),
INTENT(OUT) :: wa(m)
991 INTEGER :: i, j, jm1, k, l, minmn, np1
992 REAL (dp) :: sum, temp
993 REAL (dp),
PARAMETER :: one = 1.0_dp, zero = 0.0_dp
1028 IF (wa(k) /= zero)
THEN 1032 sum = sum + q(i,j) * wa(i)
1036 q(i,j) = q(i,j) - temp * wa(i)
1045 END SUBROUTINE qform
1048 SUBROUTINE qrfac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm, wa)
1050 INTEGER,
INTENT(IN) :: m
1051 INTEGER,
INTENT(IN) :: n
1052 INTEGER,
INTENT(IN) :: lda
1053 REAL (dp),
INTENT(IN OUT) :: a(lda,n)
1054 LOGICAL,
INTENT(IN) :: pivot
1055 INTEGER,
INTENT(IN) :: lipvt
1056 INTEGER,
INTENT(OUT) :: ipvt(lipvt)
1057 REAL (dp),
INTENT(OUT) :: rdiag(n)
1058 REAL (dp),
INTENT(OUT) :: acnorm(n)
1059 REAL (dp),
INTENT(OUT) :: wa(n)
1133 INTEGER :: i, j, jp1, k, kmax, minmn
1134 REAL (dp) :: ajnorm, epsmch, sum, temp
1135 REAL (dp),
PARAMETER :: one = 1.0_dp, p05 = 0.05_dp, zero = 0.0_dp
1139 epsmch = epsilon(1.0_dp)
1144 acnorm(j) = enorm(m, a(1:,j))
1145 rdiag(j) = acnorm(j)
1147 IF (pivot) ipvt(j) = j
1160 IF (rdiag(k) > rdiag(kmax)) kmax = k
1168 rdiag(kmax) = rdiag(j)
1171 ipvt(j) = ipvt(kmax)
1179 ajnorm = enorm(m-j+1, a(j:,j))
1180 IF (ajnorm /= zero)
THEN 1181 IF (a(j,j) < zero) ajnorm = -ajnorm
1183 a(i,j) = a(i,j) / ajnorm
1185 a(j,j) = a(j,j) + one
1194 sum = sum + a(i,j) * a(i,k)
1198 a(i,k) = a(i,k) - temp * a(i,j)
1200 IF (.NOT.(.NOT.pivot.OR.rdiag(k) == zero))
THEN 1201 temp = a(j,k) / rdiag(k)
1202 rdiag(k) = rdiag(k) * sqrt(max(zero,one-temp**2))
1203 IF (p05*(rdiag(k)/wa(k))**2 <= epsmch)
THEN 1204 rdiag(k) = enorm(m-j, a(jp1:,k))
1217 END SUBROUTINE qrfac
1221 SUBROUTINE r1mpyq(m, n, a, lda, v, w)
1223 INTEGER,
INTENT(IN) :: m
1224 INTEGER,
INTENT(IN) :: n
1225 INTEGER,
INTENT(IN) :: lda
1226 REAL (dp),
INTENT(IN OUT) :: a(lda,n)
1227 REAL (dp),
INTENT(IN) :: v(n)
1228 REAL (dp),
INTENT(IN) :: w(n)
1276 INTEGER :: i, j, nmj, nm1
1277 REAL (dp) :: cos, sin, temp
1278 REAL (dp),
PARAMETER :: one = 1.0_dp
1286 IF (abs(v(j)) > one) cos = one / v(j)
1287 IF (abs(v(j)) > one) sin = sqrt(one-cos**2)
1288 IF (abs(v(j)) <= one) sin = v(j)
1289 IF (abs(v(j)) <= one) cos = sqrt(one-sin**2)
1291 temp = cos * a(i,j) - sin * a(i,n)
1292 a(i,n) = sin * a(i,j) + cos * a(i,n)
1300 IF (abs(w(j)) > one) cos = one / w(j)
1301 IF (abs(w(j)) > one) sin = sqrt(one-cos**2)
1302 IF (abs(w(j)) <= one) sin = w(j)
1303 IF (abs(w(j)) <= one) cos = sqrt(one-sin**2)
1305 temp = cos * a(i,j) + sin * a(i,n)
1306 a(i,n) = -sin * a(i,j) + cos * a(i,n)
1315 END SUBROUTINE r1mpyq
1319 SUBROUTINE r1updt(m, n, s, ls, u, v, w, sing)
1321 INTEGER,
INTENT(IN) :: m
1322 INTEGER,
INTENT(IN) :: n
1323 INTEGER,
INTENT(IN) :: ls
1324 REAL (dp),
INTENT(IN OUT) :: s(ls)
1325 REAL (dp),
INTENT(IN) :: u(m)
1326 REAL (dp),
INTENT(IN OUT) :: v(n)
1327 REAL (dp),
INTENT(OUT) :: w(m)
1328 LOGICAL,
INTENT(OUT) :: sing
1394 INTEGER :: i, j, jj, l, nmj, nm1
1395 REAL (dp) :: cos, cotan, giant, sin, tan, tau, temp
1396 REAL (dp),
PARAMETER :: one = 1.0_dp, p5 = 0.5_dp, p25 = 0.25_dp, zero = 0.0_dp
1400 giant = huge(1.0_dp)
1404 jj = (n*(2*m-n+1)) / 2 - (m-n)
1423 IF (v(j) /= zero)
THEN 1427 IF (abs(v(n)) < abs(v(j)))
THEN 1429 sin = p5 / sqrt(p25+p25*cotan**2)
1432 IF (abs(cos)*giant > one) tau = one / cos
1435 cos = p5 / sqrt(p25+p25*tan**2)
1443 v(n) = sin * v(j) + cos * v(n)
1450 temp = cos * s(l) - sin * w(i)
1451 w(i) = sin * s(l) + cos * w(i)
1462 w(i) = w(i) + v(n) * u(i)
1470 IF (w(j) /= zero)
THEN 1475 IF (abs(s(jj)) < abs(w(j)))
THEN 1476 cotan = s(jj) / w(j)
1477 sin = p5 / sqrt(p25 + p25*cotan**2)
1480 IF (abs(cos)*giant > one) tau = one / cos
1483 cos = p5 / sqrt(p25+p25*tan**2)
1492 temp = cos * s(l) + sin * w(i)
1493 w(i) = -sin * s(l) + cos * w(i)
1505 IF (s(jj) == zero) sing = .true.
1517 IF (s(jj) == zero) sing = .true.
1522 END SUBROUTINE r1updt
1525 FUNCTION enorm(n, x)
RESULT(fn_val)
1527 INTEGER,
INTENT(IN) :: n
1528 REAL (dp),
INTENT(IN) :: x(n)
1566 REAL (dp) :: agiant, floatn, s1, s2, s3, xabs, x1max, x3max
1567 REAL (dp),
PARAMETER :: rdwarf = 1.0d-100, rgiant = 1.0d+100
1575 agiant = rgiant / floatn
1578 IF (xabs <= rdwarf .OR. xabs >= agiant)
THEN 1579 IF (xabs > rdwarf)
THEN 1583 IF (xabs > x1max)
THEN 1584 s1 = 1.0_dp + s1 * (x1max/xabs) ** 2
1587 s1 = s1 + (xabs/x1max) ** 2
1593 IF (xabs > x3max)
THEN 1594 s3 = 1.0_dp + s3 * (x3max/xabs) ** 2
1597 IF (xabs /= 0.0_dp) s3 = s3 + (xabs/x3max) ** 2
1610 IF (s1 /= 0.0_dp)
THEN 1611 fn_val = x1max * sqrt(s1 + (s2/x1max)/x1max)
1613 IF (s2 /= 0.0_dp)
THEN 1614 IF (s2 >= x3max) fn_val = sqrt(s2*(1.0_dp + (x3max/s2)*(x3max*s3)))
1615 IF (s2 < x3max) fn_val = sqrt(x3max*((s2/x3max) + (x3max*s3)))
1617 fn_val = x3max * sqrt(s3)
1623 END MODULE solve_nonlin