12 #include "../base/base_uses.f90"
16 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'powell'
21 INTEGER :: iprint = -1
23 INTEGER :: maxfun = -1
24 REAL(
dp) :: rhobeg = 0.0_dp, rhoend = 0.0_dp
25 REAL(
dp),
DIMENSION(:),
POINTER :: w => null()
26 REAL(
dp),
DIMENSION(:),
POINTER :: xopt => null()
28 INTEGER :: np = -1, nh = -1, nptm = -1, nftest = -1, idz = -1, itest = -1, nf = -1, nfm = -1, nfmm = -1, &
29 nfsav = -1, knew = -1, kopt = -1, ksave = -1, ktemp = -1
30 REAL(
dp) :: rhosq = 0.0_dp, recip = 0.0_dp, reciq = 0.0_dp, fbeg = 0.0_dp, &
31 fopt = 0.0_dp, diffa = 0.0_dp, xoptsq = 0.0_dp, &
32 rho = 0.0_dp, delta = 0.0_dp, dsq = 0.0_dp, dnorm = 0.0_dp, &
33 ratio = 0.0_dp, temp = 0.0_dp, tempq = 0.0_dp, beta = 0.0_dp, &
34 dx = 0.0_dp, vquad = 0.0_dp, diff = 0.0_dp, diffc = 0.0_dp, &
35 diffb = 0.0_dp, fsave = 0.0_dp, detrat = 0.0_dp, hdiag = 0.0_dp, &
36 distsq = 0.0_dp, gisq = 0.0_dp, gqsq = 0.0_dp, f = 0.0_dp, &
37 bstep = 0.0_dp, alpha = 0.0_dp, dstep = 0.0_dp
55 INTEGER,
INTENT(IN) :: n
56 REAL(
dp),
DIMENSION(*),
INTENT(INOUT) :: x
59 CHARACTER(len=*),
PARAMETER :: routinen =
'powell_optimize'
61 INTEGER :: handle, npt
63 CALL timeset(routinen, handle)
65 SELECT CASE (optstate%state)
68 ALLOCATE (optstate%w((npt + 13)*(npt + n) + 3*n*(n + 3)/2))
69 ALLOCATE (optstate%xopt(n))
73 CALL newuoa(n, x, optstate)
75 CALL newuoa(n, x, optstate)
77 IF (optstate%unit > 0)
THEN
78 WRITE (optstate%unit, *)
"POWELL| Exceeding maximum number of steps"
82 IF (optstate%unit > 0)
THEN
83 WRITE (optstate%unit, *)
"POWELL| Error in trust region"
87 IF (optstate%unit > 0)
THEN
88 WRITE (optstate%unit, *)
"POWELL| N out of range"
94 x(1:n) = optstate%xopt(1:n)
95 DEALLOCATE (optstate%w)
96 DEALLOCATE (optstate%xopt)
102 CALL timestop(handle)
112 SUBROUTINE newuoa(n, x, optstate)
114 INTEGER,
INTENT(IN) :: n
115 REAL(
dp),
DIMENSION(*),
INTENT(INOUT) :: x
118 INTEGER :: ibmat, id, ifv, igq, ihq, ipq, ivl, iw, &
119 ixb, ixn, ixo, ixp, izmat, maxfun, &
121 REAL(
dp) :: rhobeg, rhoend
123 maxfun = optstate%maxfun
124 rhobeg = optstate%rhobeg
125 rhoend = optstate%rhoend
164 IF (npt < n + 2 .OR. npt > ((n + 2)*np)/2)
THEN
178 izmat = ibmat + ndim*n
179 id = izmat + npt*nptm
187 CALL newuob(n, npt, x, rhobeg, rhoend, maxfun, optstate%w(ixb:), optstate%w(ixo:), &
188 optstate%w(ixn:), optstate%w(ixp:), optstate%w(ifv:), optstate%w(igq:), optstate%w(ihq:), &
189 optstate%w(ipq:), optstate%w(ibmat:), optstate%w(izmat:), ndim, optstate%w(id:), &
190 optstate%w(ivl:), optstate%w(iw:), optstate)
192 optstate%xopt(1:n) = optstate%w(ixb:ixb + n - 1) + optstate%w(ixo:ixo + n - 1)
194 END SUBROUTINE newuoa
221 SUBROUTINE newuob(n, npt, x, rhobeg, rhoend, maxfun, xbase, &
222 xopt, xnew, xpt, fval, gq, hq, pq, bmat, zmat, ndim, d, vlag, w, opt)
224 INTEGER,
INTENT(in) :: n, npt
225 REAL(
dp),
DIMENSION(1:n),
INTENT(inout) :: x
226 REAL(
dp),
INTENT(in) :: rhobeg, rhoend
227 INTEGER,
INTENT(in) :: maxfun
228 REAL(
dp),
DIMENSION(*),
INTENT(inout) :: xbase, xopt, xnew
229 REAL(
dp),
DIMENSION(npt, *), &
231 REAL(
dp),
DIMENSION(*),
INTENT(inout) :: fval, gq, hq, pq
232 INTEGER,
INTENT(in) :: ndim
233 REAL(
dp),
DIMENSION(npt, *), &
234 INTENT(inout) :: zmat
235 REAL(
dp),
DIMENSION(ndim, *), &
236 INTENT(inout) :: bmat
237 REAL(
dp),
DIMENSION(*),
INTENT(inout) :: d, vlag, w
240 INTEGER :: i, idz, ih, ip, ipt, itemp, &
241 itest, j, jp, jpt, k, knew, &
242 kopt, ksave, ktemp, nf, nfm, &
243 nfmm, nfsav, nftest, nh, np, &
245 REAL(
dp) :: alpha, beta, bstep, bsum, crvmin, delta, detrat, diff, diffa, &
246 diffb, diffc, distsq, dnorm, dsq, dstep, dx, f, fbeg, fopt, fsave, &
247 gisq, gqsq, half, hdiag, one, ratio, recip, reciq, rho, rhosq, sum, &
248 suma, sumb, sumz, temp, tempq, tenth, vquad, xipt, xjpt, xoptsq, zero
276 IF (opt%state == 1)
THEN
337 nftest = max(maxfun, 1)
339 IF (opt%state == 2)
GOTO 1000
366 rhosq = rhobeg*rhobeg
368 reciq = sqrt(half)/rhosq
374 IF (nfm >= 1 .AND. nfm <= n)
THEN
375 xpt(nf, nfm) = rhobeg
376 ELSE IF (nfm > n)
THEN
377 xpt(nf, nfmm) = -rhobeg
381 jpt = nfm - itemp*n - n
389 IF (fval(ipt + np) < fval(ipt + 1)) xipt = -xipt
391 IF (fval(jpt + np) < fval(jpt + 1)) xjpt = -xjpt
401 x(j) = xpt(nf, j) + xbase(j)
409 ELSE IF (f < fopt)
THEN
418 IF (nfm >= 1 .AND. nfm <= n)
THEN
419 gq(nfm) = (f - fbeg)/rhobeg
420 IF (npt < nf + n)
THEN
421 bmat(1, nfm) = -one/rhobeg
422 bmat(nf, nfm) = one/rhobeg
423 bmat(npt + nfm, nfm) = -half*rhosq
425 ELSE IF (nfm > n)
THEN
426 bmat(nf - n, nfmm) = half/rhobeg
427 bmat(nf, nfmm) = -half/rhobeg
428 zmat(1, nfmm) = -reciq - reciq
429 zmat(nf - n, nfmm) = reciq
430 zmat(nf, nfmm) = reciq
431 ih = (nfmm*(nfmm + 1))/2
432 temp = (fbeg - f)/rhobeg
433 hq(ih) = (gq(nfmm) - temp)/rhobeg
434 gq(nfmm) = half*(gq(nfmm) + temp)
441 ih = (ipt*(ipt - 1))/2 + jpt
442 IF (xipt < zero) ipt = ipt + n
443 IF (xjpt < zero) jpt = jpt + n
444 zmat(1, nfmm) = recip
445 zmat(nf, nfmm) = recip
446 zmat(ipt + 1, nfmm) = -recip
447 zmat(jpt + 1, nfmm) = -recip
448 hq(ih) = (fbeg - fval(ipt + 1) - fval(jpt + 1) + f)/(xipt*xjpt)
450 IF (nf < npt)
GOTO 50
462 xopt(i) = xpt(kopt, i)
463 xoptsq = xoptsq + xopt(i)**2
471 CALL trsapp(n, npt, xopt, xpt, gq, hq, pq, delta, d, w, w(np), w(np + n), w(np + 2*n), crvmin)
476 dnorm = min(delta, sqrt(dsq))
477 IF (dnorm < half*rho)
THEN
481 IF (delta <= 1.5_dp*rho) delta = rho
482 IF (nf <= nfsav + 2)
GOTO 460
483 temp = 0.125_dp*crvmin*rho*rho
484 IF (temp <= max(diffa, diffb, diffc))
GOTO 460
491 120
IF (dsq <= 1.0e-3_dp*xoptsq)
THEN
492 tempq = 0.25_dp*xoptsq
496 sum = sum + xpt(k, i)*xopt(i)
499 sum = sum - half*xoptsq
502 gq(i) = gq(i) + temp*xpt(k, i)
503 xpt(k, i) = xpt(k, i) - half*xopt(i)
505 w(i) = sum*xpt(k, i) + tempq*xopt(i)
508 bmat(ip, j) = bmat(ip, j) + vlag(i)*w(j) + w(i)*vlag(j)
518 sumz = sumz + zmat(i, k)
519 w(i) = w(npt + i)*zmat(i, k)
522 sum = tempq*sumz*xopt(j)
524 sum = sum + w(i)*xpt(i, j)
526 IF (k < idz) sum = -sum
529 bmat(i, j) = bmat(i, j) + sum*zmat(i, k)
535 IF (k < idz) temp = -temp
537 bmat(ip, j) = bmat(ip, j) + temp*vlag(j)
549 w(j) = w(j) + pq(k)*xpt(k, j)
550 xpt(k, j) = xpt(k, j) - half*xopt(j)
554 IF (i < j) gq(j) = gq(j) + hq(ih)*xopt(i)
555 gq(i) = gq(i) + hq(ih)*xopt(j)
556 hq(ih) = hq(ih) + w(i)*xopt(j) + xopt(i)*w(j)
557 bmat(npt + i, j) = bmat(npt + j, i)
561 xbase(j) = xbase(j) + xopt(j)
572 CALL biglag(n, npt, xopt, xpt, bmat, zmat, idz, ndim, knew, dstep, &
573 d, alpha, vlag, vlag(npt + 1), w, w(np), w(np + n))
584 suma = suma + xpt(k, j)*d(j)
585 sumb = sumb + xpt(k, j)*xopt(j)
586 sum = sum + bmat(k, j)*d(j)
588 w(k) = suma*(half*suma + sumb)
595 sum = sum + zmat(i, k)*w(i)
598 beta = beta + sum*sum
601 beta = beta - sum*sum
604 vlag(i) = vlag(i) + sum*zmat(i, k)
612 sum = sum + w(i)*bmat(i, j)
614 bsum = bsum + sum*d(j)
617 sum = sum + bmat(jp, k)*d(k)
620 bsum = bsum + sum*d(j)
621 dx = dx + d(j)*xopt(j)
623 beta = dx*dx + dsq*(xoptsq + dx + dx + half*dsq) + beta - bsum
624 vlag(kopt) = vlag(kopt) + one
631 temp = one + alpha*beta/vlag(knew)**2
632 IF (abs(temp) <= 0.8_dp)
THEN
633 CALL bigden(n, npt, xopt, xpt, bmat, zmat, idz, ndim, kopt, &
634 knew, d, w, vlag, beta, xnew, w(ndim + 1), w(6*ndim + 1))
641 xnew(i) = xopt(i) + d(i)
642 x(i) = xbase(i) + xnew(i)
645 310
IF (nf > nftest)
THEN
663 IF (nf <= npt)
GOTO 70
676 vquad = vquad + d(j)*gq(j)
679 temp = d(i)*xnew(j) + d(j)*xopt(i)
680 IF (i == j) temp = half*temp
681 vquad = vquad + temp*hq(ih)
685 vquad = vquad + pq(k)*w(k)
687 diff = f - fopt - vquad
691 IF (dnorm > rho) nfsav = nf
703 xoptsq = xoptsq + xopt(i)**2
707 IF (knew > 0)
GOTO 410
711 IF (vquad >= zero)
THEN
717 ratio = (f - fsave)/vquad
718 IF (ratio <= tenth)
THEN
720 ELSE IF (ratio <= 0.7_dp)
THEN
721 delta = max(half*delta, dnorm)
723 delta = max(half*delta, dnorm + dnorm)
725 IF (delta <= 1.5_dp*rho) delta = rho
729 rhosq = max(tenth*delta, rho)**2
740 IF (j < idz) temp = -one
741 hdiag = hdiag + temp*zmat(k, j)**2
743 temp = abs(beta*hdiag + vlag(k)**2)
746 distsq = distsq + (xpt(k, j) - xopt(j))**2
748 IF (distsq > rhosq) temp = temp*(distsq/rhosq)**3
749 IF (temp > detrat .AND. k /= ktemp)
THEN
754 IF (knew == 0)
GOTO 460
760 410
CALL update(n, npt, bmat, zmat, idz, ndim, vlag, beta, knew, w)
764 temp = pq(knew)*xpt(knew, i)
767 hq(ih) = hq(ih) + temp*xpt(knew, j)
776 temp = diff*zmat(knew, j)
777 IF (j < idz) temp = -temp
779 pq(k) = pq(k) + temp*zmat(k, j)
784 gq(i) = gq(i) + diff*bmat(knew, i)
785 gqsq = gqsq + gq(i)**2
786 xpt(knew, i) = xnew(i)
793 IF (ksave == 0 .AND. delta == rho)
THEN
794 IF (abs(ratio) > 1.0e-2_dp)
THEN
798 vlag(k) = fval(k) - fval(kopt)
804 sum = sum + bmat(k, i)*vlag(k)
806 gisq = gisq + sum*sum
814 IF (gqsq < 1.0e2_dp*gisq) itest = 0
825 w(j) = w(j) + vlag(k)*zmat(k, j)
827 IF (j < idz) w(j) = -w(j)
832 pq(k) = pq(k) + zmat(k, j)*w(j)
839 IF (f < fsave) kopt = knew
845 IF (f <= fsave + tenth*vquad)
GOTO 100
846 IF (ksave > 0)
GOTO 100
852 460 distsq = 4.0_dp*delta*delta
856 sum = sum + (xpt(k, j) - xopt(j))**2
858 IF (sum > distsq)
THEN
868 dstep = max(min(tenth*sqrt(distsq), half*delta), rho)
872 IF (ratio > zero)
GOTO 100
873 IF (max(delta, dnorm) > rho)
GOTO 100
878 490
IF (rho > rhoend)
THEN
881 IF (ratio <= 16.0_dp)
THEN
883 ELSE IF (ratio <= 250.0_dp)
THEN
884 rho = sqrt(ratio)*rhoend
888 delta = max(delta, rho)
895 IF (knew == -1)
GOTO 290
898 530
IF (fopt <= f)
THEN
900 x(i) = xbase(i) + xopt(i)
913 SUBROUTINE get_state()
958 END SUBROUTINE get_state
964 SUBROUTINE set_state()
1009 END SUBROUTINE set_state
1011 END SUBROUTINE newuob
1035 SUBROUTINE bigden(n, npt, xopt, xpt, bmat, zmat, idz, ndim, kopt, &
1036 knew, d, w, vlag, beta, s, wvec, prod)
1038 INTEGER,
INTENT(in) :: n, npt
1039 REAL(
dp),
DIMENSION(*),
INTENT(in) :: xopt
1040 REAL(
dp),
DIMENSION(npt, *),
INTENT(in) :: xpt
1041 INTEGER,
INTENT(in) :: ndim, idz
1042 REAL(
dp),
DIMENSION(npt, *),
INTENT(inout) :: zmat
1043 REAL(
dp),
DIMENSION(ndim, *),
INTENT(inout) :: bmat
1044 INTEGER,
INTENT(inout) :: kopt, knew
1045 REAL(
dp),
DIMENSION(*),
INTENT(inout) :: d, w, vlag
1046 REAL(
dp),
INTENT(inout) :: beta
1047 REAL(
dp),
DIMENSION(*),
INTENT(inout) :: s
1048 REAL(
dp),
DIMENSION(ndim, *),
INTENT(inout) :: wvec, prod
1050 REAL(
dp),
PARAMETER :: half = 0.5_dp, one = 1._dp, &
1051 quart = 0.25_dp, two = 2._dp, &
1054 INTEGER :: i, ip, isave, iterc, iu, j, jc, k, ksav, &
1056 REAL(
dp) :: alpha,
angle, dd, denmax, denold, densav, diff, ds, dstemp, dtest, ss, ssden, &
1057 sstemp, step, sum, sumold, tau, temp, tempa, tempb, tempc, xoptd, xopts, xoptsq
1058 REAL(
dp),
DIMENSION(9) :: den, denex, par
1094 temp = zmat(knew, j)
1095 IF (j < idz) temp = -temp
1097 w(n + k) = w(n + k) + temp*zmat(k, j)
1113 s(i) = xpt(knew, i) - xopt(i)
1116 xoptsq = xoptsq + xopt(i)**2
1118 IF (ds*ds > 0.99_dp*dd*ss)
THEN
1126 diff = xpt(k, i) - xopt(i)
1127 dstemp = dstemp + d(i)*diff
1128 sstemp = sstemp + diff*diff
1130 IF (dstemp*dstemp/sstemp < dtest)
THEN
1132 dtest = dstemp*dstemp/sstemp
1139 s(i) = xpt(ksav, i) - xopt(i)
1142 ssden = dd*ss - ds*ds
1151 temp = one/sqrt(ssden)
1155 s(i) = temp*(dd*s(i) - ds*d(i))
1156 xoptd = xoptd + xopt(i)*d(i)
1157 xopts = xopts + xopt(i)*s(i)
1162 tempa = half*xoptd*xoptd
1163 tempb = half*xopts*xopts
1164 den(1) = dd*(xoptsq + half*dd) + tempa + tempb
1165 den(2) = two*xoptd*dd
1166 den(3) = two*xopts*dd
1167 den(4) = tempa - tempb
1168 den(5) = xoptd*xopts
1180 tempa = tempa + xpt(k, i)*d(i)
1181 tempb = tempb + xpt(k, i)*s(i)
1182 tempc = tempc + xpt(k, i)*xopt(i)
1184 wvec(k, 1) = quart*(tempa*tempa + tempb*tempb)
1185 wvec(k, 2) = tempa*tempc
1186 wvec(k, 3) = tempb*tempc
1187 wvec(k, 4) = quart*(tempa*tempa - tempb*tempb)
1188 wvec(k, 5) = half*tempa*tempb
1203 IF (jc == 2 .OR. jc == 3) nw = ndim
1210 sum = sum + zmat(k, j)*wvec(k, jc)
1212 IF (j < idz) sum = -sum
1214 prod(k, jc) = prod(k, jc) + sum*zmat(k, j)
1217 IF (nw == ndim)
THEN
1221 sum = sum + bmat(k, j)*wvec(npt + j, jc)
1223 prod(k, jc) = prod(k, jc) + sum
1229 sum = sum + bmat(i, j)*wvec(i, jc)
1231 prod(npt + j, jc) = sum
1240 par(i) = half*prod(k, i)*wvec(k, i)
1243 den(1) = den(1) - par(1) - sum
1244 tempa = prod(k, 1)*wvec(k, 2) + prod(k, 2)*wvec(k, 1)
1245 tempb = prod(k, 2)*wvec(k, 4) + prod(k, 4)*wvec(k, 2)
1246 tempc = prod(k, 3)*wvec(k, 5) + prod(k, 5)*wvec(k, 3)
1247 den(2) = den(2) - tempa - half*(tempb + tempc)
1248 den(6) = den(6) - half*(tempb - tempc)
1249 tempa = prod(k, 1)*wvec(k, 3) + prod(k, 3)*wvec(k, 1)
1250 tempb = prod(k, 2)*wvec(k, 5) + prod(k, 5)*wvec(k, 2)
1251 tempc = prod(k, 3)*wvec(k, 4) + prod(k, 4)*wvec(k, 3)
1252 den(3) = den(3) - tempa - half*(tempb - tempc)
1253 den(7) = den(7) - half*(tempb + tempc)
1254 tempa = prod(k, 1)*wvec(k, 4) + prod(k, 4)*wvec(k, 1)
1255 den(4) = den(4) - tempa - par(2) + par(3)
1256 tempa = prod(k, 1)*wvec(k, 5) + prod(k, 5)*wvec(k, 1)
1257 tempb = prod(k, 2)*wvec(k, 3) + prod(k, 3)*wvec(k, 2)
1258 den(5) = den(5) - tempa - half*tempb
1259 den(8) = den(8) - par(4) + par(5)
1260 tempa = prod(k, 4)*wvec(k, 5) + prod(k, 5)*wvec(k, 4)
1261 den(9) = den(9) - half*tempa
1268 par(i) = half*prod(knew, i)**2
1271 denex(1) = alpha*den(1) + par(1) + sum
1272 tempa = two*prod(knew, 1)*prod(knew, 2)
1273 tempb = prod(knew, 2)*prod(knew, 4)
1274 tempc = prod(knew, 3)*prod(knew, 5)
1275 denex(2) = alpha*den(2) + tempa + tempb + tempc
1276 denex(6) = alpha*den(6) + tempb - tempc
1277 tempa = two*prod(knew, 1)*prod(knew, 3)
1278 tempb = prod(knew, 2)*prod(knew, 5)
1279 tempc = prod(knew, 3)*prod(knew, 4)
1280 denex(3) = alpha*den(3) + tempa + tempb - tempc
1281 denex(7) = alpha*den(7) + tempb + tempc
1282 tempa = two*prod(knew, 1)*prod(knew, 4)
1283 denex(4) = alpha*den(4) + tempa + par(2) - par(3)
1284 tempa = two*prod(knew, 1)*prod(knew, 5)
1285 denex(5) = alpha*den(5) + tempa + prod(knew, 2)*prod(knew, 3)
1286 denex(8) = alpha*den(8) + par(4) - par(5)
1287 denex(9) = alpha*den(9) + prod(knew, 4)*prod(knew, 5)
1291 sum = denex(1) + denex(2) + denex(4) + denex(6) + denex(8)
1303 par(j) = par(2)*par(j - 2) - par(3)*par(j - 1)
1304 par(j + 1) = par(2)*par(j - 1) + par(3)*par(j - 2)
1309 sum = sum + denex(j)*par(j)
1311 IF (abs(sum) > abs(denmax))
THEN
1315 ELSE IF (i == isave + 1)
THEN
1319 IF (isave == 0) tempa = sum
1320 IF (isave == iu) tempb = denold
1322 IF (tempa /= tempb)
THEN
1323 tempa = tempa - denmax
1324 tempb = tempb - denmax
1325 step = half*(tempa - tempb)/(tempa + tempb)
1327 angle = temp*(real(isave,
dp) + step)
1335 par(j) = par(2)*par(j - 2) - par(3)*par(j - 1)
1336 par(j + 1) = par(2)*par(j - 1) + par(3)*par(j - 2)
1341 beta = beta + den(j)*par(j)
1342 denmax = denmax + denex(j)*par(j)
1347 vlag(k) = vlag(k) + prod(k, j)*par(j)
1355 d(i) = par(2)*d(i) + par(3)*s(i)
1356 w(i) = xopt(i) + d(i)
1358 tempa = tempa + d(i)*w(i)
1359 tempb = tempb + w(i)*w(i)
1361 IF (iterc >= n)
EXIT mainloop
1362 IF (iterc >= 1) densav = max(densav, denold)
1363 IF (abs(denmax) <= 1.1_dp*abs(densav))
EXIT mainloop
1370 temp = tempa*xopt(i) + tempb*d(i) - vlag(npt + i)
1371 s(i) = tau*bmat(knew, i) + alpha*temp
1376 sum = sum + xpt(k, j)*w(j)
1378 temp = (tau*w(n + k) - alpha*vlag(k))*sum
1380 s(i) = s(i) + temp*xpt(k, i)
1389 ssden = dd*ss - ds*ds
1390 IF (ssden < 1.0e-8_dp*dd*ss)
EXIT mainloop
1398 w(k) = w(k) + wvec(k, j)*par(j)
1401 vlag(kopt) = vlag(kopt) + one
1403 END SUBROUTINE bigden
1426 SUBROUTINE biglag(n, npt, xopt, xpt, bmat, zmat, idz, ndim, knew, &
1427 delta, d, alpha, hcol, gc, gd, s, w)
1428 INTEGER,
INTENT(in) :: n, npt
1429 REAL(
dp),
DIMENSION(*),
INTENT(in) :: xopt
1430 REAL(
dp),
DIMENSION(npt, *),
INTENT(in) :: xpt
1431 INTEGER,
INTENT(in) :: ndim, idz
1432 REAL(
dp),
DIMENSION(npt, *),
INTENT(inout) :: zmat
1433 REAL(
dp),
DIMENSION(ndim, *),
INTENT(inout) :: bmat
1434 INTEGER,
INTENT(inout) :: knew
1435 REAL(
dp),
INTENT(inout) :: delta
1436 REAL(
dp),
DIMENSION(*),
INTENT(inout) :: d
1437 REAL(
dp),
INTENT(inout) :: alpha
1438 REAL(
dp),
DIMENSION(*),
INTENT(inout) :: hcol, gc, gd, s, w
1440 REAL(
dp),
PARAMETER :: half = 0.5_dp, one = 1._dp, zero = 0._dp
1442 INTEGER :: i, isave, iterc, iu, j, k, nptm
1443 REAL(
dp) ::
angle, cf1, cf2, cf3, cf4, cf5, cth, dd, &
1444 delsq, denom, dhd, gg, scale, sp, ss, &
1445 step, sth, sum, tau, taubeg, taumax, &
1446 tauold, temp, tempa, tempb
1479 temp = zmat(knew, j)
1480 IF (j < idz) temp = -temp
1482 hcol(k) = hcol(k) + temp*zmat(k, j)
1492 d(i) = xpt(knew, i) - xopt(i)
1493 gc(i) = bmat(knew, i)
1501 temp = temp + xpt(k, j)*xopt(j)
1502 sum = sum + xpt(k, j)*d(j)
1507 gc(i) = gc(i) + temp*xpt(k, i)
1508 gd(i) = gd(i) + sum*xpt(k, i)
1520 sp = sp + d(i)*gc(i)
1521 dhd = dhd + d(i)*gd(i)
1523 scale = delta/sqrt(dd)
1524 IF (sp*dhd < zero) scale = -scale
1526 IF (sp*sp > 0.99_dp*dd*gg) temp = one
1527 tau = scale*(abs(sp) + half*scale*abs(dhd))
1528 IF (gg*delsq < 0.01_dp*tau*tau) temp = one
1532 s(i) = gc(i) + temp*gd(i)
1549 temp = dd*ss - sp*sp
1550 IF (temp <= 1.0e-8_dp*dd*ss)
EXIT mainloop
1553 s(i) = (dd*s(i) - sp*d(i))/denom
1563 sum = sum + xpt(k, j)*s(j)
1567 w(i) = w(i) + sum*xpt(k, i)
1576 cf1 = cf1 + s(i)*w(i)
1577 cf2 = cf2 + d(i)*gc(i)
1578 cf3 = cf3 + s(i)*gc(i)
1579 cf4 = cf4 + d(i)*gd(i)
1580 cf5 = cf5 + s(i)*gd(i)
1583 cf4 = half*cf4 - cf1
1587 taubeg = cf1 + cf2 + cf4
1597 tau = cf1 + (cf2 + cf4*cth)*cth + (cf3 + cf5*cth)*sth
1598 IF (abs(tau) > abs(taumax))
THEN
1602 ELSE IF (i == isave + 1)
THEN
1607 IF (isave == 0) tempa = tau
1608 IF (isave == iu) tempb = taubeg
1610 IF (tempa /= tempb)
THEN
1611 tempa = tempa - taumax
1612 tempb = tempb - taumax
1613 step = half*(tempa - tempb)/(tempa + tempb)
1615 angle = temp*(real(isave,
dp) + step)
1621 tau = cf1 + (cf2 + cf4*cth)*cth + (cf3 + cf5*cth)*sth
1623 d(i) = cth*d(i) + sth*s(i)
1624 gd(i) = cth*gd(i) + sth*w(i)
1625 s(i) = gc(i) + gd(i)
1627 IF (abs(tau) <= 1.1_dp*abs(taubeg))
EXIT mainloop
1628 IF (iterc >= n)
EXIT mainloop
1631 END SUBROUTINE biglag
1651 SUBROUTINE trsapp(n, npt, xopt, xpt, gq, hq, pq, delta, step, d, g, hd, hs, crvmin)
1653 INTEGER,
INTENT(IN) :: n, npt
1654 REAL(
dp),
DIMENSION(*),
INTENT(IN) :: xopt
1655 REAL(
dp),
DIMENSION(npt, *), &
1657 REAL(
dp),
DIMENSION(*),
INTENT(INOUT) :: gq, hq, pq
1658 REAL(
dp),
INTENT(IN) :: delta
1659 REAL(
dp),
DIMENSION(*),
INTENT(INOUT) :: step, d, g, hd, hs
1660 REAL(
dp),
INTENT(INOUT) :: crvmin
1662 REAL(
dp),
PARAMETER :: half = 0.5_dp, zero = 0.0_dp
1664 INTEGER :: i, isave, iterc, itermax, &
1666 LOGICAL :: jump1, jump2
1667 REAL(
dp) :: alpha,
angle, angtest, bstep, cf, cth, dd, delsq, dg, dhd, &
1668 dhs, ds, gg, ggbeg, ggsav, qadd, qbeg, qmin, qnew, qred, qsav, ratio, &
1669 reduc, sg, sgk, shs, ss, sth, temp, tempa, tempb
1707 g(i) = gq(i) + hd(i)
1712 IF (dd == zero)
RETURN
1723 IF (.NOT. jump2)
THEN
1724 IF (.NOT. jump1)
THEN
1727 bstep = temp/(ds + sqrt(ds*ds + dd*temp))
1731 IF (iterc <= itersw)
THEN
1734 dhd = dhd + d(j)*hd(j)
1740 IF (dhd > zero)
THEN
1742 IF (iterc == 1) crvmin = temp
1743 crvmin = min(crvmin, temp)
1744 alpha = min(alpha, gg/dhd)
1746 qadd = alpha*(gg - half*alpha*dhd)
1754 step(i) = step(i) + alpha*d(i)
1755 hs(i) = hs(i) + alpha*hd(i)
1756 gg = gg + (g(i) + hs(i))**2
1761 IF (alpha < bstep)
THEN
1762 IF (qadd <= 0.01_dp*qred)
EXIT mainloop
1763 IF (gg <= 1.0e-4_dp*ggbeg)
EXIT mainloop
1764 IF (iterc == itermax)
EXIT mainloop
1770 d(i) = temp*d(i) - g(i) - hs(i)
1772 ds = ds + d(i)*step(i)
1773 ss = ss + step(i)**2
1775 IF (ds <= zero)
EXIT mainloop
1776 IF (ss < delsq) cycle mainloop
1781 IF (gg <= 1.0e-4_dp*ggbeg)
EXIT mainloop
1794 sg = sg + step(i)*g(i)
1795 shs = shs + step(i)*hs(i)
1798 angtest = sgk/sqrt(gg*delsq)
1799 IF (angtest <= -0.99_dp)
EXIT mainloop
1805 temp = sqrt(delsq*gg - sgk*sgk)
1809 d(i) = tempa*(g(i) + hs(i)) - tempb*step(i)
1812 IF (iterc <= itersw)
THEN
1822 dhd = dhd + hd(i)*d(i)
1823 dhs = dhs + hd(i)*step(i)
1828 cf = half*(shs - dhd)
1839 qnew = (sg + cf*cth)*cth + (dg + dhs*cth)*sth
1840 IF (qnew < qmin)
THEN
1844 ELSE IF (i == isave + 1)
THEN
1849 IF (isave == zero) tempa = qnew
1850 IF (isave == iu) tempb = qbeg
1852 IF (tempa /= tempb)
THEN
1853 tempa = tempa - qmin
1854 tempb = tempb - qmin
1855 angle = half*(tempa - tempb)/(tempa + tempb)
1863 reduc = qbeg - (sg + cf*cth)*cth - (dg + dhs*cth)*sth
1866 step(i) = cth*step(i) + sth*d(i)
1867 hs(i) = cth*hs(i) + sth*hd(i)
1868 gg = gg + (g(i) + hs(i))**2
1872 IF (iterc < itermax .AND. ratio > 0.01_dp)
THEN
1878 IF (gg <= 1.0e-4_dp*ggbeg)
EXIT mainloop
1888 INTEGER :: i, ih, j, k
1896 temp = temp + xpt(k, j)*d(j)
1900 hd(i) = hd(i) + temp*xpt(k, i)
1907 IF (i < j) hd(j) = hd(j) + hq(ih)*d(i)
1908 hd(i) = hd(i) + hq(ih)*d(j)
1911 END SUBROUTINE updatehd
1913 END SUBROUTINE trsapp
1929 SUBROUTINE update(n, npt, bmat, zmat, idz, ndim, vlag, beta, knew, w)
1931 INTEGER,
INTENT(IN) :: n, npt, ndim
1932 INTEGER,
INTENT(INOUT) :: idz
1933 REAL(
dp),
DIMENSION(npt, *),
INTENT(INOUT) :: zmat
1934 REAL(
dp),
DIMENSION(ndim, *),
INTENT(INOUT) :: bmat
1935 REAL(
dp),
DIMENSION(*),
INTENT(INOUT) :: vlag
1936 REAL(
dp),
INTENT(INOUT) :: beta
1937 INTEGER,
INTENT(INOUT) :: knew
1938 REAL(
dp),
DIMENSION(*),
INTENT(INOUT) :: w
1940 REAL(
dp),
PARAMETER :: one = 1.0_dp, zero = 0.0_dp
1942 INTEGER :: i, iflag, j, ja, jb, jl, jp, nptm
1943 REAL(
dp) :: alpha, denom, scala, scalb, tau, tausq, &
1961 ELSE IF (zmat(knew, j) /= zero)
THEN
1962 temp = sqrt(zmat(knew, jl)**2 + zmat(knew, j)**2)
1963 tempa = zmat(knew, jl)/temp
1964 tempb = zmat(knew, j)/temp
1966 temp = tempa*zmat(i, jl) + tempb*zmat(i, j)
1967 zmat(i, j) = tempa*zmat(i, j) - tempb*zmat(i, jl)
1970 zmat(knew, j) = zero
1977 tempa = zmat(knew, 1)
1978 IF (idz >= 2) tempa = -tempa
1979 IF (jl > 1) tempb = zmat(knew, jl)
1981 w(i) = tempa*zmat(i, 1)
1982 IF (jl > 1) w(i) = w(i) + tempb*zmat(i, jl)
1987 denom = alpha*beta + tausq
1988 vlag(knew) = vlag(knew) - one
1996 temp = sqrt(abs(denom))
2000 zmat(i, 1) = tempa*zmat(i, 1) - tempb*vlag(i)
2002 IF (idz == 1 .AND. temp < zero) idz = 2
2003 IF (idz >= 2 .AND. temp >= zero) iflag = 1
2009 IF (beta >= zero) ja = jl
2011 temp = zmat(knew, jb)/denom
2014 temp = zmat(knew, ja)
2015 scala = one/sqrt(abs(beta)*temp*temp + tausq)
2016 scalb = scala*sqrt(abs(denom))
2018 zmat(i, ja) = scala*(tau*zmat(i, ja) - temp*vlag(i))
2019 zmat(i, jb) = scalb*(zmat(i, jb) - tempa*w(i) - tempb*vlag(i))
2021 IF (denom <= zero)
THEN
2022 IF (beta < zero) idz = idz + 1
2023 IF (beta >= zero) iflag = 1
2030 IF (iflag == 1)
THEN
2034 zmat(i, 1) = zmat(i, idz)
2043 w(jp) = bmat(knew, j)
2044 tempa = (alpha*vlag(jp) - tau*w(jp))/denom
2045 tempb = (-beta*w(jp) - tau*vlag(jp))/denom
2047 bmat(i, j) = bmat(i, j) + tempa*vlag(i) + tempb*w(i)
2048 IF (i > npt) bmat(jp, i - npt) = bmat(i, j)
2052 END SUBROUTINE update
pure real(kind=dp) function angle(a, b)
Calculation of the angle between the vectors a and b. The angle is returned in radians.
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public twopi
subroutine, public powell_optimize(n, x, optstate)
...