97 MODULE cg_minimization
103 integer :: n, n5, n6, nf, ng
105 integer :: nrestart, nexpand, nsecant
107 real :: delta, sigma, eps
108 real :: gamma, rho, tol
109 real :: eta, fpert, f0
111 real :: wolfe_hi, wolfe_lo, awolfe_hi
112 real :: quadcutoff, stopfac, awolfefac
114 real :: psi0, psi1, psi2
115 logical :: pertrule, quadok
116 logical :: quadstep, printlevel
117 logical :: printfinal, stoprule
123 real,
dimension(:),
allocatable :: x_best
135 SUBROUTINE cg_descent (grad_tol, x, dim, cg_value, cg_grad, &
136 STATUS, gnorm, f, iter, nfunc, ngrad, d, g, xtemp, gtemp)
141 INTEGER,
INTENT(IN) :: dim
142 REAL,
INTENT(IN) :: grad_tol
143 REAL,
INTENT(IN OUT) :: x(dim)
144 INTEGER,
INTENT(OUT) :: status
145 REAL,
INTENT(IN OUT) :: gnorm
146 REAL,
INTENT(IN OUT) :: f
147 INTEGER,
INTENT(OUT) :: iter
148 INTEGER,
INTENT(OUT) :: nfunc
149 INTEGER,
INTENT(OUT) :: ngrad
150 REAL,
INTENT(IN OUT) :: d(dim)
151 REAL,
INTENT(IN OUT) :: g(dim)
152 REAL,
INTENT(IN OUT) :: xtemp(dim)
153 REAL,
INTENT(IN OUT) :: gtemp(dim)
157 SUBROUTINE cg_value (f, x, n)
159 INTEGER,
INTENT(IN) :: n
160 REAL,
INTENT(IN) :: x(:)
161 REAL,
INTENT(IN OUT) :: f
162 END SUBROUTINE cg_value
164 SUBROUTINE cg_grad (g, x, n)
166 INTEGER,
INTENT(IN) :: n
167 REAL,
INTENT(IN) :: x(:)
168 REAL,
INTENT(IN OUT) :: g(:)
169 END SUBROUTINE cg_grad
173 integer :: i, j, i1, i2, i3, i4
174 real :: delta2, eta_sq, qk
175 real :: ftemp, xnorm, gnorm2
176 real :: dnorm2, denom
177 real :: t, t1, t2, t3, t4
179 real :: alpha, talpha
180 real :: yk, yk2, ykgk, dkyk
184 allocate( x_best( dim ) )
189 CALL cg_init (grad_tol, dim)
193 delta2 = 2.0 * delta - 1.0
201 CALL cg_value (f, x, n)
202 if ( f < f_best ) x_best(1:n) = x(1:n)
203 if ( f < f_best ) f_best = f
205 CALL cg_grad (g, x, n)
212 xnorm = max(xnorm, abs(x(i)))
215 gnorm = max(gnorm, abs(t))
216 gnorm2 = gnorm2 + t*t
219 xnorm = max(xnorm, abs(x(i)))
221 gnorm = max(gnorm, abs(t))
226 gnorm = max(gnorm, abs(t1))
227 xnorm = max(xnorm, abs(x(j)))
231 gnorm = max(gnorm, abs(t2))
232 xnorm = max(xnorm, abs(x(j)))
236 gnorm = max(gnorm, abs(t3))
237 xnorm = max(xnorm, abs(x(j)))
241 gnorm = max(gnorm, abs(t4))
242 xnorm = max(xnorm, abs(x(j)))
243 gnorm2 = gnorm2 + t*t + t1*t1 + t2*t2 + t3*t3 + t4*t4
247 tol = max(gnorm*stopfac, tol)
250 IF ( printlevel )
THEN 251 WRITE (*, 10) iter, f, gnorm, awolfe
252 10
FORMAT (
'iter: ', i5,
' f= ', e14.6, &
253 ' gnorm= ', e14.6,
' AWolfe= ', l2)
256 IF ( cg_tol(f, gnorm) )
GO TO 100
259 IF ( .NOT.step )
THEN 260 alpha = psi0*xnorm/gnorm
261 IF ( xnorm == zero )
THEN 262 IF ( f /= zero )
THEN 263 alpha = psi0*abs(f)/gnorm2
282 IF ( f /= zero )
THEN 287 IF ( t > quadcutoff )
THEN 289 CALL cg_step(xtemp, x, d, talpha)
290 CALL cg_value(ftemp, xtemp, n)
291 if ( ftemp < f_best ) x_best(1:n) = xtemp(1:n)
292 if ( ftemp < f_best ) f_best = ftemp
294 IF ( ftemp < f )
THEN 295 denom = 2.0d0*(((ftemp-f)/talpha)-dphi0)
296 IF ( denom > zero )
THEN 298 alpha = -dphi0*talpha/denom
305 IF ( printlevel )
THEN 306 WRITE (*, 20) quadok, alpha, f0, dphi0
307 20
FORMAT (
'QuadOK:', l2,
' initial a:', &
308 e14.6,
' f0:', e14.6,
' dphi', e14.6)
314 ck = ck + (abs(f) - ck)/qk
322 wolfe_hi = delta*dphi0
323 wolfe_lo = sigma*dphi0
324 awolfe_hi = delta2*dphi0
326 CALL cg_line (alpha, f, dphi, dphi0, x, xtemp, d, gtemp, &
329 CALL cg_linew (alpha, f, dphi, dphi0, x, xtemp, d, gtemp, &
333 IF ( info > 0 )
GO TO 100
338 IF ( -alpha*dphi0 <= feps*abs(f) )
THEN 345 IF ( mod(iter, nrestart) /= 0 )
THEN 357 gnorm = max(gnorm, abs(t))
358 dnorm2 = dnorm2 + d(i)**2
378 yk2 = yk2 + (t1-g(i1))**2 + (t2-g(i2))**2 &
379 + (t3-g(i3))**2 + (t4-g(i4))**2
380 ykgk = ykgk + (t1-g(i1))*t1 + (t2-g(i2))*t2 &
381 + (t3-g(i3))*t3 + (t4-g(i4))*t4
383 gnorm = max(gnorm, abs(t))
385 gnorm = max(gnorm, abs(t1))
387 gnorm = max(gnorm, abs(t2))
389 gnorm = max(gnorm, abs(t3))
391 gnorm = max(gnorm, abs(t4))
392 dnorm2 = dnorm2 + d(i)**2 + d(i1)**2 + d(i2)**2 &
393 + d(i3)**2 + d(i4)**2
395 IF ( cg_tol(f, gnorm) )
GO TO 100
397 beta = (ykgk - 2.0*dphi*yk2/dkyk)/dkyk
405 beta = max(beta, -1.0/sqrt(min(eta_sq, gnorm2)*dnorm2))
412 d(i) = -t + beta*d(i)
413 gnorm2 = gnorm2 + t*t
416 d(i) = -g(i) + beta*d(i)
418 d(i1) = -g(i1) + beta*d(i1)
420 d(i2) = -g(i2) + beta*d(i2)
422 d(i3) = -g(i3) + beta*d(i3)
424 d(i4) = -g(i4) + beta*d(i4)
425 gnorm2 = gnorm2 + g(i)**2 + g(i1)**2 + g(i2)**2 &
426 + g(i3)**2 + g(i4)**2
428 dphi0 = -gnorm2 + beta*dphi
434 IF ( printlevel )
THEN 435 WRITE (*, *)
"RESTART CG" 444 gnorm = max(gnorm, abs(t))
445 gnorm2 = gnorm2 + t*t
452 gnorm = max(gnorm, abs(t))
458 gnorm = max(gnorm, abs(t1))
464 gnorm = max(gnorm, abs(t2))
470 gnorm = max(gnorm, abs(t3))
476 gnorm = max(gnorm, abs(t4))
477 gnorm2 = gnorm2 + t*t + t1*t1 + t2*t2 + t3*t3 + t4*t4
479 IF ( cg_tol(f, gnorm) )
GO TO 100
482 IF ( .NOT.awolfe )
THEN 483 IF ( abs(f-f0) < awolfefac*ck )
THEN 488 IF ( printlevel )
THEN 489 WRITE (*, 10) iter, f, gnorm, awolfe
493 IF ( f > f0 + 1.e-10*ck )
THEN 496 50
FORMAT (
' new value:', e30.16,
'old value:', e30.16)
501 IF ( dphi0 > zero )
THEN 515 gnorm = max(gnorm, abs(g(i)))
518 IF ( printfinal )
THEN 519 WRITE (6, *)
'Termination status:', status
520 IF ( status == 0 )
THEN 522 ELSE IF ( status == 1 )
THEN 524 ELSE IF ( status == 2 )
THEN 527 WRITE (6, 400) grad_tol
528 ELSE IF ( status == 3 )
THEN 533 ELSE IF ( status == 4 )
THEN 536 WRITE (6, 400) grad_tol
537 ELSE IF ( status == 5 )
THEN 539 ELSE IF ( status == 6 )
THEN 542 WRITE (6, 400) grad_tol
545 ELSE IF ( status == 7 )
THEN 548 WRITE (6, 400) grad_tol
549 ELSE IF ( status == 8 )
THEN 552 WRITE (6, 400) grad_tol
557 WRITE (6, *)
'function value:', f
558 WRITE (6, *)
'cg iterations:', iter
559 WRITE (6, *)
'function evaluations:', nfunc
560 WRITE (6, *)
'gradient evaluations:', ngrad
563 if ( f_best < f )
then 565 x(1:dim) = x_best(1:dim)
572 200
FORMAT (
' Convergence tolerance for gradient satisfied')
573 210
FORMAT (
' Terminating since change in function value <= feps*|f|')
574 220
FORMAT (
' Total number of iterations exceed max allow:', i10)
575 230
FORMAT (
' Slope always negative in line search')
576 240
FORMAT (
' Line search fails, too many secant steps')
577 250
FORMAT (
' Search direction not a descent direction')
578 260
FORMAT (
' Line search fails')
579 270
FORMAT (
' Debugger is on, function value does not improve')
580 300
FORMAT (
' Possible causes of this error message:')
581 400
FORMAT (
' - your tolerance (grad_tol = ', d11.4, &
582 ') may be too strict')
583 410
FORMAT (
' - your gradient routine has an error')
584 420
FORMAT (
' - parameter epsilon in cg_descent.parm is too small')
585 430
FORMAT (
' - your cost function has an error')
586 500
FORMAT (
' absolute largest component of gradient: ', d11.4)
588 END SUBROUTINE cg_descent
671 SUBROUTINE cg_init (grad_tol, dim)
676 real,
INTENT(IN) :: grad_tol
677 integer,
INTENT(IN) :: dim
679 real :: restart_fac, maxit_fac
684 OPEN (10, file=
'cg_descent_f.parm')
694 READ (10, *) quadcutoff
696 READ (10, *) awolfefac
697 READ (10, *) restart_fac
698 READ (10, *) maxit_fac
703 READ (10, *) pertrule
704 READ (10, *) quadstep
705 READ (10, *) printlevel
706 READ (10, *) printfinal
707 READ (10, *) stoprule
711 nrestart = n*restart_fac
721 END SUBROUTINE cg_init
734 LOGICAL FUNCTION cg_wolfe (alpha, f, dphi)
739 real,
INTENT(IN) :: alpha
740 real,
INTENT(IN) :: f
741 real,
INTENT(IN) :: dphi
744 IF ( dphi >= wolfe_lo )
THEN 748 IF ( f-f0 <= alpha*wolfe_hi )
THEN 750 IF ( printlevel )
THEN 751 WRITE (*, 10) f, f0, alpha*wolfe_hi, dphi
752 10
FORMAT (
' wolfe f:', e14.6,
' f0: ', &
753 e14.6, e14.6,
' dphi:', e14.6)
759 ELSE IF ( awolfe )
THEN 760 IF ( (f <= fpert).AND.(dphi <= awolfe_hi) )
THEN 762 IF ( printlevel )
THEN 763 WRITE (*, 20) f, fpert, dphi, awolfe_hi
764 20
FORMAT (
'f:', e14.6,
' fpert:', e14.6, &
765 ' dphi: ', e14.6,
' fappx:', e14.6)
773 END FUNCTION cg_wolfe
783 LOGICAL FUNCTION cg_tol (f, gnorm)
788 real,
INTENT(IN OUT) :: f
789 real,
INTENT(IN) :: gnorm
793 IF ( gnorm <= tol )
THEN 798 IF ( gnorm <= tol*(1.0 + abs(f)) )
THEN 815 real FUNCTION cg_dot (x, y)
820 real,
INTENT(IN) :: x(:)
821 real,
INTENT(IN) :: y(:)
832 t = t + x(i)*y(i) + x(i+1)*y(i+1) + x(i+2)*y(i+2) &
833 + x(i+3)*y(i+3) + x(i+4)*y(i+4)
849 SUBROUTINE cg_step (xtemp, x, d, alpha)
854 real,
INTENT(OUT) :: xtemp(:)
855 real,
INTENT(IN) :: x(:)
856 real,
INTENT(IN) :: d(:)
857 real,
INTENT(IN) :: alpha
863 xtemp(i) = x(i) + alpha*d(i)
866 xtemp(i) = x(i) + alpha*d(i)
868 xtemp(j) = x(j) + alpha*d(j)
870 xtemp(j) = x(j) + alpha*d(j)
872 xtemp(j) = x(j) + alpha*d(j)
874 xtemp(j) = x(j) + alpha*d(j)
876 END SUBROUTINE cg_step
893 SUBROUTINE cg_line (alpha, phi, dphi, dphi0, x, xtemp, d, gtemp, &
900 real,
INTENT(IN OUT) :: alpha
901 real,
INTENT(IN OUT) :: phi
902 real,
INTENT(OUT) :: dphi
903 real,
INTENT(IN) :: dphi0
904 real,
INTENT(IN OUT) :: x(:)
905 real,
INTENT(IN OUT) :: xtemp(:)
906 real,
INTENT(IN OUT) :: d(:)
907 real,
INTENT(IN OUT) :: gtemp(:)
911 SUBROUTINE cg_value (f, x, n)
913 INTEGER,
INTENT(IN) :: n
914 REAL,
INTENT(IN) :: x(:)
915 REAL,
INTENT(IN OUT) :: f
916 END SUBROUTINE cg_value
918 SUBROUTINE cg_grad (g, x, n)
920 INTEGER,
INTENT(IN) :: n
921 REAL,
INTENT(IN) :: x(:)
922 REAL,
INTENT(IN OUT) :: g(:)
923 END SUBROUTINE cg_grad
927 real :: a, dphia, b, dphib, c
928 real :: a0, da0, b0, db0
931 integer :: ngrow, nshrink, iter, flag
934 CALL cg_step (xtemp, x, d, alpha)
935 CALL cg_grad (gtemp, xtemp, n)
937 dphi = cg_dot(gtemp, d)
946 DO WHILE ( dphi < zero )
947 CALL cg_value (phi, xtemp, n)
953 IF ( ngrow == 0 ) fquad = min(phi, f0)
954 IF ( phi <= fquad )
THEN 955 IF ( printlevel )
THEN 956 WRITE (*, 10) alpha, phi, fquad
957 10
FORMAT (
'alpha:', e14.6,
' phi:', e14.6, &
960 IF ( cg_wolfe(alpha, phi, dphi) )
RETURN 963 IF ( phi <= fpert )
THEN 973 nshrink = nshrink + 1
974 IF ( nshrink > nexpand )
THEN 978 CALL cg_step(xtemp, x, d, alpha)
979 CALL cg_grad(gtemp, xtemp, n)
981 dphi = cg_dot(gtemp, d)
982 IF ( dphi >= zero )
GO TO 100
983 CALL cg_value(phi, xtemp, n)
985 IF ( printlevel )
THEN 986 WRITE (6, 20) a, b, alpha, phi, dphi
987 20
FORMAT (
'contract, a:', e14.6, &
988 ' b:', e14.6,
' alpha:', e14.6,
' phi:', e14.6,
' dphi:', e14.6)
990 IF ( quadok .AND. (phi <= fquad) )
THEN 991 IF ( cg_wolfe(alpha, phi, dphi) )
RETURN 993 IF ( phi <= fpert )
THEN 1005 IF ( ngrow > nexpand )
THEN 1010 CALL cg_step(xtemp, x, d, alpha)
1011 CALL cg_grad(gtemp, xtemp, n)
1013 dphi = cg_dot(gtemp, d)
1014 IF ( printlevel )
THEN 1015 WRITE (*, 30) a, alpha, phi, dphi
1016 30
FORMAT (
'expand, a:', e14.6,
' alpha:', e14.6, &
1017 ' phi:', e14.6,
' dphi:', e14.6)
1024 CALL cg_value(phi, xtemp, n)
1026 IF ( ngrow + nshrink == 0 ) fquad = min(phi, f0)
1027 IF ( phi <= fquad )
THEN 1028 IF ( cg_wolfe(alpha, phi, dphi) )
RETURN 1031 DO iter = 1, nsecant
1032 IF ( printlevel )
THEN 1033 WRITE (*, 40) a, b, dphia, dphib
1034 40
FORMAT (
'secant, a:', e14.6,
' b:', e14.6, &
1035 ' da:', e14.6,
' db:', e14.6)
1037 width = gamma*(b - a)
1038 IF ( -dphia <= dphib )
THEN 1039 alpha = a - (a-b)*(dphia/(dphia-dphib))
1041 alpha = b - (a-b)*(dphib/(dphia-dphib))
1048 flag = cg_update(a, dphia, b, dphib, alpha, phi, &
1049 dphi, x, xtemp, d, gtemp, cg_value, cg_grad)
1050 IF ( flag > 0 )
THEN 1052 ELSE IF ( flag == 0 )
THEN 1054 IF ( dphi > da0 )
THEN 1055 alpha = c - (c-a0)*(dphi/(dphi-da0))
1060 IF ( dphi < db0 )
THEN 1061 alpha = c - (c-b0)*(dphi/(dphi-db0))
1066 IF ( (alpha > a) .AND. (alpha < b) )
THEN 1067 IF ( printlevel )
WRITE (*, *)
"2nd secant" 1068 flag = cg_update(a, dphia, b, dphib, alpha, phi, &
1069 dphi, x, xtemp, d, gtemp, cg_value, cg_grad)
1070 IF ( flag > 0 )
RETURN 1076 IF ( (b-a) >= width )
THEN 1078 IF ( printlevel )
WRITE (*, *)
"bisection" 1079 flag = cg_update(a, dphia, b, dphib, alpha, phi, &
1080 dphi, x, xtemp, d, gtemp, cg_value, cg_grad)
1081 IF ( flag > 0 )
RETURN 1091 END SUBROUTINE cg_line
1111 SUBROUTINE cg_linew (alpha, phi, dphi, dphi0, x, xtemp, d, gtemp, cg_value, cg_grad)
1116 real,
INTENT(IN OUT) :: alpha
1117 real,
INTENT(IN OUT) :: phi
1118 real,
INTENT(OUT) :: dphi
1119 real,
INTENT(IN) :: dphi0
1120 real,
INTENT(IN OUT) :: x(:)
1121 real,
INTENT(IN OUT) :: xtemp(:)
1122 real,
INTENT(IN OUT) :: d(:)
1123 real,
INTENT(IN OUT) :: gtemp(:)
1127 SUBROUTINE cg_value (f, x, n)
1129 INTEGER,
INTENT(IN) :: n
1130 REAL,
INTENT(IN) :: x(:)
1131 REAL,
INTENT(IN OUT) :: f
1132 END SUBROUTINE cg_value
1134 SUBROUTINE cg_grad (g, x, n)
1136 INTEGER,
INTENT(IN) :: n
1137 REAL,
INTENT(IN) :: x(:)
1138 REAL,
INTENT(IN OUT) :: g(:)
1139 END SUBROUTINE cg_grad
1143 real :: a, dpsia, b, dpsib, c, &
1144 a0, da0, b0, db0, width, fquad, psi, dpsi
1146 INTEGER :: ngrow, nshrink, iter, flag
1150 CALL cg_step (xtemp, x, d, alpha)
1151 CALL cg_grad (gtemp, xtemp, n)
1153 dphi = cg_dot(gtemp, d)
1154 dpsi = dphi - wolfe_hi
1160 dpsia = dphi0 - wolfe_hi
1163 DO WHILE ( dpsi < zero )
1164 CALL cg_value(phi, xtemp, n)
1165 psi = phi - alpha*wolfe_hi
1172 IF ( ngrow == 0 ) fquad = min(phi, f0)
1173 IF ( phi <= fquad )
THEN 1174 IF ( printlevel )
THEN 1175 WRITE (*, 10) alpha, phi, fquad
1176 10
FORMAT (
'alpha:', e14.6,
' phi:', e14.6, &
1179 IF ( cg_wolfe(alpha, phi, dphi) )
RETURN 1182 IF ( psi <= fpert )
THEN 1192 nshrink = nshrink + 1
1193 IF ( nshrink > nexpand )
THEN 1197 CALL cg_step (xtemp, x, d, alpha)
1198 CALL cg_grad (gtemp, xtemp, n)
1200 dphi = cg_dot(gtemp, d)
1201 dpsi = dphi - wolfe_hi
1202 IF ( dpsi >= zero )
GO TO 100
1203 CALL cg_value (phi, xtemp, n)
1204 psi = phi - alpha*wolfe_hi
1206 IF ( printlevel )
THEN 1207 WRITE (6, 20) a, b, alpha, phi, dphi
1208 20
FORMAT (
'contract, a:', e14.6, &
1209 ' b:', e14.6,
' alpha:', e14.6,
' phi:', e14.6,
' dphi:', e14.6)
1211 IF ( quadok .AND. (phi <= fquad) )
THEN 1212 IF ( cg_wolfe(alpha, phi, dphi) )
RETURN 1214 IF ( psi <= fpert )
THEN 1226 IF ( ngrow > nexpand )
THEN 1231 CALL cg_step (xtemp, x, d, alpha)
1232 CALL cg_grad (gtemp, xtemp, n)
1234 dphi = cg_dot(gtemp, d)
1235 dpsi = dphi - wolfe_hi
1236 IF ( printlevel )
THEN 1237 WRITE (*, 30) a, alpha, phi, dphi
1238 30
FORMAT (
'expand, a:', e14.6,
' alpha:', e14.6, &
1239 ' phi:', e14.6,
' dphi:', e14.6)
1240 WRITE (6, *)
"expand, alpha:", alpha,
"dphi:", dphi
1247 CALL cg_value(phi, xtemp, n)
1249 IF ( ngrow + nshrink == 0 ) fquad = min(phi, f0)
1250 IF ( phi <= fquad )
THEN 1251 IF ( cg_wolfe(alpha, phi, dphi) )
RETURN 1254 DO iter = 1, nsecant
1255 IF ( printlevel )
THEN 1256 WRITE (*, 40) a, b, dpsia, dpsib
1257 40
FORMAT (
'secant, a:', e14.6,
' b:', e14.6, &
1258 ' da:', e14.6,
' db:', e14.6)
1260 width = gamma*(b - a)
1261 IF ( -dpsia <= dpsib )
THEN 1262 alpha = a - (a-b)*(dpsia/(dpsia-dpsib))
1264 alpha = b - (a-b)*(dpsib/(dpsia-dpsib))
1271 flag = cg_updatew(a, dpsia, b, dpsib, alpha, &
1272 phi, dphi, dpsi, x, xtemp, d, gtemp, cg_value, cg_grad)
1273 IF ( flag > 0 )
THEN 1275 ELSE IF ( flag == 0 )
THEN 1277 IF ( dpsi > da0 )
THEN 1278 alpha = c - (c-a0)*(dpsi/(dpsi-da0))
1283 IF ( dpsi < db0 )
THEN 1284 alpha = c -(c-b0)*(dpsi/(dpsi-db0))
1289 IF ( (alpha > a) .AND. (alpha < b) )
THEN 1290 IF ( printlevel )
WRITE (*, *)
"2nd secant" 1291 flag = cg_updatew(a, dpsia, b, dpsib, alpha, &
1292 phi, dphi, dpsi, x, xtemp, d, gtemp, cg_value, cg_grad)
1293 IF ( flag > 0 )
RETURN 1299 IF ( (b-a) >= width )
THEN 1301 IF ( printlevel )
WRITE (*, *)
"bisection" 1302 flag = cg_updatew(a, dpsia, b, dpsib, alpha, &
1303 phi, dphi, dpsi, x, xtemp, d, gtemp, cg_value, cg_grad)
1304 IF ( flag > 0 )
RETURN 1314 END SUBROUTINE cg_linew
1337 INTEGER FUNCTION cg_update (a, dphia, b, dphib, alpha, phi, &
1338 dphi, x, xtemp, d, gtemp, cg_value, cg_grad)
1343 real,
INTENT(OUT) :: a
1344 real,
INTENT(OUT) :: dphia
1345 real,
INTENT(OUT) :: b
1346 real,
INTENT(OUT) :: dphib
1347 real,
INTENT(IN OUT) :: alpha
1348 real,
INTENT(IN OUT) :: phi
1349 real,
INTENT(OUT) :: dphi
1350 real,
INTENT(IN OUT) :: x(:)
1351 real,
INTENT(IN OUT) :: xtemp(:)
1352 real,
INTENT(IN OUT) :: d(:)
1353 real,
INTENT(IN OUT) :: gtemp(:)
1357 SUBROUTINE cg_value (f, x, n)
1359 INTEGER,
INTENT(IN) :: n
1360 REAL,
INTENT(IN) :: x(:)
1361 REAL,
INTENT(IN OUT) :: f
1362 END SUBROUTINE cg_value
1364 SUBROUTINE cg_grad (g, x, n)
1366 INTEGER,
INTENT(IN) :: n
1367 REAL,
INTENT(IN) :: x(:)
1368 REAL,
INTENT(IN OUT) :: g(:)
1369 END SUBROUTINE cg_grad
1377 CALL cg_step (xtemp, x, d, alpha)
1378 CALL cg_value (phi, xtemp, n)
1380 CALL cg_grad (gtemp, xtemp, n)
1382 dphi = cg_dot(gtemp, d)
1383 IF ( printlevel )
THEN 1384 WRITE (*, 10) alpha, phi, dphi
1385 10
FORMAT (
'update alpha:', e14.6,
' phi:', e14.6,
' dphi:', e14.6)
1388 IF ( cg_wolfe(alpha, phi, dphi) )
THEN 1392 IF ( dphi >= zero )
THEN 1397 IF ( phi <= fpert )
THEN 1407 nshrink = nshrink + 1
1408 IF ( nshrink > nexpand )
THEN 1413 CALL cg_step (xtemp, x, d, alpha)
1414 CALL cg_grad (gtemp, xtemp, n)
1416 dphi = cg_dot(gtemp, d)
1417 CALL cg_value (phi, xtemp, n)
1419 IF ( printlevel )
THEN 1420 WRITE (6, 20) a, alpha, phi, dphi
1421 20
FORMAT (
'contract, a:', e14.6,
' alpha:', e14.6, &
1422 ' phi:', e14.6,
' dphi:', e14.6)
1424 IF ( cg_wolfe(alpha, phi, dphi) )
THEN 1428 IF ( dphi >= zero )
THEN 1433 IF ( phi <= fpert )
THEN 1434 IF ( printlevel )
THEN 1435 WRITE (6, *)
"update a:", alpha,
"dphia:", dphi
1446 IF ( printlevel )
THEN 1447 WRITE (*, 200) a, b, dphia, dphib, cg_update
1448 200
FORMAT (
'UP a:', e14.6,
' b:', e14.6, &
1449 ' da:', e14.6,
' db:', e14.6,
' up:', i2)
1452 END FUNCTION cg_update
1480 INTEGER FUNCTION cg_updatew (a, dpsia, b, dpsib, alpha, phi, dphi, &
1481 dpsi, x, xtemp, d, gtemp, cg_value, cg_grad)
1486 real,
INTENT(OUT) :: a
1487 real,
INTENT(OUT) :: dpsia
1488 real,
INTENT(OUT) :: b
1489 real,
INTENT(OUT) :: dpsib
1490 real,
INTENT(IN OUT) :: alpha
1491 real,
INTENT(IN OUT) :: phi
1492 real,
INTENT(OUT) :: dphi
1493 real,
INTENT(OUT) :: dpsi
1494 real,
INTENT(IN OUT) :: x(:)
1495 real,
INTENT(IN OUT) :: xtemp(:)
1496 real,
INTENT(IN OUT) :: d(:)
1497 real,
INTENT(IN OUT) :: gtemp(:)
1501 SUBROUTINE cg_value (f, x, n)
1503 INTEGER,
INTENT(IN) :: n
1504 REAL,
INTENT(IN) :: x(:)
1505 REAL,
INTENT(IN OUT) :: f
1506 END SUBROUTINE cg_value
1508 SUBROUTINE cg_grad (g, x, n)
1510 INTEGER,
INTENT(IN) :: n
1511 REAL,
INTENT(IN) :: x(:)
1512 REAL,
INTENT(IN OUT) :: g(:)
1513 END SUBROUTINE cg_grad
1522 CALL cg_step (xtemp, x, d, alpha)
1523 CALL cg_value (phi, xtemp, n)
1524 psi = phi - alpha*wolfe_hi
1526 CALL cg_grad (gtemp, xtemp, n)
1528 dphi = cg_dot(gtemp, d)
1529 dpsi = dphi - wolfe_hi
1530 IF ( printlevel )
THEN 1531 WRITE (*, 10) alpha, psi, dpsi
1532 10
FORMAT (
'update alpha:', e14.6,
' psi:', e14.6,
' dpsi:', e14.6)
1535 IF ( cg_wolfe(alpha, phi, dphi) )
THEN 1539 IF ( dpsi >= zero )
THEN 1544 IF ( psi <= fpert )
THEN 1554 nshrink = nshrink + 1
1555 IF ( nshrink > nexpand )
THEN 1560 CALL cg_step (xtemp, x, d, alpha)
1561 CALL cg_grad (gtemp, xtemp, n)
1563 dphi = cg_dot(gtemp, d)
1564 dpsi = dphi - wolfe_hi
1565 CALL cg_value (phi, xtemp, n)
1566 psi = phi - alpha*wolfe_hi
1568 IF ( printlevel )
THEN 1569 WRITE (6, 20) a, alpha, phi, dphi
1570 20
FORMAT (
'contract, a:', e14.6,
' alpha:', e14.6, &
1571 ' phi:', e14.6,
' dphi:', e14.6)
1573 IF ( cg_wolfe(alpha, phi, dphi) )
THEN 1577 IF ( dpsi >= zero )
THEN 1582 IF ( psi <= fpert )
THEN 1583 IF ( printlevel )
THEN 1584 WRITE (6, *)
"update a:", alpha,
"dpsia:", dpsi
1595 IF ( printlevel )
THEN 1596 WRITE (*, 200) a, b, dpsia, dpsib, cg_updatew
1597 200
FORMAT (
'UP a:', e14.6,
' b:', e14.6, &
1598 ' da:', e14.6,
' db:', e14.6,
' up:', i2)
1601 END FUNCTION cg_updatew
1620 END MODULE cg_minimization