18 SUBROUTINE uobyqa(n, x, rhobeg, rhoend, iprint, maxfun)
20 INTEGER,
INTENT(IN) :: n
21 REAL,
INTENT(IN OUT) :: x(:)
22 REAL,
INTENT(IN) :: rhobeg
23 REAL,
INTENT(IN) :: rhoend
24 INTEGER,
INTENT(IN) :: iprint
25 INTEGER,
INTENT(IN) :: maxfun
59 npt = (n*n + 3*n + 2) / 2
60 CALL uobyqb(n, x, rhobeg, rhoend, iprint, maxfun, npt)
66 SUBROUTINE uobyqb(n, x, rhobeg, rhoend, iprint, maxfun, npt)
68 INTEGER,
INTENT(IN) :: n
69 REAL,
INTENT(IN OUT) :: x(:)
70 REAL,
INTENT(IN) :: rhobeg
71 REAL,
INTENT(IN) :: rhoend
72 INTEGER,
INTENT(IN) :: iprint
73 INTEGER,
INTENT(IN) :: maxfun
74 INTEGER,
INTENT(IN) :: npt
77 SUBROUTINE calfun(n, x, f)
79 INTEGER,
PARAMETER :: dp = selected_real_kind(15, 307)
80 INTEGER,
INTENT(IN) :: n
81 REAL,
INTENT(IN) :: x(:)
82 REAL,
INTENT(OUT) :: f
88 REAL :: xbase(n), xopt(n), xnew(n), xpt(npt,n), pq(npt-1)
89 REAL :: pl(npt,npt-1), h(n,n), g(n), d(n), vlag(npt), w(npt)
111 REAL :: half = 0.5_dp, one = 1.0_dp, tol = 0.01_dp, two = 2.0_dp
112 REAL :: zero = 0.0_dp
113 REAL :: ddknew, delta, detrat, diff, distest, dnorm, errtol, estim
114 REAL :: evalue, f, fbase, fopt, fsave, ratio, rho, rhosq, sixthm
115 REAL :: sum, sumg, sumh, temp, tempa, tworsq, vmax, vquad, wmult
116 INTEGER :: i, ih, ip, iq, iw, j, jswitch, k, knew, kopt, ksave, ktemp
117 INTEGER :: nf, nftest, nnp, nptm
123 nftest = max(maxfun,1)
134 pl(1:npt,1:nptm) = zero
141 50 x(1:n) = xbase(1:n) + xpt(nf+1,1:n)
163 IF (jswitch > 0)
THEN 166 IF (w(j) < zero)
THEN 167 d(j) = (fsave+f-two*fbase) / rhosq
168 pq(j) = (fsave-f) / (two*rho)
169 pl(1,ih) = -two / rhosq
170 pl(nf-1,j) = half / rho
171 pl(nf-1,ih) = one / rhosq
173 pq(j) = (4.0d0*fsave-3.0d0*fbase-f) / (two*rho)
174 d(j) = (fbase+f-two*fsave) / rhosq
175 pl(1,j) = -1.5d0 / rho
176 pl(1,ih) = one / rhosq
177 pl(nf-1,j) = two / rho
178 pl(nf-1,ih) = -two / rhosq
181 pl(nf,j) = -half / rho
182 pl(nf,ih) = one / rhosq
196 xpt(nf+1,j) = two * rho
202 IF (nf < nnp)
GO TO 50
212 temp = one / (w(ip)*w(iq))
213 tempa = f - fbase - w(ip) * pq(ip) - w(iq) * pq(iq)
214 pq(ih) = (tempa - half*rhosq*(d(ip)+d(iq))) * temp
217 IF (w(ip) < zero) iw = iw + 1
220 IF (w(iq) < zero) iw = iw + 1
244 80 tworsq = (two*rho) ** 2
252 xopt(j) = xpt(kopt,j)
256 g(i) = g(i) + pq(ih) * xopt(j)
257 IF (i < j) g(j) = g(j) + pq(ih) * xopt(i)
266 CALL trstep(n, g, h, delta, tol, d, evalue)
269 temp = temp + d(i)**2
271 dnorm = min(delta,sqrt(temp))
273 IF (dnorm < half*rho)
THEN 275 errtol = half * evalue * rho * rho
276 IF (nf <= npt+9) errtol = zero
283 xnew(i) = xopt(i) + d(i)
284 x(i) = xbase(i) + xnew(i)
286 150
IF (nf >= nftest)
THEN 287 IF (iprint > 0)
WRITE(*, 5000)
292 IF (iprint == 3)
THEN 293 WRITE(*, 5100) nf, f, x(1:n)
295 IF (nf <= npt)
GO TO 70
296 IF (knew == -1)
GO TO 420
305 vquad = vquad + w(j) * pq(j)
308 w(ih) = d(i) * xnew(j) + d(j) * xopt(i)
309 IF (i == j) w(ih) = half * w(ih)
310 vquad = vquad + w(ih) * pq(ih)
316 temp = temp + w(j) * pl(k,j)
320 vlag(kopt) = vlag(kopt) + one
325 diff = f - fopt - vquad
330 temp = temp + (xpt(k,i)-xnew(i)) ** 2
333 sum = sum + abs(temp*temp*temp*vlag(k))
335 sixthm = max(sixthm, abs(diff)/sum)
343 xopt(1:n) = xnew(1:n)
350 IF (vquad >= zero)
THEN 351 IF (iprint > 0)
WRITE(*, 5200)
354 ratio = (f-fsave) / vquad
355 IF (ratio <= 0.1d0)
THEN 357 ELSE IF (ratio <= 0.7d0)
THEN 358 delta = max(half*delta,dnorm)
360 delta = max(delta, 1.25d0*dnorm, dnorm+rho)
362 IF (delta <= 1.5d0*rho) delta = rho
375 sum = sum + (xpt(k,i)-xopt(i)) ** 2
378 IF (sum > rhosq) temp = temp * (sum/rhosq) ** 1.5d0
379 IF (temp > detrat .AND. k /= ktemp)
THEN 385 IF (knew == 0)
GO TO 290
392 xpt(knew,i) = xnew(i)
394 temp = one / vlag(knew)
396 pl(knew,j) = temp * pl(knew,j)
397 pq(j) = pq(j) + diff * pl(knew,j)
403 pl(k,j) = pl(k,j) - temp * pl(knew,j)
416 IF (ksave > 0)
GO TO 90
417 IF (dnorm > two*rho)
GO TO 90
418 IF (ddknew > tworsq)
GO TO 90
426 w(k) = w(k) + (xpt(k,i)-xopt(i)) ** 2
432 IF (w(k) > distest)
THEN 450 g(j) = g(j) + temp * xopt(i)
452 g(i) = g(i) + temp * xopt(j)
453 sumh = sumh + temp * temp
457 sumh = sumh + half * temp * temp
464 IF (errtol > zero)
THEN 468 sumg = sumg + g(i) ** 2
470 estim = rho * (sqrt(sumg)+rho*sqrt(half*sumh))
471 wmult = sixthm * distest ** 1.5d0
472 IF (wmult*estim <= errtol)
GO TO 320
479 CALL lagmax(n, g, h, rho, d, xnew, vmax)
480 IF (errtol > zero)
THEN 481 IF (wmult*vmax <= errtol)
GO TO 320
485 IF (dnorm > rho)
GO TO 90
491 IF (rho > rhoend)
THEN 494 xbase(j) = xbase(j) + xopt(j)
496 xpt(k,j) = xpt(k,j) - xopt(j)
500 pq(i) = pq(i) + pq(ih) * xopt(j)
502 pq(j) = pq(j) + pq(ih) * xopt(i)
504 pl(k,j) = pl(k,j) + pl(k,ih) * xopt(i)
508 pl(k,i) = pl(k,i) + pl(k,ih) * xopt(j)
517 IF (ratio <= 16.0d0)
THEN 519 ELSE IF (ratio <= 250.0d0)
THEN 520 rho = sqrt(ratio) * rhoend
524 delta = max(delta,rho)
525 IF (iprint >= 2)
THEN 526 IF (iprint >= 3)
WRITE(*, 5300)
527 WRITE(*, 5400) rho, nf
528 WRITE(*, 5500) fopt, xbase(1:n)
536 IF (errtol >= zero)
GO TO 130
537 420
IF (fopt <= f)
THEN 539 x(i) = xbase(i) + xopt(i)
543 IF (iprint >= 1)
THEN 545 WRITE(*, 5500) f, x(1:n)
549 5000
FORMAT (/t5,
'Return from UOBYQA because CALFUN has been', &
550 ' called MAXFUN times')
551 5100
FORMAT (/t5,
'Function number',i6,
' F =', g18.10, &
552 ' The corresponding X is:'/ (t3, 5g15.6))
553 5200
FORMAT (/t5,
'Return from UOBYQA because a trust', &
554 ' region step has failed to reduce Q')
556 5400
FORMAT (/t5,
'New RHO =', g11.4,
' Number of function values =',i6)
557 5500
FORMAT (t5,
'Least value of F =', g23.15, &
558 ' The corresponding X is:'/ (t3, 5g15.6))
559 5600
FORMAT (/t5,
'At the return from UOBYQA', &
560 ' Number of function values =', i6)
561 END SUBROUTINE uobyqb
565 SUBROUTINE trstep(n, g, h, delta, tol, d, evalue)
567 INTEGER,
INTENT(IN) :: n
568 REAL,
INTENT(IN) :: g(:)
569 REAL,
INTENT(IN OUT) :: h(:,:)
570 REAL,
INTENT(IN) :: delta
571 REAL,
INTENT(IN) :: tol
572 REAL,
INTENT(OUT) :: d(:)
573 REAL,
INTENT(OUT) :: evalue
599 REAL :: gg(n), td(n), tn(n), w(n), piv(n), z(n)
601 REAL :: delsq, dhd, dnorm, dsq, dtg, dtz, gam, gnorm, gsq, hnorm
602 REAL :: par, parl, parlest, paru, paruest, phi, phil, phiu, pivksv
603 REAL :: pivot, posdef, scale, shfmax, shfmin, shift, slope, sum
604 REAL :: tdmin, temp, tempa, tempb, wsq, wwsq, wz, zsq
605 INTEGER :: i, iterc, j, jp, k, kp, kpp, ksav, ksave, nm
606 REAL :: one = 1.0_dp, two = 2.0_dp, zero = 0.0_dp
610 delsq = delta * delta
632 sum = sum + h(i,k) ** 2
635 IF (sum == zero)
THEN 640 tn(k) = sign(sqrt(sum+temp*temp),temp)
641 h(kp,k) = -sum / (temp+tn(k))
642 temp = sqrt(two/(sum+h(kp,k)**2))
652 z(i) = z(i) + h(i,j) * w(j)
653 z(j) = z(j) + h(i,j) * w(i)
655 wz = wz + w(j) * z(j)
657 wz = wz + w(n) * z(n)
659 td(j) = td(j) + w(j) * (wz*w(j)-two*z(j))
663 h(i,j) = h(i,j) - w(i) * z(j) - w(j) * (z(i)-wz*w(i))
675 gsq = gsq + g(i) ** 2
682 sum = sum + gg(i) * h(i,k)
685 gg(i) = gg(i) - sum * h(i,k)
692 hnorm = abs(td(1)) + abs(tn(1))
696 temp = abs(tn(i-1)) + abs(td(i)) + abs(tn(i))
697 hnorm = max(hnorm,temp)
698 tdmin = min(tdmin,td(i))
700 IF (hnorm == zero)
THEN 701 IF (gnorm == zero)
GO TO 420
702 scale = delta / gnorm
704 d(i) = -scale * gg(i)
711 parl = max(zero, -tdmin, gnorm/delta-hnorm)
721 160 iterc = iterc + 1
725 170
IF (piv(k) > zero)
THEN 726 piv(k+1) = td(k+1) + par - tn(k) ** 2 / piv(k)
728 IF (piv(k) < zero .OR. tn(k) /= zero)
GO TO 180
730 piv(k+1) = td(k+1) + par
734 IF (piv(k) >= zero)
THEN 735 IF (piv(k) == zero) ksav = k
740 IF (ksav == 0 .AND. gsq > zero)
GO TO 250
741 IF (gsq == zero)
THEN 742 IF (par == zero)
GO TO 380
745 IF (ksav == 0)
GO TO 210
754 IF (abs(tn(k)) <= abs(piv(k)))
THEN 759 IF (temp <= abs(piv(k)))
THEN 760 d(k+1) = sign(one,-tn(k))
761 dhd = piv(k) + temp - two * abs(tn(k))
763 d(k+1) = -tn(k) / temp
764 dhd = piv(k) + tn(k) * d(k+1)
766 dsq = one + d(k+1) ** 2
770 IF (tn(k) /= zero)
THEN 771 d(k) = -tn(k) * d(k+1) / piv(k)
772 dsq = dsq + d(k) ** 2
778 parlest = par - dhd / dsq
784 IF (gsq == zero) temp = temp * (one-tol)
785 IF (paruest > zero .AND. parlest >= temp)
THEN 786 dtg = dot_product( d(1:n), gg(1:n) )
787 scale = -sign(delta/sqrt(dsq),dtg)
788 d(1:n) = scale * d(1:n)
794 240
IF (paru == zero)
THEN 795 par = two * parlest + gnorm / delta
797 par = 0.5d0 * (parl+paru)
798 par = max(par,parlest)
800 IF (paruest > zero) par = min(par,paruest)
805 250 w(1) = -gg(1) / piv(1)
807 w(i) = (-gg(i)-tn(i-1)*w(i-1)) / piv(i)
811 d(i) = w(i) - tn(i) * d(i+1) / piv(i)
819 dsq = dsq + d(i) ** 2
820 wsq = wsq + piv(i) * w(i) ** 2
822 IF (par /= zero .OR. dsq > delsq)
THEN 827 phi = one / dnorm - one / delta
828 temp = tol * (one+par*dsq/wsq) - dsq * phi * phi
829 IF (temp >= zero)
THEN 830 scale = delta / dnorm
836 IF (iterc >= 2 .AND. par <= parl)
GO TO 380
837 IF (paru > zero .AND. par >= paru)
GO TO 380
843 IF (posdef == one)
THEN 844 IF (phi <= phil)
GO TO 380
845 slope = (phi-phil) / (par-parl)
846 parlest = par - phi / slope
849 IF (paru > zero) slope = (phiu-phi) / (paru-par)
850 temp = par - phi / slope
851 IF (paruest > zero) temp = min(temp,paruest)
861 IF (posdef == zero)
THEN 864 temp = -tn(i-1) * w(i-1)
865 w(i) = (sign(one,temp)+temp) / piv(i)
869 z(i) = w(i) - tn(i) * z(i+1) / piv(i)
875 wwsq = wwsq + piv(i) * w(i) ** 2
876 zsq = zsq + z(i) ** 2
877 dtz = dtz + d(i) * z(i)
882 tempa = abs(delsq-dsq)
883 tempb = sqrt(dtz*dtz+tempa*zsq)
884 gam = tempa / (sign(tempb,dtz)+dtz)
885 temp = tol * (wsq+par*delsq) - gam * gam * wwsq
886 IF (temp >= zero)
THEN 888 d(i) = d(i) + gam * z(i)
892 parlest = max(parlest,par-wwsq/zsq)
898 IF (paru > zero)
THEN 899 IF (phi >= phiu)
GO TO 380
900 slope = (phiu-phi) / (paru-par)
902 parlest = max(parlest,par-phi/slope)
904 IF (posdef == one)
THEN 905 slope = (phi-phil) / (par-parl)
906 paruest = par - phi / slope
920 pivot = td(k) - tn(k-1) ** 2 / pivot
921 shfmax = min(shfmax,pivot)
928 350 shift = 0.5d0 * (shfmin+shfmax)
932 360
IF (temp > zero)
THEN 935 temp = td(k+1) - shift - tn(k) ** 2 / temp
941 IF (k < ksave)
GO TO 370
943 IF (pivksv == zero)
GO TO 370
944 IF (piv(k)-temp < temp-pivksv)
THEN 949 shfmax = (shift*piv(k) - shfmin*temp) / (piv(k)-temp)
957 IF (shfmin <= 0.99d0*shfmax)
GO TO 350
967 sum = sum + d(i) * h(i,k)
970 d(i) = d(i) - sum * h(i,k)
977 END SUBROUTINE trstep
981 SUBROUTINE lagmax(n, g, h, rho, d, v, vmax)
983 INTEGER,
INTENT(IN) :: n
984 REAL,
INTENT(IN) :: g(:)
985 REAL,
INTENT(OUT) :: h(:,:)
986 REAL,
INTENT(IN) :: rho
987 REAL,
INTENT(OUT) :: d(:)
988 REAL,
INTENT(OUT) :: v(:)
989 REAL,
INTENT(OUT) :: vmax
1006 REAL :: half = 0.5_dp, one = 1.0_dp, zero = 0.0_dp
1007 REAL :: dd, dhd, dlin, dsq, gd, gg, ghg, gnorm, halfrt, hmax, ratio
1008 REAL :: scale, sum, sumv, temp, tempa, tempb, tempc, tempd, tempv
1009 REAL :: vhg, vhv, vhw, vlin, vmu, vnorm, vsq, vv, wcos, whw, wsin, wsq
1023 sum = sum + h(i,j) ** 2
1025 IF (sum > hmax)
THEN 1042 vsq = vsq + v(i) ** 2
1043 d(i) = dot_product( h(i,1:n), v(1:n) )
1044 vhv = vhv + v(i) * d(i)
1045 dsq = dsq + d(i) ** 2
1047 IF (vhv*vhv <= 0.9999d0*dsq*vsq)
THEN 1051 d(i) = d(i) - temp * v(i)
1052 wsq = wsq + d(i) ** 2
1055 ratio = sqrt(wsq/vsq)
1057 temp = dot_product( h(i,1:n), d(1:n) )
1058 whw = whw + temp * d(i)
1061 vhv = ratio * ratio * vhv
1063 temp = half * (whw-vhv)
1064 temp = temp + sign(sqrt(temp**2+vhw**2),whw+vhv)
1066 d(i) = vhw * v(i) + temp * d(i)
1079 gd = gd + g(i) * d(i)
1081 sum = dot_product( h(i,1:n), d(1:n) )
1082 dhd = dhd + sum * d(i)
1086 scale = sign(rho/sqrt(dd),gd*dhd)
1088 v(i) = d(i) - temp * g(i)
1093 IF (gnorm*dd <= 0.5d-2*rho*abs(dhd) .OR. vv/dd <= 1.0d-4)
THEN 1094 vmax = abs(scale*(gd + half*scale*dhd))
1106 sum = dot_product( h(i,1:n), g(1:n) )
1107 sumv = dot_product( h(i,1:n), v(1:n) )
1108 ghg = ghg + sum * g(i)
1109 vhg = vhg + sumv * g(i)
1110 vhv = vhv + sumv * v(i)
1114 vhg = vhg / (vnorm*gnorm)
1116 IF (abs(vhg) <= 0.01d0*max(abs(ghg),abs(vhv)))
THEN 1121 temp = half * (ghg-vhv)
1122 vmu = temp + sign(sqrt(temp**2+vhg**2),temp)
1123 temp = sqrt(vmu**2+vhg**2)
1127 tempa = wcos / gnorm
1128 tempb = wsin / vnorm
1129 tempc = wcos / vnorm
1130 tempd = wsin / gnorm
1132 d(i) = tempa * g(i) + tempb * v(i)
1133 v(i) = tempc * v(i) - tempd * g(i)
1139 dlin = wcos * gnorm / rho
1140 vlin = -wsin * gnorm / rho
1141 tempa = abs(dlin) + half * abs(vmu+vhv)
1142 tempb = abs(vlin) + half * abs(ghg-vmu)
1143 tempc = halfrt * (abs(dlin)+abs(vlin)) + 0.25d0 * abs(ghg+vhv)
1144 IF (tempa >= tempb .AND. tempa >= tempc)
THEN 1145 tempd = sign(rho,dlin*(vmu+vhv))
1147 ELSE IF (tempb >= tempc)
THEN 1149 tempv = sign(rho,vlin*(ghg-vmu))
1151 tempd = sign(halfrt*rho,dlin*(ghg+vhv))
1152 tempv = sign(halfrt*rho,vlin*(ghg+vhv))
1155 d(i) = tempd * d(i) + tempv * v(i)
1157 vmax = rho * rho * max(tempa,tempb,tempc)
1159 END SUBROUTINE lagmax
1161 END MODULE powell_optimize