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
100 INTEGER :: maxfev, ml, mode, mu, nfev, nprint
102 REAL (dp),
PARAMETER :: factor = 1.0_dp, zero = 0.0_dp
109 IF (n <= 0 .OR. epsfcn < zero .OR. tol < zero)
GO TO 20
119 CALL hybrd(fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, &
120 factor, nprint, info, nfev)
121 IF (info == 5) info = 4
130 SUBROUTINE hybrd(fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, &
131 factor, nprint, info, nfev)
133 INTEGER,
INTENT(IN) :: n
134 REAL (dp),
INTENT(IN OUT) :: x(n)
135 REAL (dp),
INTENT(IN OUT) :: fvec(n)
136 REAL (dp),
INTENT(IN) :: xtol
137 INTEGER,
INTENT(IN OUT) :: maxfev
138 INTEGER,
INTENT(IN OUT) :: ml
139 INTEGER,
INTENT(IN) :: mu
140 REAL (dp),
INTENT(IN) :: epsfcn
141 REAL (dp),
INTENT(OUT) :: diag(n)
142 INTEGER,
INTENT(IN) :: mode
143 REAL (dp),
INTENT(IN) :: factor
144 INTEGER,
INTENT(IN OUT) :: nprint
145 INTEGER,
INTENT(OUT) :: info
146 INTEGER,
INTENT(OUT) :: nfev
150 SUBROUTINE fcn(N, X, FVEC, IFLAG)
152 INTEGER,
PARAMETER :: dp = selected_real_kind(14, 60)
153 INTEGER,
INTENT(IN) :: n
154 REAL (dp),
INTENT(IN) :: x(n)
155 REAL (dp),
INTENT(OUT) :: fvec(n)
156 INTEGER,
INTENT(IN OUT) :: iflag
305 INTEGER :: i, iflag, iter, j, jm1, l, lr, msum, ncfail, ncsuc, nslow1, &
308 LOGICAL :: jeval, sing
309 REAL (dp) :: actred, delta, epsmch, fnorm, fnorm1, pnorm, prered, &
310 ratio, sum, temp, xnorm
311 REAL (dp),
PARAMETER :: one = 1.0_dp, p1 = 0.1_dp, p5 = 0.5_dp, &
312 p001 = 0.001_dp, p0001 = 0.0001_dp, zero = 0.0_dp
315 REAL (dp) :: fjac(n,n), r(n*(n+1)/2), qtf(n), wa1(n), wa2(n), &
320 epsmch = epsilon(1.0_dp)
329 IF (n > 0 .AND. xtol >= zero .AND. maxfev > 0 .AND. ml >= 0 .AND. mu >= &
330 0 .AND. factor > zero )
THEN 338 CALL fcn(n, x, fvec, iflag)
341 fnorm = enorm(n, fvec)
345 msum = min(ml+mu+1,n)
362 CALL fdjac1(fcn, n, x, fvec, fjac, n, iflag, ml, mu, epsfcn, wa1, wa2)
368 CALL qrfac(n, n, fjac, n, .false., iwa, 1, wa1, wa2, wa3)
377 IF (wa2(j) == zero) diag(j) = one
384 wa3(1:n) = diag(1:n) * x(1:n)
385 xnorm = enorm(n, wa3)
386 delta = factor * xnorm
387 IF (delta == zero) delta = factor
394 IF (fjac(j,j) /= zero)
THEN 397 sum = sum + fjac(i,j) * qtf(i)
399 temp = -sum / fjac(j,j)
401 qtf(i) = qtf(i) + fjac(i,j) * temp
419 IF (wa1(j) == zero) sing = .true.
424 CALL qform(n, n, fjac, n, wa1)
430 diag(j) = max(diag(j), wa2(j))
438 120
IF (nprint > 0)
THEN 440 IF (mod(iter-1, nprint) == 0)
CALL fcn(n, x, fvec, iflag)
441 IF (iflag < 0)
GO TO 190
446 CALL dogleg(n, r, lr, diag, qtf, delta, wa1, wa2, wa3)
452 wa2(j) = x(j) + wa1(j)
453 wa3(j) = diag(j) * wa1(j)
455 pnorm = enorm(n, wa3)
459 IF (iter == 1) delta = min(delta, pnorm)
464 CALL fcn(n, wa2, wa4, iflag)
467 fnorm1 = enorm(n, wa4)
472 IF (fnorm1 < fnorm) actred = one - (fnorm1/fnorm) ** 2
480 sum = sum + r(l) * wa1(j)
483 wa3(i) = qtf(i) + sum
487 IF (temp < fnorm) prered = one - (temp/fnorm) ** 2
492 IF (prered > zero) ratio = actred / prered
503 IF (ratio >= p5 .OR. ncsuc > 1) delta = max(delta,pnorm/p5)
504 IF (abs(ratio-one) <= p1) delta = pnorm / p5
509 IF (ratio >= p0001)
THEN 515 wa2(j) = diag(j) * x(j)
518 xnorm = enorm(n, wa2)
526 IF (actred >= p001) nslow1 = 0
527 IF (jeval) nslow2 = nslow2 + 1
528 IF (actred >= p1) nslow2 = 0
532 IF (delta <= xtol*xnorm .OR. fnorm == zero) info = 1
537 IF (nfev >= maxfev) info = 2
538 IF (p1*max(p1*delta, pnorm) <= epsmch*xnorm) info = 3
539 IF (nslow2 == 5) info = 4
540 IF (nslow1 == 10) info = 5
546 IF (ncfail /= 2)
THEN 554 sum = sum + fjac(i,j) * wa4(i)
556 wa2(j) = (sum-wa3(j)) / pnorm
557 wa1(j) = diag(j) * ((diag(j)*wa1(j))/pnorm)
558 IF (ratio >= p0001) qtf(j) = sum
563 CALL r1updt(n, n, r, lr, wa1, wa2, wa3, sing)
564 CALL r1mpyq(n, n, fjac, n, wa2, wa3)
565 CALL r1mpyq(1, n, qtf, 1, wa2, wa3)
585 190
IF (iflag < 0) info = iflag
587 IF (nprint > 0)
CALL fcn(n, x, fvec, iflag)
596 SUBROUTINE dogleg(n, r, lr, diag, qtb, delta, x, wa1, wa2)
598 INTEGER,
INTENT(IN) :: n
599 INTEGER,
INTENT(IN) :: lr
600 REAL (dp),
INTENT(IN) :: r(lr)
601 REAL (dp),
INTENT(IN) :: diag(n)
602 REAL (dp),
INTENT(IN) :: qtb(n)
603 REAL (dp),
INTENT(IN) :: delta
604 REAL (dp),
INTENT(IN OUT) :: x(n)
605 REAL (dp),
INTENT(OUT) :: wa1(n)
606 REAL (dp),
INTENT(OUT) :: wa2(n)
666 INTEGER :: i, j, jj, jp1, k, l
667 REAL (dp) :: alpha, bnorm, epsmch, gnorm, qnorm, sgnorm, sum, temp
671 epsmch = epsilon(1.0_dp)
675 jj = (n*(n+1)) / 2 + 1
684 sum = sum + r(l) * x(i)
689 IF (temp == 0.0)
THEN 692 temp = max(temp,abs(r(l)))
696 IF (temp == 0.0) temp = epsmch
698 x(j) = (qtb(j)-sum) / temp
705 wa2(j) = diag(j) * x(j)
707 qnorm = enorm(n, wa2)
708 IF (qnorm > delta)
THEN 717 wa1(i) = wa1(i) + r(l) * temp
720 wa1(j) = wa1(j) / diag(j)
726 gnorm = enorm(n, wa1)
728 alpha = delta / qnorm
729 IF (gnorm /= 0.0)
THEN 735 wa1(j) = (wa1(j)/gnorm) / diag(j)
741 sum = sum + r(l) * wa1(i)
747 sgnorm = (gnorm/temp) / temp
752 IF (sgnorm < delta)
THEN 758 bnorm = enorm(n, qtb)
759 temp = (bnorm/gnorm) * (bnorm/qnorm) * (sgnorm/delta)
760 temp = temp - (delta/qnorm) * (sgnorm/delta) ** 2 + sqrt(( &
761 temp-(delta/qnorm))**2+(1.0-(delta/qnorm)**2)*(1.0-( sgnorm/delta)**2))
762 alpha = ((delta/qnorm)*(1.0-(sgnorm/delta)**2)) / temp
769 temp = (1.0-alpha) * min(sgnorm,delta)
771 x(j) = temp * wa1(j) + alpha * x(j)
775 END SUBROUTINE dogleg
778 SUBROUTINE fdjac1(fcn, n, x, fvec, fjac, ldfjac, iflag, ml, mu, epsfcn, &
781 INTEGER,
INTENT(IN) :: n
782 REAL (dp),
INTENT(IN OUT) :: x(n)
783 REAL (dp),
INTENT(IN) :: fvec(n)
784 INTEGER,
INTENT(IN) :: ldfjac
785 REAL (dp),
INTENT(OUT) :: fjac(ldfjac,n)
786 INTEGER,
INTENT(IN OUT) :: iflag
787 INTEGER,
INTENT(IN) :: ml
788 INTEGER,
INTENT(IN) :: mu
789 REAL (dp),
INTENT(IN) :: epsfcn
790 REAL (dp),
INTENT(IN OUT) :: wa1(n)
791 REAL (dp),
INTENT(OUT) :: wa2(n)
795 SUBROUTINE fcn(N, X, FVEC, IFLAG)
797 INTEGER,
PARAMETER :: dp = selected_real_kind(14, 60)
798 INTEGER,
INTENT(IN) :: n
799 REAL (dp),
INTENT(IN) :: x(n)
800 REAL (dp),
INTENT(OUT) :: fvec(n)
801 INTEGER,
INTENT(IN OUT) :: iflag
887 INTEGER :: i, j, k, msum
888 REAL (dp) :: eps, epsmch, h, temp
889 REAL (dp),
PARAMETER :: zero = 0.0_dp
893 epsmch = epsilon(1.0_dp)
895 eps = sqrt(max(epsfcn, epsmch))
904 IF (h == zero) h = eps
906 CALL fcn(n, x, wa1, iflag)
910 fjac(i,j) = (wa1(i)-fvec(i)) / h
920 h = eps * abs(wa2(j))
921 IF (h == zero) h = eps
924 CALL fcn(n, x, wa1, iflag)
928 h = eps * abs(wa2(j))
929 IF (h == zero) h = eps
932 IF (i >= j-mu .AND. i <= j+ml) fjac(i,j) = (wa1(i)-fvec(i)) / h
941 END SUBROUTINE fdjac1
945 SUBROUTINE qform(m, n, q, ldq, wa)
947 INTEGER,
INTENT(IN) :: m
948 INTEGER,
INTENT(IN) :: n
949 INTEGER,
INTENT(IN) :: ldq
950 REAL (dp),
INTENT(OUT) :: q(ldq,m)
951 REAL (dp),
INTENT(OUT) :: wa(m)
989 INTEGER :: i, j, jm1, k, l, minmn, np1
990 REAL (dp) :: sum, temp
991 REAL (dp),
PARAMETER :: one = 1.0_dp, zero = 0.0_dp
1026 IF (wa(k) /= zero)
THEN 1030 sum = sum + q(i,j) * wa(i)
1034 q(i,j) = q(i,j) - temp * wa(i)
1043 END SUBROUTINE qform
1046 SUBROUTINE qrfac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm, wa)
1048 INTEGER,
INTENT(IN) :: m
1049 INTEGER,
INTENT(IN) :: n
1050 INTEGER,
INTENT(IN) :: lda
1051 REAL (dp),
INTENT(IN OUT) :: a(lda,n)
1052 LOGICAL,
INTENT(IN) :: pivot
1053 INTEGER,
INTENT(IN) :: lipvt
1054 INTEGER,
INTENT(OUT) :: ipvt(lipvt)
1055 REAL (dp),
INTENT(OUT) :: rdiag(n)
1056 REAL (dp),
INTENT(OUT) :: acnorm(n)
1057 REAL (dp),
INTENT(OUT) :: wa(n)
1131 INTEGER :: i, j, jp1, k, kmax, minmn
1132 REAL (dp) :: ajnorm, epsmch, sum, temp
1133 REAL (dp),
PARAMETER :: one = 1.0_dp, p05 = 0.05_dp, zero = 0.0_dp
1137 epsmch = epsilon(1.0_dp)
1142 acnorm(j) = enorm(m, a(1:,j))
1143 rdiag(j) = acnorm(j)
1145 IF (pivot) ipvt(j) = j
1158 IF (rdiag(k) > rdiag(kmax)) kmax = k
1166 rdiag(kmax) = rdiag(j)
1169 ipvt(j) = ipvt(kmax)
1177 ajnorm = enorm(m-j+1, a(j:,j))
1178 IF (ajnorm /= zero)
THEN 1179 IF (a(j,j) < zero) ajnorm = -ajnorm
1181 a(i,j) = a(i,j) / ajnorm
1183 a(j,j) = a(j,j) + one
1192 sum = sum + a(i,j) * a(i,k)
1196 a(i,k) = a(i,k) - temp * a(i,j)
1198 IF (.NOT.(.NOT.pivot.OR.rdiag(k) == zero))
THEN 1199 temp = a(j,k) / rdiag(k)
1200 rdiag(k) = rdiag(k) * sqrt(max(zero,one-temp**2))
1201 IF (p05*(rdiag(k)/wa(k))**2 <= epsmch)
THEN 1202 rdiag(k) = enorm(m-j, a(jp1:,k))
1215 END SUBROUTINE qrfac
1219 SUBROUTINE r1mpyq(m, n, a, lda, v, w)
1221 INTEGER,
INTENT(IN) :: m
1222 INTEGER,
INTENT(IN) :: n
1223 INTEGER,
INTENT(IN) :: lda
1224 REAL (dp),
INTENT(IN OUT) :: a(lda,n)
1225 REAL (dp),
INTENT(IN) :: v(n)
1226 REAL (dp),
INTENT(IN) :: w(n)
1274 INTEGER :: i, j, nmj, nm1
1275 REAL (dp) :: cos, sin, temp
1276 REAL (dp),
PARAMETER :: one = 1.0_dp
1284 IF (abs(v(j)) > one) cos = one / v(j)
1285 IF (abs(v(j)) > one) sin = sqrt(one-cos**2)
1286 IF (abs(v(j)) <= one) sin = v(j)
1287 IF (abs(v(j)) <= one) cos = sqrt(one-sin**2)
1289 temp = cos * a(i,j) - sin * a(i,n)
1290 a(i,n) = sin * a(i,j) + cos * a(i,n)
1298 IF (abs(w(j)) > one) cos = one / w(j)
1299 IF (abs(w(j)) > one) sin = sqrt(one-cos**2)
1300 IF (abs(w(j)) <= one) sin = w(j)
1301 IF (abs(w(j)) <= one) cos = sqrt(one-sin**2)
1303 temp = cos * a(i,j) + sin * a(i,n)
1304 a(i,n) = -sin * a(i,j) + cos * a(i,n)
1313 END SUBROUTINE r1mpyq
1317 SUBROUTINE r1updt(m, n, s, ls, u, v, w, sing)
1319 INTEGER,
INTENT(IN) :: m
1320 INTEGER,
INTENT(IN) :: n
1321 INTEGER,
INTENT(IN) :: ls
1322 REAL (dp),
INTENT(IN OUT) :: s(ls)
1323 REAL (dp),
INTENT(IN) :: u(m)
1324 REAL (dp),
INTENT(IN OUT) :: v(n)
1325 REAL (dp),
INTENT(OUT) :: w(m)
1326 LOGICAL,
INTENT(OUT) :: sing
1392 INTEGER :: i, j, jj, l, nmj, nm1
1393 REAL (dp) :: cos, cotan, giant, sin, tan, tau, temp
1394 REAL (dp),
PARAMETER :: one = 1.0_dp, p5 = 0.5_dp, p25 = 0.25_dp, zero = 0.0_dp
1398 giant = huge(1.0_dp)
1402 jj = (n*(2*m-n+1)) / 2 - (m-n)
1421 IF (v(j) /= zero)
THEN 1425 IF (abs(v(n)) < abs(v(j)))
THEN 1427 sin = p5 / sqrt(p25+p25*cotan**2)
1430 IF (abs(cos)*giant > one) tau = one / cos
1433 cos = p5 / sqrt(p25+p25*tan**2)
1441 v(n) = sin * v(j) + cos * v(n)
1448 temp = cos * s(l) - sin * w(i)
1449 w(i) = sin * s(l) + cos * w(i)
1460 w(i) = w(i) + v(n) * u(i)
1468 IF (w(j) /= zero)
THEN 1473 IF (abs(s(jj)) < abs(w(j)))
THEN 1474 cotan = s(jj) / w(j)
1475 sin = p5 / sqrt(p25 + p25*cotan**2)
1478 IF (abs(cos)*giant > one) tau = one / cos
1481 cos = p5 / sqrt(p25+p25*tan**2)
1490 temp = cos * s(l) + sin * w(i)
1491 w(i) = -sin * s(l) + cos * w(i)
1503 IF (s(jj) == zero) sing = .true.
1515 IF (s(jj) == zero) sing = .true.
1520 END SUBROUTINE r1updt
1523 FUNCTION enorm(n, x)
RESULT(fn_val)
1525 INTEGER,
INTENT(IN) :: n
1526 REAL (dp),
INTENT(IN) :: x(n)
1564 REAL (dp) :: agiant, floatn, s1, s2, s3, xabs, x1max, x3max
1565 REAL (dp),
PARAMETER :: rdwarf = 1.0d-100, rgiant = 1.0d+100
1573 agiant = rgiant / floatn
1576 IF (xabs <= rdwarf .OR. xabs >= agiant)
THEN 1577 IF (xabs > rdwarf)
THEN 1581 IF (xabs > x1max)
THEN 1582 s1 = 1.0_dp + s1 * (x1max/xabs) ** 2
1585 s1 = s1 + (xabs/x1max) ** 2
1591 IF (xabs > x3max)
THEN 1592 s3 = 1.0_dp + s3 * (x3max/xabs) ** 2
1595 IF (xabs /= 0.0_dp) s3 = s3 + (xabs/x3max) ** 2
1608 IF (s1 /= 0.0_dp)
THEN 1609 fn_val = x1max * sqrt(s1 + (s2/x1max)/x1max)
1611 IF (s2 /= 0.0_dp)
THEN 1612 IF (s2 >= x3max) fn_val = sqrt(s2*(1.0_dp + (x3max/s2)*(x3max*s3)))
1613 IF (s2 < x3max) fn_val = sqrt(x3max*((s2/x3max) + (x3max*s3)))
1615 fn_val = x3max * sqrt(s3)
1621 END MODULE solve_nonlin