MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
cg_descent.f90
1 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2 ! ________________________________________________________________
3 ! | A conjugate gradient method with guaranteed descent |
4 ! | |
5 ! | Version 1.1 (December 10, 2004) |
6 ! | Version 1.2 (June 4, 2005) |
7 ! | Version 1.3 (October 6, 2005) |
8 ! | Version 1.4 (November 14, 2005) |
9 ! | |
10 ! | William W. Hager and Hongchao Zhang |
11 ! | hager@math.ufl.edu hzhang@math.ufl.edu |
12 ! | Department of Mathematics |
13 ! | University of Florida |
14 ! | Gainesville, Florida 32611 USA |
15 ! | 352-392-0281 x 244 |
16 ! | |
17 ! | Copyright 2004 by William W. Hager |
18 ! | |
19 ! |This program is free software; you can redistribute it and/or |
20 ! |modify it under the terms of the GNU General Public License as |
21 ! |published by the Free Software Foundation; either version 2 of |
22 ! |the License, or (at your option) any later version. |
23 ! |This program is distributed in the hope that it will be useful, |
24 ! |but WITHOUT ANY WARRANTY; without even the implied warranty of |
25 ! |MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
26 ! |GNU General Public License for more details. |
27 ! | |
28 ! |You should have received a copy of the GNU General Public |
29 ! |License along with this program; if not, write to the Free |
30 ! |Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, |
31 ! |MA 02110-1301 USA |
32 ! | |
33 ! | http://www.math.ufl.edu/~hager/papers/CG |
34 ! | |
35 ! | INPUT: |
36 ! | |
37 ! |(double) grad_tol-- StopRule = T: |g|_infty <= max (grad_tol, |
38 ! | StopFac*initial |g|_infty) [default] |
39 ! | StopRule = F: |g|_infty <= grad_tol(1+|f|) |
40 ! | |
41 ! |(double) x --starting guess (length n) |
42 ! | |
43 ! |(int) dim --problem dimension (also denoted n) |
44 ! | |
45 ! | cg_value--name of cost evaluation subroutine |
46 ! | (external in main program, cg_value(f, x, n) |
47 ! | puts value of cost function at x in f |
48 ! | f is double precision scalar and x is |
49 ! | double precision array of length n) |
50 ! | |
51 ! | cg_grad --name gradient evaluation subroutine |
52 ! | (external in main program, cg_grad (g, x, n) |
53 ! | puts gradient at x in g, g and x are |
54 ! | double precision arrays of length n) |
55 ! | |
56 ! |(double) gnorm --if the parameter Step in cg_descent.parm is |
57 ! | .true., then gnorm contains the initial step |
58 ! | used at iteration 0 in the line search |
59 ! | |
60 ! |(double) d --direction (work array, length n) |
61 ! | |
62 ! |(double) g --gradient (work array, length n) |
63 ! | |
64 ! |(double) xtemp --work array (work array, length n) |
65 ! | |
66 ! |(double) gtemp --work array (work array, length n) |
67 ! | |
68 ! | OUTPUT: |
69 ! | |
70 ! |(int) status -- 0 (convergence tolerance satisfied) |
71 ! | 1 (change in func <= feps*|f|) |
72 ! | 2 (total iterations exceeded maxit) |
73 ! | 3 (slope always negative in line search) |
74 ! | 4 (number secant iterations exceed nsecant) |
75 ! | 5 (search direction not a descent direction)|
76 ! | 6 (line search fails in initial interval) |
77 ! | 7 (line search fails during bisection) |
78 ! | 8 (line search fails during interval update)|
79 ! | |
80 ! |(double) gnorm --max abs component of gradient |
81 ! | |
82 ! |(double) f --function value at solution |
83 ! | |
84 ! |(double) x --solution (length n) |
85 ! | |
86 ! |(int) iter --number of iterations |
87 ! | |
88 ! |(int) nfunc --number of function evaluations |
89 ! | |
90 ! |(int) ngrad --number of gradient evaluations |
91 ! | |
92 ! |Note: The file cg_descent.parm must be placed in the directory |
93 ! | where the code is run |
94 ! |________________________________________________________________|
95 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
96 
97 MODULE cg_minimization
98 
99 
100  implicit none
101  save
102 
103  integer :: n, n5, n6, nf, ng
104  integer :: info
105  integer :: nrestart, nexpand, nsecant
106  integer :: maxit
107  real :: delta, sigma, eps
108  real :: gamma, rho, tol
109  real :: eta, fpert, f0
110  real :: ck, qdecay
111  real :: wolfe_hi, wolfe_lo, awolfe_hi
112  real :: quadcutoff, stopfac, awolfefac
113  real :: zero, feps
114  real :: psi0, psi1, psi2
115  logical :: pertrule, quadok
116  logical :: quadstep, printlevel
117  logical :: printfinal, stoprule
118  logical :: awolfe
119  logical :: step
120  logical :: debug
121 
122  real :: f_best
123  real, dimension(:), allocatable :: x_best
124 
125  PRIVATE
126  PUBLIC :: cg_descent
127 
128 CONTAINS
129 
130 
131  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
132  !
133  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
134 
135  SUBROUTINE cg_descent (grad_tol, x, dim, cg_value, cg_grad, &
136  STATUS, gnorm, f, iter, nfunc, ngrad, d, g, xtemp, gtemp)
137 
138  implicit none
139 
140 
141  INTEGER, INTENT(IN) :: dim ! (= n) dimensionality of problem
142  REAL, INTENT(IN) :: grad_tol
143  REAL, INTENT(IN OUT) :: x(dim) ! initial guess / for output: solution (length n)
144  INTEGER, INTENT(OUT) :: status
145  REAL, INTENT(IN OUT) :: gnorm ! max abs component of gradient
146  REAL, INTENT(IN OUT) :: f ! function value at solution
147  INTEGER, INTENT(OUT) :: iter ! number of iterations
148  INTEGER, INTENT(OUT) :: nfunc ! number of function evaluations
149  INTEGER, INTENT(OUT) :: ngrad ! number of gradient evaluations
150  REAL, INTENT(IN OUT) :: d(dim) ! direction (work array, length n)
151  REAL, INTENT(IN OUT) :: g(dim) ! gradient (work array, length n)
152  REAL, INTENT(IN OUT) :: xtemp(dim) ! work array (length n)
153  REAL, INTENT(IN OUT) :: gtemp(dim) ! work array (length n)
154 
155  INTERFACE
156 
157  SUBROUTINE cg_value (f, x, n)
158  IMPLICIT NONE
159  INTEGER, INTENT(IN) :: n
160  REAL, INTENT(IN) :: x(:)
161  REAL, INTENT(IN OUT) :: f
162  END SUBROUTINE cg_value
163 
164  SUBROUTINE cg_grad (g, x, n)
165  IMPLICIT NONE
166  INTEGER, INTENT(IN) :: n
167  REAL, INTENT(IN) :: x(:)
168  REAL, INTENT(IN OUT) :: g(:)
169  END SUBROUTINE cg_grad
170 
171  END INTERFACE
172 
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
178  real :: dphi, dphi0
179  real :: alpha, talpha
180  real :: yk, yk2, ykgk, dkyk
181  real :: beta
182  !---------------------------------------------------------------------------
183 
184  allocate( x_best( dim ) )
185  f_best = 1.e30
186 
187  ! initialize the parameters
188 
189  CALL cg_init (grad_tol, dim)
190  IF ( step ) THEN
191  alpha = gnorm
192  END IF
193  delta2 = 2.0 * delta - 1.0
194  eta_sq = eta * eta
195  iter = 0
196  ck = 0
197  qk = 0
198 
199  ! initial function and gradient evaluations, initial direction
200 
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
204  nf = nf + 1
205  CALL cg_grad (g, x, n)
206  ng = ng + 1
207  f0 = f + f
208  gnorm = zero
209  xnorm = zero
210  gnorm2 = zero
211  DO i = 1, n5
212  xnorm = max(xnorm, abs(x(i)))
213  t = g(i)
214  d(i) = -t
215  gnorm = max(gnorm, abs(t))
216  gnorm2 = gnorm2 + t*t
217  END DO
218  DO i = n6, n, 5
219  xnorm = max(xnorm, abs(x(i)))
220  t = g(i)
221  gnorm = max(gnorm, abs(t))
222  d(i) = -t
223  j = i + 1
224  t1 = g(j)
225  d(j) = -t1
226  gnorm = max(gnorm, abs(t1))
227  xnorm = max(xnorm, abs(x(j)))
228  j = i + 2
229  t2 = g(j)
230  d(j) = -t2
231  gnorm = max(gnorm, abs(t2))
232  xnorm = max(xnorm, abs(x(j)))
233  j = i + 3
234  t3 = g(j)
235  d(j) = -t3
236  gnorm = max(gnorm, abs(t3))
237  xnorm = max(xnorm, abs(x(j)))
238  j = i + 4
239  t4 = g(j)
240  d(j) = -t4
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
244  END DO
245 
246  IF ( stoprule ) THEN
247  tol = max(gnorm*stopfac, tol)
248  END IF
249 
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)
254  END IF
255 
256  IF ( cg_tol(f, gnorm) ) GO TO 100
257 
258  dphi0 = -gnorm2
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
264  ELSE
265  alpha = 1.0
266  END IF
267  END IF
268  END IF
269 
270  ! start the conjugate gradient iteration
271 
272 
273  ! alpha starts as old step, ends as initial step for next iteration
274  ! f is function value for alpha = 0
275  ! QuadOK = .true. means that a quadratic step was taken
276 
277 
278  DO iter = 1, maxit
279  quadok = .false.
280  alpha = psi2*alpha
281  IF ( quadstep ) THEN
282  IF ( f /= zero ) THEN
283  t = abs((f-f0)/f)
284  ELSE
285  t = 1.0
286  END IF
287  IF ( t > quadcutoff ) THEN
288  talpha = psi1*alpha
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
293  nf = nf + 1
294  IF ( ftemp < f ) THEN
295  denom = 2.0d0*(((ftemp-f)/talpha)-dphi0)
296  IF ( denom > zero ) THEN
297  quadok = .true.
298  alpha = -dphi0*talpha/denom
299  END IF
300  END IF
301  END IF
302  END IF
303  f0 = f
304 
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)
309  END IF
310 
311  ! parameters in Wolfe and approximiate Wolfe conditions, and in update
312 
313  qk = qdecay*qk + 1.
314  ck = ck + (abs(f) - ck)/qk
315 
316  IF ( pertrule ) THEN
317  fpert = f + eps*ck
318  ELSE
319  fpert = f + eps
320  END IF
321 
322  wolfe_hi = delta*dphi0
323  wolfe_lo = sigma*dphi0
324  awolfe_hi = delta2*dphi0
325  IF ( awolfe ) THEN
326  CALL cg_line (alpha, f, dphi, dphi0, x, xtemp, d, gtemp, &
327  cg_value, cg_grad)
328  ELSE
329  CALL cg_linew (alpha, f, dphi, dphi0, x, xtemp, d, gtemp, &
330  cg_value, cg_grad)
331  END IF
332 
333  IF ( info > 0 ) GO TO 100
334 
335  ! Test for convergence to within machine epsilon
336  ! (set feps to zero to remove this test)
337 
338  IF ( -alpha*dphi0 <= feps*abs(f) ) THEN
339  info = 1
340  GO TO 100
341  END IF
342 
343  ! compute beta, yk2, gnorm, gnorm2, dnorm2, update x and g,
344 
345  IF ( mod(iter, nrestart) /= 0 ) THEN
346  gnorm = zero
347  dnorm2 = zero
348  yk2 = zero
349  ykgk = zero
350  DO i = 1, n5
351  x(i) = xtemp(i)
352  t = gtemp(i)
353  yk = t - g(i)
354  yk2 = yk2 + yk**2
355  ykgk = ykgk + yk*t
356  g(i) = t
357  gnorm = max(gnorm, abs(t))
358  dnorm2 = dnorm2 + d(i)**2
359  END DO
360  DO i = n6, n, 5
361  x(i) = xtemp(i)
362  t = gtemp(i)
363  yk = t - g(i)
364  yk2 = yk2 + yk**2
365  ykgk = ykgk + yk*t
366  i1 = i + 1
367  x(i1) = xtemp(i1)
368  t1 = gtemp(i1)
369  i2 = i + 2
370  x(i2) = xtemp(i2)
371  t2 = gtemp(i2)
372  i3 = i + 3
373  x(i3) = xtemp(i3)
374  t3 = gtemp(i3)
375  i4 = i + 4
376  x(i4) = xtemp(i4)
377  t4 = gtemp(i4)
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
382  g(i) = t
383  gnorm = max(gnorm, abs(t))
384  g(i1) = t1
385  gnorm = max(gnorm, abs(t1))
386  g(i2) = t2
387  gnorm = max(gnorm, abs(t2))
388  g(i3) = t3
389  gnorm = max(gnorm, abs(t3))
390  g(i4) = t4
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
394  END DO
395  IF ( cg_tol(f, gnorm) ) GO TO 100
396  dkyk = dphi - dphi0
397  beta = (ykgk - 2.0*dphi*yk2/dkyk)/dkyk
398 
399  ! faster: initialize dnorm2 = gnorm2 at start, then
400  ! dnorm2 = gnorm2 + beta**2*dnorm2 - 2.0*beta*dphi
401  ! gnorm2 = ||g_{k+1}||^2
402  ! dnorm2 = ||d_{k+1}||^2
403  ! dpi = g_{k+1}' d_k
404 
405  beta = max(beta, -1.0/sqrt(min(eta_sq, gnorm2)*dnorm2))
406 
407  ! update search direction d = -g + beta*dold
408 
409  gnorm2 = zero
410  DO i = 1, n5
411  t = g(i)
412  d(i) = -t + beta*d(i)
413  gnorm2 = gnorm2 + t*t
414  END DO
415  DO i = n6, n, 5
416  d(i) = -g(i) + beta*d(i)
417  i1 = i + 1
418  d(i1) = -g(i1) + beta*d(i1)
419  i2 = i + 2
420  d(i2) = -g(i2) + beta*d(i2)
421  i3 = i + 3
422  d(i3) = -g(i3) + beta*d(i3)
423  i4 = i + 4
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
427  END DO
428  dphi0 = -gnorm2 + beta*dphi
429 
430  ELSE
431 
432  ! search direction d = -g
433 
434  IF ( printlevel ) THEN
435  WRITE (*, *) "RESTART CG"
436  END IF
437  gnorm = zero
438  gnorm2 = zero
439  DO i = 1, n5
440  x(i) = xtemp(i)
441  t = gtemp(i)
442  g(i) = t
443  d(i) = -t
444  gnorm = max(gnorm, abs(t))
445  gnorm2 = gnorm2 + t*t
446  END DO
447  DO i = n6, n, 5
448  x(i) = xtemp(i)
449  t = gtemp(i)
450  g(i) = t
451  d(i) = -t
452  gnorm = max(gnorm, abs(t))
453  j = i + 1
454  x(j) = xtemp(j)
455  t1 = gtemp(j)
456  g(j) = t1
457  d(j) = -t1
458  gnorm = max(gnorm, abs(t1))
459  j = i + 2
460  x(j) = xtemp(j)
461  t2 = gtemp(j)
462  g(j) = t2
463  d(j) = -t2
464  gnorm = max(gnorm, abs(t2))
465  j = i + 3
466  x(j) = xtemp(j)
467  t3 = gtemp(j)
468  g(j) = t3
469  d(j) = -t3
470  gnorm = max(gnorm, abs(t3))
471  j = i + 4
472  x(j) = xtemp(j)
473  t4 = gtemp(j)
474  g(j) = t4
475  d(j) = -t4
476  gnorm = max(gnorm, abs(t4))
477  gnorm2 = gnorm2 + t*t + t1*t1 + t2*t2 + t3*t3 + t4*t4
478  END DO
479  IF ( cg_tol(f, gnorm) ) GO TO 100
480  dphi0 = -gnorm2
481  END IF
482  IF ( .NOT.awolfe ) THEN
483  IF ( abs(f-f0) < awolfefac*ck ) THEN
484  awolfe = .true.
485  END IF
486  END IF
487 
488  IF ( printlevel ) THEN
489  WRITE (*, 10) iter, f, gnorm, awolfe
490  END IF
491 
492  IF ( debug ) THEN
493  IF ( f > f0 + 1.e-10*ck ) THEN
494  WRITE (*, 270)
495  WRITE (*, 50) f, f0
496 50 FORMAT (' new value:', e30.16, 'old value:', e30.16)
497  stop
498  END IF
499  END IF
500 
501  IF ( dphi0 > zero ) THEN
502  info = 5
503  GO TO 100
504  END IF
505  END DO
506  info = 2
507 100 nfunc = nf
508  ngrad = ng
509  status = info
510  IF ( info > 2 ) THEN
511  gnorm = zero
512  DO i = 1, n
513  x(i) = xtemp(i)
514  g(i) = gtemp(i)
515  gnorm = max(gnorm, abs(g(i)))
516  END DO
517  END IF
518  IF ( printfinal ) THEN
519  WRITE (6, *) 'Termination status:', status
520  IF ( status == 0 ) THEN
521  WRITE (6, 200)
522  ELSE IF ( status == 1 ) THEN
523  WRITE (6, 210)
524  ELSE IF ( status == 2 ) THEN
525  WRITE (6, 220) maxit
526  WRITE (6, 300)
527  WRITE (6, 400) grad_tol
528  ELSE IF ( status == 3 ) THEN
529  WRITE (6, 230)
530  WRITE (6, 300)
531  WRITE (6, 430)
532  WRITE (6, 410)
533  ELSE IF ( status == 4 ) THEN
534  WRITE (6, 240)
535  WRITE (6, 300)
536  WRITE (6, 400) grad_tol
537  ELSE IF ( status == 5 ) THEN
538  WRITE (6, 250)
539  ELSE IF ( status == 6 ) THEN
540  WRITE (6, 260)
541  WRITE (6, 300)
542  WRITE (6, 400) grad_tol
543  WRITE (6, 410)
544  WRITE (6, 420)
545  ELSE IF ( status == 7 ) THEN
546  WRITE (6, 260)
547  WRITE (6, 300)
548  WRITE (6, 400) grad_tol
549  ELSE IF ( status == 8 ) THEN
550  WRITE (6, 260)
551  WRITE (6, 300)
552  WRITE (6, 400) grad_tol
553  WRITE (6, 410)
554  WRITE (6, 420)
555  END IF
556  WRITE (6, 500) gnorm
557  WRITE (6, *) 'function value:', f
558  WRITE (6, *) 'cg iterations:', iter
559  WRITE (6, *) 'function evaluations:', nfunc
560  WRITE (6, *) 'gradient evaluations:', ngrad
561  END IF
562 
563  if ( f_best < f ) then
564  f = f_best
565  x(1:dim) = x_best(1:dim)
566  end if
567 
568  deallocate( x_best )
569 
570  RETURN
571 
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)
587 
588  END SUBROUTINE cg_descent
589 
590 
591 
592 
593  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
594  ! PARAMETERS:
595 
596  ! delta - range (0, .5), used in the Wolfe conditions
597  ! sigma - range [delta, 1), used in the Wolfe conditions
598  ! eps - range [0, infty), used to compute line search perturbation
599  ! gamma - range (0,1), determines when to perform bisection step
600  ! rho - range (1, infty), growth factor when finding initial interval
601  ! eta - range (0, infty), used in lower bound for beta
602  ! psi0 - range (0, 1), factor used in very initial starting guess
603  ! psi1 - range (0, 1), factor previous step multiplied by in QuadStep
604  ! psi2 - range (1, infty), factor previous step is multipled by for startup
605  ! QuadCutOff - perform QuadStep if relative change in f > QuadCutOff
606  ! StopFac - used in StopRule
607  ! AWolfeFac - used to decide when to switch from Wolfe to AWolfe
608  ! restart_fac - range (0, infty) restart cg when iter = n*restart
609  ! maxit_fac - range (0, infty) terminate in maxit = maxit_fac*n iterations
610  ! feps - stop when -alpha*dphi0 (est. change in value) <= feps*|f|
611  ! (feps = 0 removes this test, example: feps = eps*1.e-5
612  ! where eps is machine epsilon)
613  ! tol - range (0, infty), convergence tolerance
614  ! nexpand - range [0, infty), number of grow/shrink allowed in bracket
615  ! nsecant - range [0, infty), maximum number of secant steps
616  ! PertRule - gives the rule used for the perturbation in f
617  ! F => fpert = eps
618  ! T => fpert = eps*Ck, Ck is an average of prior |f|
619  ! Ck is an average of prior |f|
620  ! QuadStep- .true. (use quadratic step) .false. (no quadratic step)
621  ! PrintLevel- .false. (no printout) .true. (print intermediate results)
622  ! PrintFinal- .false. (no printout) .true. (print messages, final error)
623  ! StopRule - .true. (max abs grad <= max (tol, StopFac*initial abs grad))
624  ! .false. (... <= tol*(1+|f|))
625  ! AWolfe - .false. (use standard Wolfe initially)
626  ! - .true. (use approximate + standard Wolfe)
627  ! Step - .false. (program computing starting step at iteration 0)
628  ! - .true. (user provides starting step in gnorm argument of cg_descent
629  ! debug - .false. (no debugging)
630  ! - .true. (check that function values do not increase)
631  ! info - same as status
632 
633  ! DEFAULT PARAMETER VALUES:
634 
635  ! delta : 0.1
636  ! sigma : 0.9
637  ! eps : 1.e-6
638  ! gamma : 0.66
639  ! rho : 5.0
640  ! restart: 1.0
641  ! eta : 0.01
642  ! psi0 : 0.01
643  ! psi1 : 0.1
644  ! psi2 : 2.0
645  ! QuadCutOff: 1.E-12
646  ! StopFac: 0.0
647  ! AWolfeFac: 1.E-3
648  ! tol : grad_tol
649  ! nrestart: n (restart_fac = 1)
650  ! maxit : 500*n (maxit_fac = 500)
651  ! feps : 0.0
652  ! Qdecay : 0.7
653  ! nexpand: 50
654  ! nsecant: 50
655  ! PertRule: .true.
656  ! QuadStep: .true.
657  ! PrintLevel: .false.
658  ! PrintFinal: .true.
659  ! StopRule: .true.
660  ! AWolfe: .false.
661  ! Step: .false.
662  ! debug: .false.
663  ! info : 0
664  ! feps : 0.0
665 
666 
667  ! (double) grad_tol-- used in stopping rule
668  ! (int) dim --problem dimension (also denoted n)
669  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
670 
671  SUBROUTINE cg_init (grad_tol, dim)
672 
673  ! use define_cg_variables
674  implicit none
675 
676  real, INTENT(IN) :: grad_tol
677  integer, INTENT(IN) :: dim
678 
679  real :: restart_fac, maxit_fac
680  !---------------------------------------------------------------------------
681 
682  n = dim
683  tol = grad_tol
684  OPEN (10, file='cg_descent_f.parm')
685  READ (10, *) delta
686  READ (10, *) sigma
687  READ (10, *) eps
688  READ (10, *) gamma
689  READ (10, *) rho
690  READ (10, *) eta
691  READ (10, *) psi0
692  READ (10, *) psi1
693  READ (10, *) psi2
694  READ (10, *) quadcutoff
695  READ (10, *) stopfac
696  READ (10, *) awolfefac
697  READ (10, *) restart_fac
698  READ (10, *) maxit_fac
699  READ (10, *) feps
700  READ (10, *) qdecay
701  READ (10, *) nexpand
702  READ (10, *) nsecant
703  READ (10, *) pertrule
704  READ (10, *) quadstep
705  READ (10, *) printlevel
706  READ (10, *) printfinal
707  READ (10, *) stoprule
708  READ (10, *) awolfe
709  READ (10, *) step
710  READ (10, *) debug
711  nrestart = n*restart_fac
712  maxit = n*maxit_fac
713  zero = 0.0
714  info = 0
715  n5 = mod(n, 5)
716  n6 = n5 + 1
717  nf = 0
718  ng = 0
719  CLOSE (10)
720 
721  END SUBROUTINE cg_init
722 
723 
724 
725  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
726  ! check whether the Wolfe or the approximate Wolfe conditions
727  ! are satisfied
728 
729  ! (double) alpha -- stepsize
730  ! (double) f -- function value associated with stepsize alpha
731  ! (double) dphi -- derivative value associated with stepsize alpha
732  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
733 
734  LOGICAL FUNCTION cg_wolfe (alpha, f, dphi)
735 
736  ! use define_cg_variables
737  implicit none
738 
739  real, INTENT(IN) :: alpha
740  real, INTENT(IN) :: f
741  real, INTENT(IN) :: dphi
742  !---------------------------------------------------------------------------
743 
744  IF ( dphi >= wolfe_lo ) THEN
745 
746  ! test original Wolfe conditions
747 
748  IF ( f-f0 <= alpha*wolfe_hi ) THEN
749  cg_wolfe = .true.
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)
754  END IF
755  RETURN
756 
757  ! test approximate Wolfe conditions
758 
759  ELSE IF ( awolfe ) THEN
760  IF ( (f <= fpert).AND.(dphi <= awolfe_hi) ) THEN
761  cg_wolfe = .true.
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)
766  END IF
767  RETURN
768  END IF
769  END IF
770  END IF
771  cg_wolfe = .false.
772 
773  END FUNCTION cg_wolfe
774 
775 
776 
777  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
778  ! check for convergence of the cg iterations
779  ! (double) f -- function value associated with stepsize
780  ! (double) gnorm -- gradient (infinity) norm
781  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
782 
783  LOGICAL FUNCTION cg_tol (f, gnorm)
784 
785  ! use define_cg_variables
786  implicit none
787 
788  real, INTENT(IN OUT) :: f
789  real, INTENT(IN) :: gnorm
790  !---------------------------------------------------------------------------
791 
792  IF ( stoprule ) THEN
793  IF ( gnorm <= tol ) THEN
794  cg_tol = .true.
795  RETURN
796  END IF
797  ELSE
798  IF ( gnorm <= tol*(1.0 + abs(f)) ) THEN
799  cg_tol = .true.
800  RETURN
801  END IF
802  END IF
803  cg_tol = .false.
804 
805  END FUNCTION cg_tol
806 
807 
808 
809  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
810  ! compute dot product of x and y, vectors of length n
811  ! (double) x -- first vector
812  ! (double) y -- second vector
813  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
814 
815  real FUNCTION cg_dot (x, y)
816 
817  ! use define_cg_variables
818  implicit none
819 
820  real, INTENT(IN) :: x(:)
821  real, INTENT(IN) :: y(:)
822  real :: t
823 
824  integer :: i
825  !---------------------------------------------------------------------------
826 
827  t = zero
828  DO i = 1, n5
829  t = t + x(i)*y(i)
830  END DO
831  DO i = n6, n, 5
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)
834  END DO
835  cg_dot = t
836 
837  END FUNCTION cg_dot
838 
839 
840  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
841  ! compute xtemp = x + alpha d
842 
843  ! (double) xtemp -- output vector
844  ! (double) x -- initial vector
845  ! (double) d -- search direction vector
846  ! (double) alpha -- stepsize along search direction vector
847  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
848 
849  SUBROUTINE cg_step (xtemp, x, d, alpha)
850 
851  ! use define_cg_variables
852  implicit none
853 
854  real, INTENT(OUT) :: xtemp(:)
855  real, INTENT(IN) :: x(:)
856  real, INTENT(IN) :: d(:)
857  real, INTENT(IN) :: alpha
858 
859  integer :: i, j
860  !---------------------------------------------------------------------------
861 
862  DO i = 1, n5
863  xtemp(i) = x(i) + alpha*d(i)
864  END DO
865  DO i = n6, n, 5
866  xtemp(i) = x(i) + alpha*d(i)
867  j = i + 1
868  xtemp(j) = x(j) + alpha*d(j)
869  j = i + 2
870  xtemp(j) = x(j) + alpha*d(j)
871  j = i + 3
872  xtemp(j) = x(j) + alpha*d(j)
873  j = i + 4
874  xtemp(j) = x(j) + alpha*d(j)
875  END DO
876  END SUBROUTINE cg_step
877 
878 
879 
880  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
881  ! (double) alpha -- stepsize along search direction vector
882  ! (double) phi -- function value for step alpha
883  ! (double) dphi -- function derivative for step alpha
884  ! (double) dphi0 -- function derivative at starting point (alpha = 0)
885  ! (double) x -- current iterate
886  ! (double) xtemp -- x + alpha*d
887  ! (double) d -- current search direction
888  ! (double) gtemp -- gradient at x + alpha*d
889  ! (external) cg_value -- routine to evaluate function value
890  ! (external) cg_grad -- routine to evaluate function gradient
891  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
892 
893  SUBROUTINE cg_line (alpha, phi, dphi, dphi0, x, xtemp, d, gtemp, &
894  cg_value, cg_grad)
895 
896  ! use define_cg_variables
897  implicit none
898 
899 
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(:)
908 
909  INTERFACE
910 
911  SUBROUTINE cg_value (f, x, n)
912  IMPLICIT NONE
913  INTEGER, INTENT(IN) :: n
914  REAL, INTENT(IN) :: x(:)
915  REAL, INTENT(IN OUT) :: f
916  END SUBROUTINE cg_value
917 
918  SUBROUTINE cg_grad (g, x, n)
919  IMPLICIT NONE
920  INTEGER, INTENT(IN) :: n
921  REAL, INTENT(IN) :: x(:)
922  REAL, INTENT(IN OUT) :: g(:)
923  END SUBROUTINE cg_grad
924 
925  END INTERFACE
926 
927  real :: a, dphia, b, dphib, c
928  real :: a0, da0, b0, db0
929  real :: width, fquad
930 
931  integer :: ngrow, nshrink, iter, flag
932  !---------------------------------------------------------------------------
933 
934  CALL cg_step (xtemp, x, d, alpha)
935  CALL cg_grad (gtemp, xtemp, n)
936  ng = ng + 1
937  dphi = cg_dot(gtemp, d)
938 
939  ! Find initial interval [a,b] such that dphia < 0, dphib >= 0,
940  ! and phia <= phi0 + feps*dabs (phi0)
941 
942  a = zero
943  dphia = dphi0
944  ngrow = 0
945  nshrink = 0
946  DO WHILE ( dphi < zero )
947  CALL cg_value (phi, xtemp, n)
948  nf = nf + 1
949 
950  ! if quadstep in effect and quadratic conditions hold, check wolfe condition
951 
952  IF ( quadok ) THEN
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, &
958  ' fquad:', e14.6)
959  END IF
960  IF ( cg_wolfe(alpha, phi, dphi) ) RETURN
961  END IF
962  END IF
963  IF ( phi <= fpert ) THEN
964  a = alpha
965  dphia = dphi
966  ELSE
967 
968  ! contraction phase
969 
970  b = alpha
971  DO WHILE ( .true. )
972  alpha = .5d0*(a+b)
973  nshrink = nshrink + 1
974  IF ( nshrink > nexpand ) THEN
975  info = 6
976  RETURN
977  END IF
978  CALL cg_step(xtemp, x, d, alpha)
979  CALL cg_grad(gtemp, xtemp, n)
980  ng = ng + 1
981  dphi = cg_dot(gtemp, d)
982  IF ( dphi >= zero ) GO TO 100
983  CALL cg_value(phi, xtemp, n)
984  nf = nf + 1
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)
989  END IF
990  IF ( quadok .AND. (phi <= fquad) ) THEN
991  IF ( cg_wolfe(alpha, phi, dphi) ) RETURN
992  END IF
993  IF ( phi <= fpert ) THEN
994  a = alpha
995  dphia = dphi
996  ELSE
997  b = alpha
998  END IF
999  END DO
1000  END IF
1001 
1002  ! expansion phase
1003 
1004  ngrow = ngrow + 1
1005  IF ( ngrow > nexpand ) THEN
1006  info = 3
1007  RETURN
1008  END IF
1009  alpha = rho*alpha
1010  CALL cg_step(xtemp, x, d, alpha)
1011  CALL cg_grad(gtemp, xtemp, n)
1012  ng = ng + 1
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)
1018  END IF
1019  END DO
1020 100 CONTINUE
1021  b = alpha
1022  dphib = dphi
1023  IF ( quadok ) THEN
1024  CALL cg_value(phi, xtemp, n)
1025  nf = nf + 1
1026  IF ( ngrow + nshrink == 0 ) fquad = min(phi, f0)
1027  IF ( phi <= fquad ) THEN
1028  IF ( cg_wolfe(alpha, phi, dphi) ) RETURN
1029  END IF
1030  END IF
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)
1036  END IF
1037  width = gamma*(b - a)
1038  IF ( -dphia <= dphib ) THEN
1039  alpha = a - (a-b)*(dphia/(dphia-dphib))
1040  ELSE
1041  alpha = b - (a-b)*(dphib/(dphia-dphib))
1042  END IF
1043  c = alpha
1044  a0 = a
1045  b0 = b
1046  da0 = dphia
1047  db0 = 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
1051  RETURN
1052  ELSE IF ( flag == 0 ) THEN
1053  IF ( c == a ) THEN
1054  IF ( dphi > da0 ) THEN
1055  alpha = c - (c-a0)*(dphi/(dphi-da0))
1056  ELSE
1057  alpha = a
1058  END IF
1059  ELSE
1060  IF ( dphi < db0 ) THEN
1061  alpha = c - (c-b0)*(dphi/(dphi-db0))
1062  ELSE
1063  alpha = b
1064  END IF
1065  END IF
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
1071  END IF
1072  END IF
1073 
1074  ! bisection iteration
1075 
1076  IF ( (b-a) >= width ) THEN
1077  alpha = .5d0*(b+a)
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
1082  ELSE
1083  IF ( b <= a ) THEN
1084  info = 7
1085  RETURN
1086  END IF
1087  END IF
1088  END DO
1089  info = 4
1090 
1091  END SUBROUTINE cg_line
1092 
1093 
1094  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1095  ! This routine is identical to cg_line except that the function
1096  ! psi (a) = phi (a) - phi (0) - a*delta*dphi(0) is miniminized instead of
1097  ! the function phi
1098 
1099  ! (double) alpha -- stepsize along search direction vector
1100  ! (double) phi -- function value for step alpha
1101  ! (double) dphi -- function derivative for step alpha
1102  ! (double) dphi0 -- function derivative at starting point (alpha = 0)
1103  ! (double) x -- current iterate
1104  ! (double) xtemp -- x + alpha*d
1105  ! (double) d -- current search direction
1106  ! (double) gtemp -- gradient at x + alpha*d
1107  ! (external) cg_value -- routine to evaluate function value
1108  ! (external) cg_grad -- routine to evaluate function gradient
1109  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1110 
1111  SUBROUTINE cg_linew (alpha, phi, dphi, dphi0, x, xtemp, d, gtemp, cg_value, cg_grad)
1112 
1113  ! use define_cg_variables
1114  implicit none
1115 
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(:)
1124 
1125  INTERFACE
1126 
1127  SUBROUTINE cg_value (f, x, n)
1128  IMPLICIT NONE
1129  INTEGER, INTENT(IN) :: n
1130  REAL, INTENT(IN) :: x(:)
1131  REAL, INTENT(IN OUT) :: f
1132  END SUBROUTINE cg_value
1133 
1134  SUBROUTINE cg_grad (g, x, n)
1135  IMPLICIT NONE
1136  INTEGER, INTENT(IN) :: n
1137  REAL, INTENT(IN) :: x(:)
1138  REAL, INTENT(IN OUT) :: g(:)
1139  END SUBROUTINE cg_grad
1140 
1141  END INTERFACE
1142 
1143  real :: a, dpsia, b, dpsib, c, &
1144  a0, da0, b0, db0, width, fquad, psi, dpsi
1145 
1146  INTEGER :: ngrow, nshrink, iter, flag
1147  !---------------------------------------------------------------------------
1148 
1149 
1150  CALL cg_step (xtemp, x, d, alpha)
1151  CALL cg_grad (gtemp, xtemp, n)
1152  ng = ng + 1
1153  dphi = cg_dot(gtemp, d)
1154  dpsi = dphi - wolfe_hi
1155 
1156  ! Find initial interval [a,b] such that dpsia < 0, dpsib >= 0,
1157  ! and psia <= phi0 + feps*dabs(phi0)
1158 
1159  a = zero
1160  dpsia = dphi0 - wolfe_hi
1161  ngrow = 0
1162  nshrink = 0
1163  DO WHILE ( dpsi < zero )
1164  CALL cg_value(phi, xtemp, n)
1165  psi = phi - alpha*wolfe_hi
1166 
1167  nf = nf + 1
1168 
1169  ! if quadstep in effect and quadratic conditions hold, check wolfe condition
1170 
1171  IF ( quadok ) THEN
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, &
1177  ' fquad:', e14.6)
1178  END IF
1179  IF ( cg_wolfe(alpha, phi, dphi) ) RETURN
1180  END IF
1181  END IF
1182  IF ( psi <= fpert ) THEN
1183  a = alpha
1184  dpsia = dpsi
1185  ELSE
1186 
1187  ! contraction phase
1188 
1189  b = alpha
1190  DO WHILE ( .true. )
1191  alpha = .5d0*(a+b)
1192  nshrink = nshrink + 1
1193  IF ( nshrink > nexpand ) THEN
1194  info = 6
1195  RETURN
1196  END IF
1197  CALL cg_step (xtemp, x, d, alpha)
1198  CALL cg_grad (gtemp, xtemp, n)
1199  ng = ng + 1
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
1205  nf = nf + 1
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)
1210  END IF
1211  IF ( quadok .AND. (phi <= fquad) ) THEN
1212  IF ( cg_wolfe(alpha, phi, dphi) ) RETURN
1213  END IF
1214  IF ( psi <= fpert ) THEN
1215  a = alpha
1216  dpsia = dpsi
1217  ELSE
1218  b = alpha
1219  END IF
1220  END DO
1221  END IF
1222 
1223  ! expansion phase
1224 
1225  ngrow = ngrow + 1
1226  IF ( ngrow > nexpand ) THEN
1227  info = 3
1228  RETURN
1229  END IF
1230  alpha = rho*alpha
1231  CALL cg_step (xtemp, x, d, alpha)
1232  CALL cg_grad (gtemp, xtemp, n)
1233  ng = ng + 1
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
1241  END IF
1242  END DO
1243 100 CONTINUE
1244  b = alpha
1245  dpsib = dpsi
1246  IF ( quadok ) THEN
1247  CALL cg_value(phi, xtemp, n)
1248  nf = nf + 1
1249  IF ( ngrow + nshrink == 0 ) fquad = min(phi, f0)
1250  IF ( phi <= fquad ) THEN
1251  IF ( cg_wolfe(alpha, phi, dphi) ) RETURN
1252  END IF
1253  END IF
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)
1259  END IF
1260  width = gamma*(b - a)
1261  IF ( -dpsia <= dpsib ) THEN
1262  alpha = a - (a-b)*(dpsia/(dpsia-dpsib))
1263  ELSE
1264  alpha = b - (a-b)*(dpsib/(dpsia-dpsib))
1265  END IF
1266  c = alpha
1267  a0 = a
1268  b0 = b
1269  da0 = dpsia
1270  db0 = 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
1274  RETURN
1275  ELSE IF ( flag == 0 ) THEN
1276  IF ( c == a ) THEN
1277  IF ( dpsi > da0 ) THEN
1278  alpha = c - (c-a0)*(dpsi/(dpsi-da0))
1279  ELSE
1280  alpha = a
1281  END IF
1282  ELSE
1283  IF ( dpsi < db0 ) THEN
1284  alpha = c -(c-b0)*(dpsi/(dpsi-db0))
1285  ELSE
1286  alpha = b
1287  END IF
1288  END IF
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
1294  END IF
1295  END IF
1296 
1297  ! bisection iteration
1298 
1299  IF ( (b-a) >= width ) THEN
1300  alpha = .5d0*(b+a)
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
1305  ELSE
1306  IF ( b <= a ) THEN
1307  info = 7
1308  RETURN
1309  END IF
1310  END IF
1311  END DO
1312  info = 4
1313 
1314  END SUBROUTINE cg_linew
1315 
1316 
1317  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1318  ! update returns 1 if Wolfe condition is satisfied or too many iterations
1319  ! returns 0 if the interval updated successfully
1320  ! returns -1 if search done
1321 
1322  ! (double) a -- left side of bracketting interval
1323  ! (double) dphia -- derivative at a
1324  ! (double) b -- right side of bracketting interval
1325  ! (double) dphib -- derivative at b
1326  ! (double) alpha -- trial step (between a and b)
1327  ! (double) phi -- function value at alpha (returned)
1328  ! (double) dphi -- function derivative at alpha (returned)
1329  ! (double) x -- current iterate
1330  ! (double) xtemp -- x + alpha*d
1331  ! (double) d -- current search direction
1332  ! (double) gtemp -- gradient at x + alpha*d
1333  ! (external) cg_value -- routine to evaluate function value
1334  ! (external) cg_grad -- routine to evaluate function gradient
1335  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1336 
1337  INTEGER FUNCTION cg_update (a, dphia, b, dphib, alpha, phi, &
1338  dphi, x, xtemp, d, gtemp, cg_value, cg_grad)
1339 
1340  ! use define_cg_variables
1341  implicit none
1342 
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(:)
1354 
1355  INTERFACE
1356 
1357  SUBROUTINE cg_value (f, x, n)
1358  IMPLICIT NONE
1359  INTEGER, INTENT(IN) :: n
1360  REAL, INTENT(IN) :: x(:)
1361  REAL, INTENT(IN OUT) :: f
1362  END SUBROUTINE cg_value
1363 
1364  SUBROUTINE cg_grad (g, x, n)
1365  IMPLICIT NONE
1366  INTEGER, INTENT(IN) :: n
1367  REAL, INTENT(IN) :: x(:)
1368  REAL, INTENT(IN OUT) :: g(:)
1369  END SUBROUTINE cg_grad
1370 
1371  END INTERFACE
1372 
1373 
1374  INTEGER :: nshrink
1375  !---------------------------------------------------------------------------
1376 
1377  CALL cg_step (xtemp, x, d, alpha)
1378  CALL cg_value (phi, xtemp, n)
1379  nf = nf + 1
1380  CALL cg_grad (gtemp, xtemp, n)
1381  ng = ng + 1
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)
1386  END IF
1387  cg_update = 0
1388  IF ( cg_wolfe(alpha, phi, dphi) ) THEN
1389  cg_update = 1
1390  GO TO 110
1391  END IF
1392  IF ( dphi >= zero ) THEN
1393  b = alpha
1394  dphib = dphi
1395  GO TO 110
1396  ELSE
1397  IF ( phi <= fpert ) THEN
1398  a = alpha
1399  dphia = dphi
1400  GO TO 110
1401  END IF
1402  END IF
1403  nshrink = 0
1404  b = alpha
1405  DO WHILE ( .true. )
1406  alpha = .5d0*(a+b)
1407  nshrink = nshrink + 1
1408  IF ( nshrink > nexpand ) THEN
1409  info = 8
1410  cg_update = 1
1411  GO TO 110
1412  END IF
1413  CALL cg_step (xtemp, x, d, alpha)
1414  CALL cg_grad (gtemp, xtemp, n)
1415  ng = ng + 1
1416  dphi = cg_dot(gtemp, d)
1417  CALL cg_value (phi, xtemp, n)
1418  nf = nf + 1
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)
1423  END IF
1424  IF ( cg_wolfe(alpha, phi, dphi) ) THEN
1425  cg_update = 1
1426  GO TO 110
1427  END IF
1428  IF ( dphi >= zero ) THEN
1429  b = alpha
1430  dphib = dphi
1431  GO TO 100
1432  END IF
1433  IF ( phi <= fpert ) THEN
1434  IF ( printlevel ) THEN
1435  WRITE (6, *) "update a:", alpha, "dphia:", dphi
1436  END IF
1437  a = alpha
1438  dphia = dphi
1439  ELSE
1440  b = alpha
1441  END IF
1442  END DO
1443 100 CONTINUE
1444  cg_update = -1
1445 110 CONTINUE
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)
1450  END IF
1451 
1452  END FUNCTION cg_update
1453 
1454 
1455  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1456  ! This routine is identical to cg_update except that the function
1457  ! psi (a) = phi(a) - phi(0) - a*delta*dphi(0) is miniminized instead of
1458  ! the function phi
1459 
1460  ! update returns 1 if Wolfe condition is satisfied or too many iterations
1461  ! returns 0 if the interval updated successfully
1462  ! returns -1 if search done
1463 
1464  ! (double) a -- left side of bracketting interval
1465  ! (double) dpsia -- derivative at a
1466  ! (double) b -- right side of bracketting interval
1467  ! (double) dpsib -- derivative at b
1468  ! (double) alpha -- trial step (between a and b)
1469  ! (double) phi -- function value at alpha(returned)
1470  ! (double) dphi -- derivative of phi at alpha(returned)
1471  ! (double) dpsi -- derivative of psi at alpha(returned)
1472  ! (double) x -- current iterate
1473  ! (double) xtemp -- x + alpha*d
1474  ! (double) d -- current search direction
1475  ! (double) gtemp -- gradient at x + alpha*d
1476  ! (external) cg_value -- routine to evaluate function value
1477  ! (external) cg_grad -- routine to evaluate function gradient
1478  !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1479 
1480  INTEGER FUNCTION cg_updatew (a, dpsia, b, dpsib, alpha, phi, dphi, &
1481  dpsi, x, xtemp, d, gtemp, cg_value, cg_grad)
1482 
1483  ! use define_cg_variables
1484  implicit none
1485 
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(:)
1498 
1499  INTERFACE
1500 
1501  SUBROUTINE cg_value (f, x, n)
1502  IMPLICIT NONE
1503  INTEGER, INTENT(IN) :: n
1504  REAL, INTENT(IN) :: x(:)
1505  REAL, INTENT(IN OUT) :: f
1506  END SUBROUTINE cg_value
1507 
1508  SUBROUTINE cg_grad (g, x, n)
1509  IMPLICIT NONE
1510  INTEGER, INTENT(IN) :: n
1511  REAL, INTENT(IN) :: x(:)
1512  REAL, INTENT(IN OUT) :: g(:)
1513  END SUBROUTINE cg_grad
1514 
1515  END INTERFACE
1516 
1517  real :: psi
1518 
1519  INTEGER :: nshrink
1520  !---------------------------------------------------------------------------
1521 
1522  CALL cg_step (xtemp, x, d, alpha)
1523  CALL cg_value (phi, xtemp, n)
1524  psi = phi - alpha*wolfe_hi
1525  nf = nf + 1
1526  CALL cg_grad (gtemp, xtemp, n)
1527  ng = ng + 1
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)
1533  END IF
1534  cg_updatew = 0
1535  IF ( cg_wolfe(alpha, phi, dphi) ) THEN
1536  cg_updatew = 1
1537  GO TO 110
1538  END IF
1539  IF ( dpsi >= zero ) THEN
1540  b = alpha
1541  dpsib = dpsi
1542  GO TO 110
1543  ELSE
1544  IF ( psi <= fpert ) THEN
1545  a = alpha
1546  dpsia = dpsi
1547  GO TO 110
1548  END IF
1549  END IF
1550  nshrink = 0
1551  b = alpha
1552  DO WHILE ( .true. )
1553  alpha = .5d0*(a+b)
1554  nshrink = nshrink + 1
1555  IF ( nshrink > nexpand ) THEN
1556  info = 8
1557  cg_updatew = 1
1558  GO TO 110
1559  END IF
1560  CALL cg_step (xtemp, x, d, alpha)
1561  CALL cg_grad (gtemp, xtemp, n)
1562  ng = ng + 1
1563  dphi = cg_dot(gtemp, d)
1564  dpsi = dphi - wolfe_hi
1565  CALL cg_value (phi, xtemp, n)
1566  psi = phi - alpha*wolfe_hi
1567  nf = nf + 1
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)
1572  END IF
1573  IF ( cg_wolfe(alpha, phi, dphi) ) THEN
1574  cg_updatew = 1
1575  GO TO 110
1576  END IF
1577  IF ( dpsi >= zero ) THEN
1578  b = alpha
1579  dpsib = dpsi
1580  GO TO 100
1581  END IF
1582  IF ( psi <= fpert ) THEN
1583  IF ( printlevel ) THEN
1584  WRITE (6, *) "update a:", alpha, "dpsia:", dpsi
1585  END IF
1586  a = alpha
1587  dpsia = dpsi
1588  ELSE
1589  b = alpha
1590  END IF
1591  END DO
1592 100 CONTINUE
1593  cg_updatew = -1
1594 110 CONTINUE
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)
1599  END IF
1600 
1601  END FUNCTION cg_updatew
1602  ! Version 1.2 Changes:
1603 
1604  ! 1. Fix problem with user specified initial step (overwriting step)
1605  ! 2. Change dphi to dpsi at lines 1228 and 1234 in cg_lineW
1606  ! 3. Add comment about how to compute dnorm2 by an update of previous dnorm2
1607  ! 4. In comment statements for cg_lineW and cg_updateW, insert "delta"
1608  ! in definition of psi (a)
1609  ! 5. In dimension statements, change "(1)" to "(*)"
1610 
1611  ! Version 1.3 Changes:
1612  ! 1. Remove extraneous write in line 985 (same thing written out twice)
1613  ! 2. Remove the parameter theta from cg_descent.parm and from the code
1614  ! (we use theta = .5 in the cg_update)
1615 
1616  ! Version 1.4 Change:
1617  ! 1. The variable dpsi needs to be included in the argument list for
1618  ! subroutine updateW (update of a Wolfe line search)
1619 
1620 END MODULE cg_minimization