28 #include "../base/base_uses.f90"
34 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'cp_lbfgs'
183 SUBROUTINE setulb(n, m, x, lower_bound, upper_bound, nbd, f, g, factr, pgtol, wa, iwa, &
184 task, iprint, csave, lsave, isave, dsave, trust_radius, spgr)
186 INTEGER,
INTENT(in) :: n, m
187 REAL(kind=
dp),
INTENT(inout) :: x(n)
188 REAL(kind=
dp) :: lower_bound(n), upper_bound(n)
190 REAL(kind=
dp) :: f, g(n)
191 REAL(kind=
dp),
INTENT(in) :: factr, pgtol
192 REAL(kind=
dp) :: wa(2*m*n + 5*n + 11*m*m + 8*m)
194 CHARACTER(LEN=60) :: task
196 CHARACTER(LEN=60) :: csave
199 REAL(kind=
dp) :: dsave(29)
200 REAL(kind=
dp),
INTENT(in) :: trust_radius
201 TYPE(spgr_type),
OPTIONAL,
POINTER :: spgr
203 INTEGER :: i, ld, lr, lsnd, lss, lsy, lt, lwa, lwn, &
204 lws, lwt, lwy, lxp, lz
222 IF (task ==
'START')
THEN
230 isave(5) = isave(4) + isave(1)
232 isave(6) = isave(5) + isave(1)
234 isave(7) = isave(6) + isave(2)
236 isave(8) = isave(7) + isave(2)
238 isave(9) = isave(8) + isave(2)
240 isave(10) = isave(9) + isave(3)
242 isave(11) = isave(10) + isave(3)
244 isave(12) = isave(11) + n
246 isave(13) = isave(12) + n
248 isave(14) = isave(13) + n
250 isave(15) = isave(14) + n
252 isave(16) = isave(15) + n
272 IF (trust_radius >= 0)
THEN
274 lower_bound(i) = x(i) - trust_radius
275 upper_bound(i) = x(i) + trust_radius
281 CALL mainlb(n, m, x, lower_bound, upper_bound, nbd, f, g, factr, pgtol, &
282 wa(lws), wa(lwy), wa(lsy), wa(lss), wa(lwt), &
283 wa(lwn), wa(lsnd), wa(lz), wa(lr), wa(ld), wa(lt), wa(lxp), &
285 iwa(1), iwa(n + 1), iwa(2*n + 1), task, iprint, &
286 csave, lsave, isave(22), dsave, spgr=spgr)
390 SUBROUTINE mainlb(n, m, x, lower_bound, upper_bound, nbd, f, g, factr, pgtol, ws, wy, &
391 sy, ss, wt, wn, snd, z, r, d, t, xp, wa, &
392 index, iwhere, indx2, task, &
393 iprint, csave, lsave, isave, dsave, spgr)
394 INTEGER,
INTENT(in) :: n, m
395 REAL(kind=
dp),
INTENT(inout) :: x(n)
396 REAL(kind=
dp),
INTENT(in) :: lower_bound(n), upper_bound(n)
398 REAL(kind=
dp) :: f, g(n), factr, pgtol, ws(n, m), wy(n, m), sy(m, m), ss(m, m), wt(m, m), &
399 wn(2*m, 2*m), snd(2*m, 2*m), z(n), r(n), d(n), t(n), xp(n), wa(8*m)
400 INTEGER :: index(n), iwhere(n), indx2(n)
401 CHARACTER(LEN=60) :: task
403 CHARACTER(LEN=60) :: csave
406 REAL(kind=
dp) :: dsave(29)
407 TYPE(spgr_type),
OPTIONAL,
POINTER :: spgr
409 REAL(kind=
dp),
PARAMETER :: one = 1.0_dp, zero = 0.0_dp
411 CHARACTER(LEN=3) :: word
412 INTEGER :: col, head, i, iback, ifun, ileave, info, &
413 itail, iter, itfile, iupdat, iword, k, &
414 nact, nenter, nfgv, nfree, nintol, &
416 LOGICAL :: boxed, constrained, first, &
417 keep_space_group, updatd, wrk, &
419 REAL(kind=
dp) :: cachyt, cpu1, cpu2, ddot, ddum, dnorm, dr, dtd, epsmch, fold, g_inf_norm, &
420 gd, gdold, lnscht, rr, sbtime, step_max, stp, theta, time, time1, time2, tol, xstep
442 keep_space_group = .false.
443 IF (
PRESENT(spgr))
THEN
444 IF (
ASSOCIATED(spgr)) keep_space_group = spgr%keep_space_group
447 IF (task ==
'START')
THEN
449 epsmch = epsilon(one)
500 IF (iprint >= 1)
THEN
502 CALL open_file(file_name=
'iterate.dat', unit_number=itfile, file_action=
'WRITE', file_status=
'UNKNOWN')
507 CALL errclb(n, m, factr, lower_bound, upper_bound, nbd, task, info, k)
508 IF (task(1:5) ==
'ERROR')
THEN
509 CALL prn3lb(n, x, f, task, iprint, info, itfile, &
510 iter, nfgv, nintol, nskip, nact, g_inf_norm, &
511 zero, nseg, word, iback, stp, xstep, k, &
512 cachyt, sbtime, lnscht)
516 CALL prn1lb(n, m, lower_bound, upper_bound, x, iprint, itfile, epsmch)
520 CALL active(n, lower_bound, upper_bound, nbd, x, iwhere, iprint, x_projected, constrained, boxed)
522 IF (keep_space_group)
THEN
529 CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
530 iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
531 cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
535 IF (keep_space_group)
THEN
542 x_projected = lsave(1)
543 constrained = lsave(2)
578 g_inf_norm = dsave(13)
586 IF (task(1:4) ==
'STOP')
THEN
587 IF (task(7:9) ==
'CPU')
THEN
589 CALL dcopy(n, t, 1, x, 1)
590 CALL dcopy(n, r, 1, g, 1)
595 CALL prn3lb(n, x, f, task, iprint, info, itfile, &
596 iter, nfgv, nintol, nskip, nact, g_inf_norm, &
597 time, nseg, word, iback, stp, xstep, k, &
598 cachyt, sbtime, lnscht)
599 CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
600 iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
601 cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
606 IF (.NOT. (task(1:5) ==
'FG_LN' .OR. task(1:5) ==
'NEW_X'))
THEN
613 CALL projgr(n, lower_bound, upper_bound, nbd, x, g, g_inf_norm)
615 IF (iprint >= 1)
THEN
616 WRITE (*, 1002) iter, f, g_inf_norm
617 WRITE (itfile, 1003) iter, nfgv, g_inf_norm, f
619 IF (g_inf_norm <= pgtol)
THEN
621 task =
'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
624 CALL prn3lb(n, x, f, task, iprint, info, itfile, &
625 iter, nfgv, nintol, nskip, nact, g_inf_norm, &
626 time, nseg, word, iback, stp, xstep, k, &
627 cachyt, sbtime, lnscht)
628 CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
629 iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
630 cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
637 IF (.NOT. first .OR. .NOT. (task(1:5) ==
'FG_LN' .OR. task(1:5) ==
'NEW_X'))
THEN
638 IF (iprint >= 99)
WRITE (*, 1001) iter + 1
641 IF (.NOT. constrained .AND. col > 0)
THEN
643 CALL dcopy(n, x, 1, z, 1)
655 CALL cauchy(n, x, lower_bound, upper_bound, nbd, g, indx2, iwhere, t, d, z, &
656 m, wy, ws, sy, wt, theta, col, head, &
657 wa(1), wa(2*m + 1), wa(4*m + 1), wa(6*m + 1), nseg, &
658 iprint, g_inf_norm, info, epsmch)
660 IF (keep_space_group)
THEN
665 IF (iprint >= 1)
WRITE (*, 1005)
673 cachyt = cachyt + cpu2 - cpu1
678 cachyt = cachyt + cpu2 - cpu1
679 nintol = nintol + nseg
684 CALL freev(n, nfree, index, nenter, ileave, indx2, &
685 iwhere, wrk, updatd, constrained, iprint, iter)
693 IF (.NOT. (nfree == 0 .OR. col == 0))
THEN
709 IF (wrk)
CALL formk(n, nfree, index, nenter, ileave, indx2, iupdat, &
710 updatd, wn, snd, m, ws, wy, sy, theta, col, head, info)
714 IF (iprint >= 1)
WRITE (*, 1006)
722 sbtime = sbtime + cpu2 - cpu1
729 CALL cmprlb(n, m, x, g, ws, wy, sy, wt, z, r, wa, index, &
730 theta, col, head, nfree, constrained, info)
732 IF (keep_space_group)
THEN
739 CALL subsm(n, m, nfree, index, lower_bound, upper_bound, nbd, z, r, xp, ws, wy, &
740 theta, x, g, col, head, iword, wa, wn, iprint, info)
742 IF (keep_space_group)
THEN
750 IF (iprint >= 1)
WRITE (*, 1005)
758 sbtime = sbtime + cpu2 - cpu1
764 sbtime = sbtime + cpu2 - cpu1
775 IF (keep_space_group)
THEN
784 IF (.NOT. first .OR. .NOT. (task(1:5) ==
'NEW_X'))
THEN
786 IF (keep_space_group)
THEN
793 CALL lnsrlb(n, lower_bound, upper_bound, nbd, x, f, fold, gd, gdold, g, d, r, t, z, stp, dnorm, &
794 dtd, xstep, step_max, iter, ifun, iback, nfgv, info, task, &
795 boxed, constrained, csave, isave(22), dsave(17))
797 IF (keep_space_group)
THEN
801 IF (info /= 0 .OR. iback >= 20)
THEN
803 CALL dcopy(n, t, 1, x, 1)
804 CALL dcopy(n, r, 1, g, 1)
815 task =
'ABNORMAL_TERMINATION_IN_LNSRCH'
819 CALL prn3lb(n, x, f, task, iprint, info, itfile, &
820 iter, nfgv, nintol, nskip, nact, g_inf_norm, &
821 time, nseg, word, iback, stp, xstep, k, &
822 cachyt, sbtime, lnscht)
823 CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
824 iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
825 cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
829 IF (iprint >= 1)
WRITE (*, 1008)
830 IF (info == 0) nfgv = nfgv - 1
837 task =
'RESTART_FROM_LNSRCH'
839 lnscht = lnscht + cpu2 - cpu1
843 ELSE IF (task(1:5) ==
'FG_LN')
THEN
845 CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
846 iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
847 cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
852 lnscht = lnscht + cpu2 - cpu1
857 CALL projgr(n, lower_bound, upper_bound, nbd, x, g, g_inf_norm)
861 CALL prn2lb(n, x, f, g, iprint, itfile, iter, nfgv, nact, &
862 g_inf_norm, nseg, word, iword, iback, stp, xstep)
863 CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
864 iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
865 cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
872 IF (g_inf_norm <= pgtol)
THEN
874 task =
'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
877 CALL prn3lb(n, x, f, task, iprint, info, itfile, &
878 iter, nfgv, nintol, nskip, nact, g_inf_norm, &
879 time, nseg, word, iback, stp, xstep, k, &
880 cachyt, sbtime, lnscht)
881 CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
882 iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
883 cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
887 ddum = max(abs(fold), abs(f), one)
888 IF ((fold - f) <= tol*ddum)
THEN
890 task =
'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
891 IF (iback >= 10) info = -5
895 CALL prn3lb(n, x, f, task, iprint, info, itfile, &
896 iter, nfgv, nintol, nskip, nact, g_inf_norm, &
897 time, nseg, word, iback, stp, xstep, k, &
898 cachyt, sbtime, lnscht)
899 CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
900 iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
901 cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
906 IF (keep_space_group)
THEN
913 rr = ddot(n, r, 1, r, 1)
918 dr = (gd - gdold)*stp
919 CALL dscal(n, stp, d, 1)
923 IF (dr <= epsmch*ddum)
THEN
927 IF (iprint >= 1)
WRITE (*, 1004) dr, ddum
943 CALL matupd(n, m, ws, wy, sy, ss, d, r, itail, &
944 iupdat, col, head, theta, rr, dr, stp, dtd)
951 CALL formt(m, wt, sy, ss, col, theta, info)
956 IF (iprint >= 1)
WRITE (*, 1007)
973 1001
FORMAT(//,
'ITERATION ', i5)
975 (/,
'At iterate', i5, 4x,
'f= ', 1p, d12.5, 4x,
'|proj g|= ', 1p, d12.5)
976 1003
FORMAT(2(1x, i4), 5x,
'-', 5x,
'-', 3x,
'-', 5x,
'-', 5x,
'-', 8x,
'-', 3x, &
978 1004
FORMAT(
' ys=', 1p, e10.3,
' -gs=', 1p, e10.3,
' BFGS update SKIPPED')
980 ' Singular triangular system detected;', /, &
981 ' refresh the lbfgs memory and restart the iteration.')
983 ' Nonpositive definiteness in Cholesky factorization in formk;', /, &
984 ' refresh the lbfgs memory and restart the iteration.')
986 ' Nonpositive definiteness in Cholesky factorization in formt;', /, &
987 ' refresh the lbfgs memory and restart the iteration.')
989 ' Bad direction in the line search;', /, &
990 ' refresh the lbfgs memory and restart the iteration.')
994 END SUBROUTINE mainlb
1018 SUBROUTINE active(n, lower_bound, upper_bound, nbd, x, iwhere, iprint, &
1019 x_projected, constrained, boxed)
1021 INTEGER,
INTENT(in) :: n
1022 REAL(kind=
dp),
INTENT(in) :: lower_bound(n), upper_bound(n)
1024 REAL(kind=
dp) :: x(n)
1025 INTEGER,
INTENT(out) :: iwhere(n)
1027 LOGICAL :: x_projected, constrained, boxed
1029 REAL(kind=
dp),
PARAMETER :: zero = 0.0_dp
1037 x_projected = .false.
1038 constrained = .false.
1044 IF (nbd(i) > 0)
THEN
1045 IF (nbd(i) <= 2 .AND. x(i) <= lower_bound(i))
THEN
1046 IF (x(i) < lower_bound(i))
THEN
1047 x_projected = .true.
1048 x(i) = lower_bound(i)
1051 ELSE IF (nbd(i) >= 2 .AND. x(i) >= upper_bound(i))
THEN
1052 IF (x(i) > upper_bound(i))
THEN
1053 x_projected = .true.
1054 x(i) = upper_bound(i)
1064 IF (nbd(i) /= 2) boxed = .false.
1065 IF (nbd(i) == 0)
THEN
1071 constrained = .true.
1072 IF (nbd(i) == 2 .AND. upper_bound(i) - lower_bound(i) <= zero)
THEN
1081 IF (iprint >= 0)
THEN
1082 IF (x_projected)
WRITE (*, *) &
1083 &
'The initial X is infeasible. Restart with its projection.'
1084 IF (.NOT. constrained) &
1085 WRITE (*, *)
'This problem is unconstrained.'
1088 IF (iprint > 0)
WRITE (*, 1001) nbdd
1090 1001
FORMAT(/,
'At X0 ', i9,
' variables are exactly at the bounds')
1094 END SUBROUTINE active
1118 SUBROUTINE bmv(m, sy, wt, col, v, p, info)
1121 REAL(kind=
dp) :: sy(m, m), wt(m, m)
1123 REAL(kind=
dp),
INTENT(in) :: v(2*col)
1124 REAL(kind=
dp),
INTENT(out) :: p(2*col)
1125 INTEGER,
INTENT(out) :: info
1128 REAL(kind=
dp) :: sum
1130 IF (col == 0)
RETURN
1136 p(col + 1) = v(col + 1)
1141 sum = sum + sy(i, k)*v(k)/sy(k, k)
1146 CALL dtrsl(wt, m, col, p(col + 1), 11, info)
1147 IF (info /= 0)
RETURN
1151 p(i) = v(i)/sqrt(sy(i, i))
1158 CALL dtrsl(wt, m, col, p(col + 1), 01, info)
1159 IF (info /= 0)
RETURN
1164 p(i) = -p(i)/sqrt(sy(i, i))
1169 sum = sum + sy(k, i)*p(col + k)/sy(i, i)
1257 SUBROUTINE cauchy(n, x, lower_bound, upper_bound, nbd, g, iorder, iwhere, t, d, xcp, &
1258 m, wy, ws, sy, wt, theta, col, head, p, c, wbp, &
1259 v, nseg, iprint, g_inf_norm, info, epsmch)
1260 INTEGER,
INTENT(in) :: n
1261 REAL(kind=
dp),
INTENT(in) :: x(n), lower_bound(n), upper_bound(n)
1262 INTEGER,
INTENT(in) :: nbd(n)
1263 REAL(kind=
dp),
INTENT(in) :: g(n)
1264 INTEGER :: iorder(n)
1265 INTEGER,
INTENT(inout) :: iwhere(n)
1266 REAL(kind=
dp) :: t(n), d(n), xcp(n)
1267 INTEGER,
INTENT(in) :: m
1268 REAL(kind=
dp),
INTENT(in) :: sy(m, m), wt(m, m), theta
1269 INTEGER,
INTENT(in) :: col
1270 REAL(kind=
dp),
INTENT(in) :: ws(n, col), wy(n, col)
1271 INTEGER,
INTENT(in) :: head
1272 REAL(kind=
dp) :: p(2*m), c(2*m), wbp(2*m), v(2*m)
1273 INTEGER :: nseg, iprint
1274 REAL(kind=
dp),
INTENT(in) :: g_inf_norm
1275 INTEGER,
INTENT(inout) :: info
1276 REAL(kind=
dp) :: epsmch
1278 REAL(kind=
dp),
PARAMETER :: one = 1.0_dp, zero = 0.0_dp
1280 INTEGER :: col2, i, ibkmin, ibp, iter, j, nbreak, &
1281 nfree, nleft, pointr
1282 LOGICAL :: bnded, xlower, xupper
1283 REAL(kind=
dp) :: bkmin, ddot, dibp, dibp2, dt, dtm, f1, &
1284 f2, f2_org, neggi, tj, tj0, tl, tsum, &
1285 tu, wmc, wmp, wmw, zibp
1306 IF (g_inf_norm <= zero)
THEN
1307 IF (iprint >= 0)
WRITE (*, *)
'Subgnorm = 0. GCP = X.'
1308 CALL dcopy(n, x, 1, xcp, 1)
1318 IF (iprint >= 99)
WRITE (*, 3010)
1332 IF (iwhere(i) /= 3 .AND. iwhere(i) /= -1)
THEN
1335 IF (nbd(i) <= 2) tl = x(i) - lower_bound(i)
1336 IF (nbd(i) >= 2) tu = upper_bound(i) - x(i)
1340 xlower = nbd(i) <= 2 .AND. tl <= zero
1341 xupper = nbd(i) >= 2 .AND. tu <= zero
1346 IF (neggi <= zero) iwhere(i) = 1
1347 ELSE IF (xupper)
THEN
1348 IF (neggi >= zero) iwhere(i) = 2
1350 IF (abs(neggi) <= zero) iwhere(i) = -3
1354 IF (iwhere(i) /= 0 .AND. iwhere(i) /= -1)
THEN
1358 f1 = f1 - neggi*neggi
1361 p(j) = p(j) + wy(i, pointr)*neggi
1362 p(col + j) = p(col + j) + ws(i, pointr)*neggi
1363 pointr = mod(pointr, m) + 1
1365 IF (nbd(i) <= 2 .AND. nbd(i) /= 0 &
1366 & .AND. neggi < zero)
THEN
1370 t(nbreak) = tl/(-neggi)
1371 IF (nbreak == 1 .OR. t(nbreak) < bkmin)
THEN
1375 ELSE IF (nbd(i) >= 2 .AND. neggi > zero)
THEN
1379 t(nbreak) = tu/neggi
1380 IF (nbreak == 1 .OR. t(nbreak) < bkmin)
THEN
1388 IF (abs(neggi) > zero) bnded = .false.
1397 IF (theta /= one)
THEN
1399 CALL dscal(col, theta, p(col + 1), 1)
1404 CALL dcopy(n, x, 1, xcp, 1)
1406 IF (nbreak == 0 .AND. nfree == n + 1)
THEN
1408 IF (iprint > 100)
WRITE (*, 1010) (xcp(i), i=1, n)
1423 CALL bmv(m, sy, wt, col, p, v, info)
1424 IF (info /= 0)
RETURN
1425 f2 = f2 - ddot(col2, v, 1, p, 1)
1431 WRITE (*, *)
'There are ', nbreak,
' breakpoints '
1440 IF (nleft == 0)
THEN
1441 IF (iprint >= 99)
THEN
1443 WRITE (*, *)
'GCP found in this segment'
1444 WRITE (*, 4010) nseg, f1, f2
1447 IF (dtm <= zero) dtm = zero
1453 CALL daxpy(n, tsum, d, 1, xcp, 1)
1456 DO WHILE (nleft > 0)
1467 ibp = iorder(ibkmin)
1472 IF (ibkmin /= nbreak)
THEN
1473 t(ibkmin) = t(nbreak)
1474 iorder(ibkmin) = iorder(nbreak)
1479 CALL hpsolb(nleft, t, iorder, iter - 2)
1486 IF (dt /= zero .AND. iprint >= 100)
THEN
1487 WRITE (*, 4011) nseg, f1, f2
1495 IF (iprint >= 99)
THEN
1497 WRITE (*, *)
'GCP found in this segment'
1498 WRITE (*, 4010) nseg, f1, f2
1501 IF (dtm <= zero) dtm = zero
1507 CALL daxpy(n, tsum, d, 1, xcp, 1)
1519 IF (dibp > zero)
THEN
1520 zibp = upper_bound(ibp) - x(ibp)
1521 xcp(ibp) = upper_bound(ibp)
1524 zibp = lower_bound(ibp) - x(ibp)
1525 xcp(ibp) = lower_bound(ibp)
1528 IF (iprint >= 100)
WRITE (*, *)
'Variable ', ibp,
' is fixed.'
1529 IF (nleft == 0 .AND. nbreak == n)
THEN
1544 f1 = f1 + dt*f2 + dibp2 - theta*dibp*zibp
1545 f2 = f2 - theta*dibp2
1549 CALL daxpy(col2, dt, p, 1, c, 1)
1555 wbp(j) = wy(ibp, pointr)
1556 wbp(col + j) = theta*ws(ibp, pointr)
1557 pointr = mod(pointr, m) + 1
1561 CALL bmv(m, sy, wt, col, wbp, v, info)
1562 IF (info /= 0)
RETURN
1563 wmc = ddot(col2, c, 1, v, 1)
1564 wmp = ddot(col2, p, 1, v, 1)
1565 wmw = ddot(col2, wbp, 1, v, 1)
1568 CALL daxpy(col2, -dibp, wbp, 1, p, 1)
1572 f2 = f2 + 2.0_dp*dibp*wmp - dibp2*wmw
1575 f2 = max(epsmch*f2_org, f2)
1588 IF (iprint >= 99)
THEN
1590 WRITE (*, *)
'GCP found in this segment'
1591 WRITE (*, 4010) nseg, f1, f2
1594 IF (dtm <= zero) dtm = zero
1600 CALL daxpy(n, tsum, d, 1, xcp, 1)
1608 IF (col > 0)
CALL daxpy(col2, dtm, p, 1, c, 1)
1609 IF (iprint > 100)
WRITE (*, 1010) (xcp(i), i=1, n)
1610 IF (iprint >= 99)
WRITE (*, 2010)
1612 1010
FORMAT(
'Cauchy X = ', /, (4x, 1p, 6(1x, d11.4)))
1613 2010
FORMAT(/,
'---------------- exit CAUCHY----------------------',/)
1614 3010
FORMAT(/,
'---------------- CAUCHY entered-------------------')
1615 4010
FORMAT(
'Piece ', i3,
' --f1, f2 at start point ', 1p, 2(1x, d11.4))
1616 4011
FORMAT(/,
'Piece ', i3,
' --f1, f2 at start point ', &
1618 5010
FORMAT(
'Distance to the next break point = ', 1p, d11.4)
1619 6010
FORMAT(
'Distance to the stationary point = ', 1p, d11.4)
1623 END SUBROUTINE cauchy
1653 SUBROUTINE cmprlb(n, m, x, g, ws, wy, sy, wt, z, r, wa, index, &
1654 theta, col, head, nfree, constrained, info)
1656 INTEGER,
INTENT(in) :: n, m
1657 REAL(kind=
dp),
INTENT(in) :: x(n), g(n), ws(n, m), wy(n, m), &
1658 sy(m, m), wt(m, m), z(n)
1659 REAL(kind=
dp),
INTENT(out) :: r(n), wa(4*m)
1660 INTEGER,
INTENT(in) :: index(n)
1661 REAL(kind=
dp),
INTENT(in) :: theta
1662 INTEGER,
INTENT(in) :: col, head, nfree
1663 LOGICAL,
INTENT(in) :: constrained
1666 INTEGER :: i, j, k, pointr
1667 REAL(kind=
dp) :: a1, a2
1669 IF (.NOT. constrained .AND. col > 0)
THEN
1676 r(i) = -theta*(z(k) - x(k)) - g(k)
1678 CALL bmv(m, sy, wt, col, wa(2*m + 1), wa(1), info)
1686 a2 = theta*wa(col + j)
1689 r(i) = r(i) + wy(k, pointr)*a1 + ws(k, pointr)*a2
1691 pointr = mod(pointr, m) + 1
1697 END SUBROUTINE cmprlb
1717 SUBROUTINE errclb(n, m, factr, lower_bound, upper_bound, nbd, task, info, k)
1719 INTEGER,
INTENT(in) :: n, m
1720 REAL(kind=
dp),
INTENT(in) :: factr, lower_bound(n), upper_bound(n)
1722 CHARACTER(LEN=60) :: task
1725 REAL(kind=
dp),
PARAMETER :: zero = 0.0_dp
1731 IF (n <= 0) task =
'ERROR: N <= 0'
1732 IF (m <= 0) task =
'ERROR: M <= 0'
1733 IF (factr < zero) task =
'ERROR: FACTR < 0'
1738 IF (nbd(i) < 0 .OR. nbd(i) > 3)
THEN
1740 task =
'ERROR: INVALID NBD'
1744 IF (nbd(i) == 2)
THEN
1745 IF (lower_bound(i) > upper_bound(i))
THEN
1747 task =
'ERROR: NO FEASIBLE SOLUTION'
1756 END SUBROUTINE errclb
1805 SUBROUTINE formk(n, nsub, ind, nenter, ileave, indx2, iupdat, &
1806 updatd, wn, wn1, m, ws, wy, sy, theta, col, &
1809 INTEGER,
INTENT(in) :: n, nsub, ind(n), nenter, ileave, &
1812 INTEGER,
INTENT(in) :: m
1813 REAL(kind=
dp) :: wn1(2*m, 2*m)
1814 REAL(kind=
dp),
INTENT(out) :: wn(2*m, 2*m)
1815 REAL(kind=
dp),
INTENT(in) :: ws(n, m), wy(n, m), sy(m, m), theta
1816 INTEGER,
INTENT(in) :: col, head
1817 INTEGER,
INTENT(out) :: info
1819 REAL(kind=
dp),
PARAMETER :: zero = 0.0_dp
1821 INTEGER :: col2, dbegin, dend, i, ipntr, is, is1, &
1822 iy, jpntr, js, js1, jy, k, k1, m2, &
1824 REAL(kind=
dp) :: ddot, temp1, temp2, temp3, temp4
1847 IF (iupdat > m)
THEN
1851 CALL dcopy(m - jy, wn1(jy + 1, jy + 1), 1, wn1(jy, jy), 1)
1852 CALL dcopy(m - jy, wn1(js + 1, js + 1), 1, wn1(js, js), 1)
1853 CALL dcopy(m - 1, wn1(m + 2, jy + 1), 1, wn1(m + 1, jy), 1)
1864 ipntr = head + col - 1
1865 IF (ipntr > m) ipntr = ipntr - m
1875 temp1 = temp1 + wy(k1, ipntr)*wy(k1, jpntr)
1880 temp2 = temp2 + ws(k1, ipntr)*ws(k1, jpntr)
1881 temp3 = temp3 + ws(k1, ipntr)*wy(k1, jpntr)
1886 jpntr = mod(jpntr, m) + 1
1891 jpntr = head + col - 1
1892 IF (jpntr > m) jpntr = jpntr - m
1900 temp3 = temp3 + ws(k1, ipntr)*wy(k1, jpntr)
1902 ipntr = mod(ipntr, m) + 1
1924 temp1 = temp1 + wy(k1, ipntr)*wy(k1, jpntr)
1925 temp2 = temp2 + ws(k1, ipntr)*ws(k1, jpntr)
1929 temp3 = temp3 + wy(k1, ipntr)*wy(k1, jpntr)
1930 temp4 = temp4 + ws(k1, ipntr)*ws(k1, jpntr)
1932 wn1(iy, jy) = wn1(iy, jy) + temp1 - temp3
1933 wn1(is, js) = wn1(is, js) - temp2 + temp4
1934 jpntr = mod(jpntr, m) + 1
1936 ipntr = mod(ipntr, m) + 1
1941 DO is = m + 1, m + upcl
1948 temp1 = temp1 + ws(k1, ipntr)*wy(k1, jpntr)
1952 temp3 = temp3 + ws(k1, ipntr)*wy(k1, jpntr)
1954 IF (is <= jy + m)
THEN
1955 wn1(is, jy) = wn1(is, jy) + temp1 - temp3
1957 wn1(is, jy) = wn1(is, jy) - temp1 + temp3
1959 jpntr = mod(jpntr, m) + 1
1961 ipntr = mod(ipntr, m) + 1
1974 wn(jy, iy) = wn1(iy, jy)/theta
1975 wn(js, is) = wn1(is1, js1)*theta
1978 wn(jy, is) = -wn1(is1, jy)
1981 wn(jy, is) = wn1(is1, jy)
1983 wn(iy, iy) = wn(iy, iy) + sy(iy, iy)
1991 CALL dpofa(wn, m2, col, info)
1998 DO js = col + 1, col2
1999 CALL dtrsl(wn, m2, col, wn(1, js), 11, info)
2005 DO is = col + 1, col2
2007 wn(is, js) = wn(is, js) + ddot(col, wn(1, is), 1, wn(1, js), 1)
2013 CALL dpofa(wn(col + 1, col + 1), m2, col, info)
2021 END SUBROUTINE formk
2042 SUBROUTINE formt(m, wt, sy, ss, col, theta, info)
2045 REAL(kind=
dp) :: wt(m, m), sy(m, m), ss(m, m)
2047 REAL(kind=
dp) :: theta
2050 REAL(kind=
dp),
PARAMETER :: zero = 0.0_dp
2052 INTEGER :: i, j, k, k1
2053 REAL(kind=
dp) :: ddum
2059 wt(1, j) = theta*ss(1, j)
2066 ddum = ddum + sy(i, k)*sy(j, k)/sy(k, k)
2068 wt(i, j) = ddum + theta*ss(i, j)
2075 CALL dpofa(wt, m, col, info)
2082 END SUBROUTINE formt
2115 SUBROUTINE freev(n, nfree, index, nenter, ileave, indx2, &
2116 iwhere, wrk, updatd, constrained, iprint, iter)
2119 INTEGER,
INTENT(inout) :: index(n)
2120 INTEGER :: nenter, ileave
2121 INTEGER,
INTENT(out) :: indx2(n)
2122 INTEGER :: iwhere(n)
2123 LOGICAL :: wrk, updatd, constrained
2124 INTEGER :: iprint, iter
2126 INTEGER :: i, iact, k
2130 IF (iter > 0 .AND. constrained)
THEN
2138 IF (iwhere(k) > 0)
THEN
2141 IF (iprint >= 100)
WRITE (*, *) &
2142 &
'Variable ', k,
' leaves the set of free variables'
2147 IF (iwhere(k) <= 0)
THEN
2150 IF (iprint >= 100)
WRITE (*, *) &
2151 &
'Variable ', k,
' enters the set of free variables'
2154 IF (iprint >= 99)
WRITE (*, *) &
2155 n + 1 - ileave,
' variables leave; ', nenter,
' variables enter'
2157 wrk = (ileave < n + 1) .OR. (nenter > 0) .OR. updatd
2164 IF (iwhere(i) <= 0)
THEN
2172 IF (iprint >= 99)
WRITE (*, *) &
2173 nfree,
' variables are free at GCP ', iter + 1
2177 END SUBROUTINE freev
2199 SUBROUTINE hpsolb(n, t, iorder, iheap)
2200 INTEGER,
INTENT(in) :: n
2201 REAL(kind=
dp),
INTENT(inout) :: t(n)
2202 INTEGER,
INTENT(inout) :: iorder(n)
2203 INTEGER,
INTENT(in) :: iheap
2205 INTEGER :: i, indxin, indxou, j, k
2206 REAL(kind=
dp) :: ddum, out
2214 IF (iheap == 0)
THEN
2226 IF (ddum < t(j))
THEN
2228 iorder(i) = iorder(j)
2252 DO WHILE (j <= n - 1)
2253 IF (t(j + 1) < t(j)) j = j + 1
2254 IF (t(j) < ddum)
THEN
2256 iorder(i) = iorder(j)
2274 END SUBROUTINE hpsolb
2317 SUBROUTINE lnsrlb(n, lower_bound, upper_bound, nbd, x, f, fold, gd, gdold, g, d, r, t, &
2318 z, stp, dnorm, dtd, xstep, step_max, iter, ifun, &
2319 iback, nfgv, info, task, boxed, constrained, csave, &
2322 INTEGER,
INTENT(in) :: n
2323 REAL(kind=
dp),
INTENT(in) :: lower_bound(n), upper_bound(n)
2325 REAL(kind=
dp) :: x(n), f, fold, gd, gdold, g(n), d(n), &
2326 r(n), t(n), z(n), stp, dnorm, dtd, &
2328 INTEGER :: iter, ifun, iback, nfgv, info
2329 CHARACTER(LEN=60) :: task
2330 LOGICAL :: boxed, constrained
2331 CHARACTER(LEN=60) :: csave
2333 REAL(kind=
dp) :: dsave(13)
2335 REAL(kind=
dp),
PARAMETER :: big = 1.0e10_dp, ftol = 1.0e-3_dp, &
2336 gtol = 0.9_dp, one = 1.0_dp, &
2337 xtol = 0.1_dp, zero = 0.0_dp
2340 REAL(kind=
dp) :: a1, a2, ddot
2342 IF (.NOT. (task(1:5) ==
'FG_LN'))
THEN
2344 dtd = ddot(n, d, 1, d, 1)
2350 IF (constrained)
THEN
2356 IF (nbd(i) /= 0)
THEN
2357 IF (a1 < zero .AND. nbd(i) <= 2)
THEN
2358 a2 = lower_bound(i) - x(i)
2359 IF (a2 >= zero)
THEN
2361 ELSE IF (a1*step_max < a2)
THEN
2364 ELSE IF (a1 > zero .AND. nbd(i) >= 2)
THEN
2365 a2 = upper_bound(i) - x(i)
2366 IF (a2 <= zero)
THEN
2368 ELSE IF (a1*step_max > a2)
THEN
2377 IF (iter == 0 .AND. .NOT. boxed)
THEN
2378 stp = min(one/dnorm, step_max)
2383 CALL dcopy(n, x, 1, t, 1)
2384 CALL dcopy(n, g, 1, r, 1)
2390 gd = ddot(n, g, 1, d, 1)
2393 IF (gd >= zero)
THEN
2396 WRITE (*, *)
' ascent direction in projection gd = ', gd
2402 CALL dcsrch(f, gd, stp, ftol, gtol, xtol, zero, step_max, csave, isave, dsave)
2405 IF (csave(1:4) /=
'CONV' .AND. csave(1:4) /=
'WARN')
THEN
2410 IF (stp == one)
THEN
2411 CALL dcopy(n, z, 1, x, 1)
2414 x(i) = stp*d(i) + t(i)
2423 END SUBROUTINE lnsrlb
2451 SUBROUTINE matupd(n, m, ws, wy, sy, ss, d, r, itail, &
2452 iupdat, col, head, theta, rr, dr, stp, dtd)
2455 REAL(kind=
dp) :: ws(n, m), wy(n, m), sy(m, m), ss(m, m), &
2457 INTEGER :: itail, iupdat, col, head
2458 REAL(kind=
dp) :: theta, rr, dr, stp, dtd
2460 REAL(kind=
dp),
PARAMETER :: one = 1.0_dp
2462 INTEGER :: j, pointr
2463 REAL(kind=
dp) :: ddot
2468 IF (iupdat <= m)
THEN
2470 itail = mod(head + iupdat - 2, m) + 1
2472 itail = mod(itail, m) + 1
2473 head = mod(head, m) + 1
2478 CALL dcopy(n, d, 1, ws(1, itail), 1)
2479 CALL dcopy(n, r, 1, wy(1, itail), 1)
2489 IF (iupdat > m)
THEN
2492 CALL dcopy(j, ss(2, j + 1), 1, ss(1, j), 1)
2493 CALL dcopy(col - j, sy(j + 1, j + 1), 1, sy(j, j), 1)
2500 sy(col, j) = ddot(n, d, 1, wy(1, pointr), 1)
2501 ss(j, col) = ddot(n, ws(1, pointr), 1, d, 1)
2502 pointr = mod(pointr, m) + 1
2504 IF (stp == one)
THEN
2507 ss(col, col) = stp*stp*dtd
2513 END SUBROUTINE matupd
2535 SUBROUTINE prn1lb(n, m, lower_bound, upper_bound, x, iprint, itfile, epsmch)
2537 INTEGER,
INTENT(in) :: n, m
2538 REAL(kind=
dp),
INTENT(in) :: lower_bound(n), upper_bound(n), x(n)
2539 INTEGER :: iprint, itfile
2540 REAL(kind=
dp) :: epsmch
2544 IF (iprint >= 0)
THEN
2545 WRITE (*, 7001) epsmch
2546 WRITE (*, *)
'N = ', n,
' M = ', m
2547 IF (iprint >= 1)
THEN
2548 WRITE (itfile, 2001) epsmch
2549 WRITE (itfile, *)
'N = ', n,
' M = ', m
2550 WRITE (itfile, 9001)
2551 IF (iprint > 100)
THEN
2552 WRITE (*, 1004)
'L =', (lower_bound(i), i=1, n)
2553 WRITE (*, 1004)
'X0 =', (x(i), i=1, n)
2554 WRITE (*, 1004)
'U =', (upper_bound(i), i=1, n)
2559 1004
FORMAT(/, a4, 1p, 6(1x, d11.4), /, (4x, 1p, 6(1x, d11.4)))
2560 2001
FORMAT(
'RUNNING THE L-BFGS-B CODE', /, /, &
2561 'it = iteration number', /, &
2562 'nf = number of function evaluations', /, &
2563 'nseg = number of segments explored during the Cauchy search', /, &
2564 'nact = number of active bounds at the generalized Cauchy point' &
2566 'sub = manner in which the subspace minimization terminated:' &
2567 , /,
' con = converged, bnd = a bound was reached', /, &
2568 'itls = number of iterations performed in the line search', /, &
2569 'stepl = step length used', /, &
2570 'tstep = norm of the displacement (total step)', /, &
2571 'projg = norm of the projected gradient', /, &
2572 'f = function value', /, /, &
2574 'Machine precision =', 1p, d10.3)
2575 7001
FORMAT(
'RUNNING THE L-BFGS-B CODE', /, /, &
2577 'Machine precision =', 1p, d10.3)
2578 9001
FORMAT(/, 3x,
'it', 3x,
'nf', 2x,
'nseg', 2x,
'nact', 2x,
'sub', 2x,
'itls', &
2579 2x,
'stepl', 4x,
'tstep', 5x,
'projg', 8x,
'f')
2583 END SUBROUTINE prn1lb
2610 SUBROUTINE prn2lb(n, x, f, g, iprint, itfile, iter, nfgv, nact, &
2611 g_inf_norm, nseg, word, iword, iback, stp, xstep)
2613 INTEGER,
INTENT(in) :: n
2614 REAL(kind=
dp),
INTENT(in) :: x(n), f, g(n)
2615 INTEGER,
INTENT(in) :: iprint, itfile, iter, nfgv, nact
2616 REAL(kind=
dp),
INTENT(in) :: g_inf_norm
2617 INTEGER,
INTENT(in) :: nseg
2618 CHARACTER(LEN=3) :: word
2619 INTEGER :: iword, iback
2620 REAL(kind=
dp) :: stp, xstep
2626 IF (iword == 0)
THEN
2629 ELSE IF (iword == 1)
THEN
2632 ELSE IF (iword == 5)
THEN
2638 IF (iprint >= 99)
THEN
2639 WRITE (*, *)
'LINE SEARCH', iback,
' times; norm of step = ', xstep
2640 WRITE (*, 2001) iter, f, g_inf_norm
2641 IF (iprint > 100)
THEN
2642 WRITE (*, 1004)
'X =', (x(i), i=1, n)
2643 WRITE (*, 1004)
'G =', (g(i), i=1, n)
2645 ELSE IF (iprint > 0)
THEN
2646 imod = mod(iter, iprint)
2647 IF (imod == 0)
WRITE (*, 2001) iter, f, g_inf_norm
2649 IF (iprint >= 1)
WRITE (itfile, 3001) &
2650 iter, nfgv, nseg, nact, word, iback, stp, xstep, g_inf_norm, f
2652 1004
FORMAT(/, a4, 1p, 6(1x, d11.4), /, (4x, 1p, 6(1x, d11.4)))
2654 (/,
'At iterate', i5, 4x,
'f= ', 1p, d12.5, 4x,
'|proj g|= ', 1p, d12.5)
2655 3001
FORMAT(2(1x, i4), 2(1x, i5), 2x, a3, 1x, i4, 1p, 2(2x, d7.1), 1p, 2(1x, d10.3))
2659 END SUBROUTINE prn2lb
2695 SUBROUTINE prn3lb(n, x, f, task, iprint, info, itfile, &
2696 iter, nfgv, nintol, nskip, nact, g_inf_norm, &
2697 time, nseg, word, iback, stp, xstep, k, &
2698 cachyt, sbtime, lnscht)
2700 INTEGER,
INTENT(in) :: n
2701 REAL(kind=
dp),
INTENT(in) :: x(n), f
2702 CHARACTER(LEN=60),
INTENT(in) :: task
2703 INTEGER,
INTENT(in) :: iprint, info, itfile, iter, nfgv, &
2705 REAL(kind=
dp),
INTENT(in) :: g_inf_norm, time
2706 INTEGER,
INTENT(in) :: nseg
2707 CHARACTER(LEN=3) :: word
2709 REAL(kind=
dp) :: stp, xstep
2711 REAL(kind=
dp) :: cachyt, sbtime, lnscht
2715 IF (iprint >= 0 .AND. .NOT. (task(1:5) ==
'ERROR'))
THEN
2718 WRITE (*, 3005) n, iter, nfgv, nintol, nskip, nact, g_inf_norm, f
2719 IF (iprint >= 100)
THEN
2720 WRITE (*, 1004)
'X =', (x(i), i=1, n)
2722 IF (iprint >= 1)
WRITE (*, *)
' F =', f
2724 IF (iprint >= 0)
THEN
2725 WRITE (*, 3009) task
2727 IF (info == -1)
WRITE (*, 9011)
2728 IF (info == -2)
WRITE (*, 9012)
2729 IF (info == -3)
WRITE (*, 9013)
2730 IF (info == -4)
WRITE (*, 9014)
2731 IF (info == -5)
WRITE (*, 9015)
2732 IF (info == -6)
WRITE (*, *)
' Input nbd(', k,
') is invalid.'
2734 WRITE (*, *)
' l(', k,
') > u(', k,
'). No feasible solution.'
2735 IF (info == -8)
WRITE (*, 9018)
2736 IF (info == -9)
WRITE (*, 9019)
2738 IF (iprint >= 1)
WRITE (*, 3007) cachyt, sbtime, lnscht
2739 WRITE (*, 3008) time
2740 IF (iprint >= 1)
THEN
2741 IF (info == -4 .OR. info == -9)
THEN
2742 WRITE (itfile, 3002) &
2743 iter, nfgv, nseg, nact, word, iback, stp, xstep
2745 WRITE (itfile, 3009) task
2747 IF (info == -1)
WRITE (itfile, 9011)
2748 IF (info == -2)
WRITE (itfile, 9012)
2749 IF (info == -3)
WRITE (itfile, 9013)
2750 IF (info == -4)
WRITE (itfile, 9014)
2751 IF (info == -5)
WRITE (itfile, 9015)
2752 IF (info == -8)
WRITE (itfile, 9018)
2753 IF (info == -9)
WRITE (itfile, 9019)
2755 WRITE (itfile, 3008) time
2759 1004
FORMAT(/, a4, 1p, 6(1x, d11.4), /, (4x, 1p, 6(1x, d11.4)))
2760 3002
FORMAT(2(1x, i4), 2(1x, i5), 2x, a3, 1x, i4, 1p, 2(2x, d7.1), 6x,
'-', 10x,
'-')
2763 'Tit = total number of iterations', /, &
2764 'Tnf = total number of function evaluations', /, &
2765 'Tnint = total number of segments explored during', &
2766 ' Cauchy searches', /, &
2767 'Skip = number of BFGS updates skipped', /, &
2768 'Nact = number of active bounds at final generalized', &
2769 ' Cauchy point', /, &
2770 'Projg = norm of the final projected gradient', /, &
2771 'F = final function value', /, /, &
2773 3004
FORMAT(/, 3x,
'N', 4x,
'Tit', 5x,
'Tnf', 2x,
'Tnint', 2x, &
2774 'Skip', 2x,
'Nact', 5x,
'Projg', 8x,
'F')
2775 3005
FORMAT(i5, 2(1x, i6), (1x, i6), (2x, i4), (1x, i5), 1p, 2(2x, d10.3))
2776 3007
FORMAT(/,
' Cauchy time', 1p, e10.3,
' seconds.', / &
2777 ' Subspace minimization time', 1p, e10.3,
' seconds.', / &
2778 ' Line search time', 1p, e10.3,
' seconds.')
2779 3008
FORMAT(/,
' Total User time', 1p, e10.3,
' seconds.',/)
2782 ' Matrix in 1st Cholesky factorization in formk is not Pos. Def.')
2784 ' Matrix in 2st Cholesky factorization in formk is not Pos. Def.')
2786 ' Matrix in the Cholesky factorization in formt is not Pos. Def.')
2788 ' Derivative >= 0, backtracking line search impossible.', /, &
2789 ' Previous x, f and g restored.', /, &
2790 ' Possible causes: 1 error in function or gradient evaluation;', /, &
2791 ' 2 rounding errors dominate computation.')
2793 ' Warning: more than 10 function and gradient', /, &
2794 ' evaluations in the last line search. Termination', /, &
2795 ' may possibly be caused by a bad search direction.')
2796 9018
FORMAT(/,
' The triangular system is singular.')
2798 ' Line search cannot locate an adequate point after 20 function', /, &
2799 ' and gradient evaluations. Previous x, f and g restored.', /, &
2800 ' Possible causes: 1 error in function or gradient evaluation;', /, &
2801 ' 2 rounding error dominate computation.')
2805 END SUBROUTINE prn3lb
2823 SUBROUTINE projgr(n, lower_bound, upper_bound, nbd, x, g, g_inf_norm)
2825 INTEGER,
INTENT(in) :: n
2826 REAL(kind=
dp),
INTENT(in) :: lower_bound(n), upper_bound(n)
2827 INTEGER,
INTENT(in) :: nbd(n)
2828 REAL(kind=
dp),
INTENT(in) :: x(n), g(n)
2829 REAL(kind=
dp) :: g_inf_norm
2831 REAL(kind=
dp),
PARAMETER :: zero = 0.0_dp
2839 IF (nbd(i) /= 0)
THEN
2841 IF (nbd(i) >= 2) gi = max((x(i) - upper_bound(i)), gi)
2843 IF (nbd(i) <= 2) gi = min((x(i) - lower_bound(i)), gi)
2846 g_inf_norm = max(g_inf_norm, abs(gi))
2851 END SUBROUTINE projgr
2958 SUBROUTINE subsm(n, m, nsub, ind, lower_bound, upper_bound, nbd, x, d, xp, ws, wy, &
2960 col, head, iword, wv, wn, iprint, info)
2961 INTEGER,
INTENT(in) :: n, m, nsub, ind(nsub)
2962 REAL(kind=
dp),
INTENT(in) :: lower_bound(n), upper_bound(n)
2963 INTEGER,
INTENT(in) :: nbd(n)
2964 REAL(kind=
dp),
INTENT(inout) :: x(n), d(n)
2965 REAL(kind=
dp) :: xp(n)
2966 REAL(kind=
dp),
INTENT(in) :: ws(n, m), wy(n, m), theta, xx(n), gg(n)
2967 INTEGER,
INTENT(in) :: col, head
2968 INTEGER,
INTENT(out) :: iword
2969 REAL(kind=
dp) :: wv(2*m)
2970 REAL(kind=
dp),
INTENT(in) :: wn(2*m, 2*m)
2972 INTEGER,
INTENT(out) :: info
2974 REAL(kind=
dp),
PARAMETER :: one = 1.0_dp, zero = 0.0_dp
2976 INTEGER :: col2, i, ibd, j, js, jy, k, m2, pointr
2977 REAL(kind=
dp) :: alpha, dd_p, dk, temp1, temp2, xk
2990 IF (nsub <= 0)
RETURN
2991 IF (iprint >= 99)
WRITE (*, 1001)
3001 temp1 = temp1 + wy(k, pointr)*d(j)
3002 temp2 = temp2 + ws(k, pointr)*d(j)
3005 wv(col + i) = theta*temp2
3006 pointr = mod(pointr, m) + 1
3013 CALL dtrsl(wn, m2, col2, wv, 11, info)
3014 IF (info /= 0)
RETURN
3018 CALL dtrsl(wn, m2, col2, wv, 01, info)
3019 IF (info /= 0)
RETURN
3028 d(i) = d(i) + wy(k, pointr)*wv(jy)/theta &
3029 & + ws(k, pointr)*wv(js)
3031 pointr = mod(pointr, m) + 1
3034 CALL dscal(nsub, one/theta, d, 1)
3041 CALL dcopy(n, x, 1, xp, 1)
3047 IF (nbd(k) /= 0)
THEN
3050 IF (nbd(k) .EQ. 1)
THEN
3051 x(k) = max(lower_bound(k), xk + dk)
3052 IF (x(k) .EQ. lower_bound(k)) iword = 1
3056 IF (nbd(k) .EQ. 2)
THEN
3057 xk = max(lower_bound(k), xk + dk)
3058 x(k) = min(upper_bound(k), xk)
3059 IF (x(k) .EQ. lower_bound(k) .OR. x(k) .EQ. upper_bound(k)) iword = 1
3063 IF (nbd(k) .EQ. 3)
THEN
3064 x(k) = min(upper_bound(k), xk + dk)
3065 IF (x(k) .EQ. upper_bound(k)) iword = 1
3076 IF (.NOT. (iword .EQ. 0))
THEN
3082 dd_p = dd_p + (x(i) - xx(i))*gg(i)
3084 IF (dd_p .GT. zero)
THEN
3085 CALL dcopy(n, xp, 1, x, 1)
3086 IF (iprint > 0)
WRITE (*, *)
' Positive dir derivative in projection '
3087 IF (iprint > 0)
WRITE (*, *)
' Using the backtracking step '
3094 IF (nbd(k) /= 0)
THEN
3095 IF (dk < zero .AND. nbd(k) <= 2)
THEN
3096 temp2 = lower_bound(k) - x(k)
3097 IF (temp2 >= zero)
THEN
3099 ELSE IF (dk*alpha < temp2)
THEN
3102 ELSE IF (dk > zero .AND. nbd(k) >= 2)
THEN
3103 temp2 = upper_bound(k) - x(k)
3104 IF (temp2 <= zero)
THEN
3106 ELSE IF (dk*alpha > temp2)
THEN
3110 IF (temp1 < alpha)
THEN
3117 IF (alpha < one)
THEN
3121 x(k) = upper_bound(k)
3123 ELSE IF (dk < zero)
THEN
3124 x(k) = lower_bound(k)
3130 x(k) = x(k) + alpha*d(i)
3135 IF (iprint >= 99)
WRITE (*, 1004)
3137 1001
FORMAT(/,
'----------------SUBSM entered-----------------',/)
3138 1004
FORMAT(/,
'----------------exit SUBSM --------------------',/)
3142 END SUBROUTINE subsm
3230 SUBROUTINE dcsrch(f, g, stp, ftol, gtol, xtol, stpmin, stpmax, &
3232 REAL(kind=
dp) :: f, g
3233 REAL(kind=
dp),
INTENT(inout) :: stp
3234 REAL(kind=
dp) :: ftol, gtol, xtol, stpmin, stpmax
3235 CHARACTER(LEN=*) :: task
3237 REAL(kind=
dp) :: dsave(13)
3239 REAL(kind=
dp),
PARAMETER :: p5 = 0.5_dp, p66 = 0.66_dp, &
3240 xtrapl = 1.1_dp, xtrapu = 4.0_dp, &
3245 REAL(kind=
dp) :: finit, fm, ftest, fx, fxm, fy, fym, &
3246 ginit, gm, gtest, gx, gxm, gy, gym, &
3247 stmax, stmin, stx, sty, width, width1
3264 IF (task(1:5) ==
'START')
THEN
3268 IF (stp < stpmin) task =
'ERROR: STP < STPMIN'
3269 IF (stp > stpmax) task =
'ERROR: STP > STPMAX'
3270 IF (g >= zero) task =
'ERROR: INITIAL G >= ZERO'
3271 IF (ftol < zero) task =
'ERROR: FTOL < ZERO'
3272 IF (gtol < zero) task =
'ERROR: GTOL < ZERO'
3273 IF (xtol < zero) task =
'ERROR: XTOL < ZERO'
3274 IF (stpmin < zero) task =
'ERROR: STPMIN < ZERO'
3275 IF (stpmax < stpmin) task =
'ERROR: STPMAX < STPMIN'
3279 IF (task(1:5) ==
'ERROR')
RETURN
3288 width = stpmax - stpmin
3305 stmax = stp + xtrapu*stp
3312 IF (isave(1) == 1)
THEN
3335 ftest = finit + stp*gtest
3336 IF (stage == 1 .AND. f <= ftest .AND. g >= zero) &
3341 IF (brackt .AND. (stp <= stmin .OR. stp >= stmax)) &
3342 task =
'WARNING: ROUNDING ERRORS PREVENT PROGRESS'
3343 IF (brackt .AND. stmax - stmin <= xtol*stmax) &
3344 task =
'WARNING: XTOL TEST SATISFIED'
3345 IF (stp == stpmax .AND. f <= ftest .AND. g <= gtest) &
3346 task =
'WARNING: STP = STPMAX'
3347 IF (stp == stpmin .AND. (f > ftest .OR. g >= gtest)) &
3348 task =
'WARNING: STP = STPMIN'
3352 IF (f <= ftest .AND. abs(g) <= gtol*(-ginit)) &
3353 task =
'CONVERGENCE'
3357 IF (.NOT. (task(1:4) ==
'WARN' .OR. task(1:4) ==
'CONV'))
THEN
3363 IF (stage == 1 .AND. f <= fx .AND. f > ftest)
THEN
3368 fxm = fx - stx*gtest
3369 fym = fy - sty*gtest
3376 CALL dcstep(stx, fxm, gxm, sty, fym, gym, stp, fm, gm, &
3377 brackt, stmin, stmax)
3381 fx = fxm + stx*gtest
3382 fy = fym + sty*gtest
3390 CALL dcstep(stx, fx, gx, sty, fy, gy, stp, f, g, &
3391 brackt, stmin, stmax)
3398 IF (abs(sty - stx) >= p66*width1) stp = stx + p5*(sty - stx)
3400 width = abs(sty - stx)
3406 stmin = min(stx, sty)
3407 stmax = max(stx, sty)
3409 stmin = stp + xtrapl*(stp - stx)
3410 stmax = stp + xtrapu*(stp - stx)
3415 stp = max(stp, stpmin)
3416 stp = min(stp, stpmax)
3421 IF (brackt .AND. (stp <= stmin .OR. stp >= stmax) &
3422 .OR. (brackt .AND. stmax - stmin <= xtol*stmax)) stp = stx
3454 END SUBROUTINE dcsrch
3499 SUBROUTINE dcstep(stx, fx, dx, sty, fy, dy, stp, fp, dp_loc, brackt, &
3501 REAL(kind=
dp),
INTENT(inout) :: stx, fx, dx, sty, fy, dy, stp
3502 REAL(kind=
dp),
INTENT(in) :: fp, dp_loc
3503 LOGICAL,
INTENT(inout) :: brackt
3504 REAL(kind=
dp),
INTENT(in) :: stpmin, stpmax
3506 REAL(kind=
dp),
PARAMETER :: p66 = 0.66_dp, three = 3.0_dp, &
3507 two = 2.0_dp, zero = 0.0_dp
3509 REAL(kind=
dp) ::
gamma, p, q, r, s, sgnd, stpc, stpf, &
3523 sgnd = dp_loc*sign(1.0_dp, dx)
3531 theta = three*(fx - fp)/(stp - stx) + dx + dp_loc
3532 s = max(abs(theta), abs(dx), abs(dp_loc))
3533 gamma = s*sqrt((theta/s)**2 - (dx/s)*(dp_loc/s))
3535 p = (
gamma - dx) + theta
3538 stpc = stx + r*(stp - stx)
3539 stpq = stx + ((dx/((fx - fp)/(stp - stx) + dx))/two)* &
3541 IF (abs(stpc - stx) < abs(stpq - stx))
THEN
3544 stpf = stpc + (stpq - stpc)/two
3553 ELSE IF (sgnd < zero)
THEN
3554 theta = three*(fx - fp)/(stp - stx) + dx + dp_loc
3555 s = max(abs(theta), abs(dx), abs(dp_loc))
3556 gamma = s*sqrt((theta/s)**2 - (dx/s)*(dp_loc/s))
3558 p = (
gamma - dp_loc) + theta
3561 stpc = stp + r*(stx - stp)
3562 stpq = stp + (dp_loc/(dp_loc - dx))*(stx - stp)
3563 IF (abs(stpc - stp) > abs(stpq - stp))
THEN
3573 ELSE IF (abs(dp_loc) < abs(dx))
THEN
3580 theta = three*(fx - fp)/(stp - stx) + dx + dp_loc
3581 s = max(abs(theta), abs(dx), abs(dp_loc))
3586 gamma = s*sqrt(max(zero, (theta/s)**2 - (dx/s)*(dp_loc/s)))
3588 p = (
gamma - dp_loc) + theta
3591 IF (r < zero .AND.
gamma /= zero)
THEN
3592 stpc = stp + r*(stx - stp)
3593 ELSE IF (stp > stx)
THEN
3598 stpq = stp + (dp_loc/(dp_loc - dx))*(stx - stp)
3606 IF (abs(stpc - stp) < abs(stpq - stp))
THEN
3612 stpf = min(stp + p66*(sty - stp), stpf)
3614 stpf = max(stp + p66*(sty - stp), stpf)
3622 IF (abs(stpc - stp) > abs(stpq - stp))
THEN
3627 stpf = min(stpmax, stpf)
3628 stpf = max(stpmin, stpf)
3638 theta = three*(fp - fy)/(sty - stp) + dy + dp_loc
3639 s = max(abs(theta), abs(dy), abs(dp_loc))
3640 gamma = s*sqrt((theta/s)**2 - (dy/s)*(dp_loc/s))
3642 p = (
gamma - dp_loc) + theta
3645 stpc = stp + r*(sty - stp)
3647 ELSE IF (stp > stx)
THEN
3661 IF (sgnd < zero)
THEN
3676 END SUBROUTINE dcstep
3700 SUBROUTINE dpofa(a, lda, n, info)
3701 INTEGER,
INTENT(in) :: lda
3702 REAL(kind=
dp) :: a(lda, *)
3703 INTEGER,
INTENT(in) :: n
3706 INTEGER :: j, jm1, k
3707 REAL(kind=
dp) :: ddot, s, t
3721 IF (.NOT. (jm1 < 1))
THEN
3723 t = a(k, j) - ddot(k - 1, a(1, k), 1, a(1, j), 1)
3731 IF (s <= 0.0_dp)
EXIT
3736 END SUBROUTINE dpofa
3768 SUBROUTINE dtrsl(t, ldt, n, b, job, info)
3769 INTEGER,
INTENT(in) :: ldt
3770 REAL(kind=
dp),
INTENT(in) :: t(ldt, *)
3771 INTEGER,
INTENT(in) :: n
3772 REAL(kind=
dp),
INTENT(inout) :: b(*)
3773 INTEGER,
INTENT(in) :: job
3774 INTEGER,
INTENT(out) :: info
3776 INTEGER :: case, j, jj
3777 REAL(kind=
dp) :: ddot, temp
3789 IF (t(info, info) == 0.0_dp)
RETURN
3796 IF (mod(job, 10) /= 0)
CASE = 2
3797 IF (mod(job, 100)/10 /= 0)
CASE =
CASE + 2
3808 CALL daxpy(n - j + 1, temp, t(j, j - 1), 1, b(j), 1)
3821 CALL daxpy(j, temp, t(1, j + 1), 1, b(1), 1)
3833 b(j) = b(j) - ddot(jj - 1, t(j + 1, j), 1, b(j + 1), 1)
3842 IF (.NOT. (n < 2))
THEN
3844 b(j) = b(j) - ddot(j - 1, t(1, j), 1, b(1), 1)
3849 cpabort(
"unexpected case")
3853 END SUBROUTINE dtrsl
3863 SUBROUTINE timer(ttime)
3864 REAL(kind=
dp) :: ttime
3885 END SUBROUTINE timer
3972 SUBROUTINE save_local(lsave,isave,dsave,x_projected,constrained,boxed,updatd,nintol,itfile,iback,nskip,head,col,itail,&
3973 iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, cpu1, &
3974 cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
3975 LOGICAL,
INTENT(out) :: lsave(4)
3976 INTEGER,
INTENT(out) :: isave(23)
3977 REAL(kind=
dp),
INTENT(out) :: dsave(29)
3978 LOGICAL,
INTENT(in) :: x_projected, constrained, boxed, updatd
3979 INTEGER,
INTENT(in) :: nintol, itfile, iback, nskip, head, col, &
3980 itail, iter, iupdat, nseg, nfgv, info, &
3981 ifun, iword, nfree, nact, ileave, &
3983 REAL(kind=
dp),
INTENT(in) :: theta, fold, tol, dnorm, epsmch, cpu1, &
3984 cachyt, sbtime, lnscht, time1, gd, &
3985 step_max, g_inf_norm, stp, gdold, dtd
3987 lsave(1) = x_projected
3988 lsave(2) = constrained
4022 dsave(12) = step_max
4023 dsave(13) = g_inf_norm
4028 END SUBROUTINE save_local
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public byrd1995
Utility routines to open and close files. Tracking of preconnections.
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
LBFGS-B routine (version 3.0, April 25, 2011)
subroutine, public setulb(n, m, x, lower_bound, upper_bound, nbd, f, g, factr, pgtol, wa, iwa, task, iprint, csave, lsave, isave, dsave, trust_radius, spgr)
This subroutine partitions the working arrays wa and iwa, and then uses the limited memory BFGS metho...
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Defines the basic variable types.
integer, parameter, public dp
Machine interface based on Fortran 2003 and POSIX.
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Space Group Symmetry Type Module (version 1.0, Ferbruary 12, 2021)
Space Group Symmetry Module (version 1.0, January 16, 2020)
subroutine, public spgr_apply_rotations_coord(spgr, coord)
routine applies the rotation matrices to the coordinates.
subroutine, public spgr_apply_rotations_force(spgr, force)
routine applies the rotation matrices to the forces.