1      subroutine ga_pcg_minim(n, iter, e, g_grad, g_work, g_s, step,
2     $     ls_tol, ls_max, eacc, conjugacy, oprint, oprint_ls,
3     $     iwork, dwork, mode)
4C$Id$
5      implicit none
6#include "global.fh"
7#include "mafdecls.fh"
8#include "util.fh"
9      integer n                 ! No. of parameters [input]
10      integer iter              ! Current macro iteration no. [input/output]
11      double precision e        ! Energy of current point [input]
12      integer g_grad            ! GA with current gradient [input]
13      integer g_work            ! GA for work space [input]
14      integer g_s               ! GA with current search direction [output]
15      double precision step     ! Step to take along the line search direction [output]
16      double precision ls_tol   ! Tolerance for current line search [input]
17      integer ls_max            ! Max. no. of steps for line search [input]
18      double precision eacc     ! Absolution precision in energy [input]
19      logical conjugacy         ! True if conjugacy is employed [input]
20      logical oprint            ! True if should print convergence [input]
21      logical oprint_ls         ! True if should print line search [input]
22      integer iwork(1000)       ! Integer work space [input]
23      double precision dwork(3*ls_max) ! Double precision work space [input]
24      character*16 mode         ! Used to communicate with calling routine
25c
26c     Use preconditioned conjugate gradient to minimize the value of
27c     a function.  It is assumed that the preconditioning is pretty
28c     good so that steps are generally in the vicinity of unity
29c     (this assumption is from the line-search).
30c
31c     Have a look in testsolve.F (routine testpcg) for sample usage.
32c
33c     g_grad holds the current gradient
34c     g_s    holds the current search direction
35c     g_work holds the previous preconditioned gradient
36c
37      integer action, act_startup, act_ls, act_precond, act_conj,
38     $     act_accept
39      parameter (act_startup=0, act_precond=1, act_conj=2, act_ls=3,
40     $     act_accept=4)
41c
42c     Pointers into iwork and dwork for info
43c
44      integer iwork_action      ! Action
45      integer iwork_ls_nstep    ! Current no. of line search steps
46      parameter (iwork_action = 1, iwork_ls_nstep = 2)
47c
48      integer ls_e, ls_g, ls_s, sk_gk, sk_gk1, gk1_c_gk, gk1_c_gk1
49      parameter (sk_gk=1, sk_gk1=2, gk1_c_gk=3, gk1_c_gk1 = 4, ls_e = 5)
50c
51      logical ls_print          ! oprint_ls restricted to node 0
52      integer nstep
53      logical line_search
54      double precision beta, numerator, denominator, gnorm, gmax
55      double precision mone
56      parameter (mone = -1.0d0)
57      external line_search
58c
59      if (mode .eq. 'startup') iwork(iwork_action) = act_startup
60c
61c     Set up pointers and look inside iwork to see what we are
62c     currently doing
63c
64      action = iwork(iwork_action)
65      nstep  = iwork(iwork_ls_nstep)
66      ls_g = ls_e + ls_max
67      ls_s = ls_g + ls_max
68      ls_print = .false.
69c
70      if (action .eq. act_startup) then
71c
72c     Startup calculation ... get energy and gradient and the
73c     next thing needed is a preconditioning
74c
75         iter = 1
76         step = 0.0d0
77         call dfill((3*ls_max), 0.0d0, dwork, 1)
78         call ifill(2, 0, iwork, 1)
79         call ga_zero(g_work)
80         call ga_zero(g_s)
81         mode = 'energy+gradient'
82         action = act_accept
83      else if (action .eq. act_precond) then
84c
85c     Preconditioning required.  Next thing to do is the conjugacy
86c
87         dwork(sk_gk1) = ga_ddot(g_s, g_grad)
88         dwork(gk1_c_gk) = ga_ddot(g_grad, g_work)
89         mode = 'precondition'
90         action = act_conj
91      else if (action .eq. act_conj) then
92c
93c     Apply the conjugacy condition and then commence the line search
94c
95         dwork(gk1_c_gk1) = ga_ddot(g_grad, g_work)
96         if (.not. conjugacy) then
97            beta = 0.0d0
98         else
99            numerator = (dwork(gk1_c_gk1) - dwork(gk1_c_gk))
100            denominator = (dwork(sk_gk1)  - dwork(sk_gk))
101            if (denominator .ne. 0.0d0) then
102               beta = numerator/denominator
103            else
104               beta = 0.0d0
105            endif
106         endif
107         if (beta .lt. 0.0d0) then
108*            if (ga_nodeid() .eq. 0) then
109*               write(6,*) ' pcg: negative beta = ', beta
110*               call util_flush(6)
111*            endif
112            beta = 0.0d0
113         endif
114         call ga_dadd(mone, g_work, beta, g_s, g_s)
115         dwork(sk_gk) = ga_ddot(g_s, g_grad) ! For the next iteration
116c
117         step = 0.0d0
118         nstep = 1
119         mode = 'energy+gradient'
120         action = act_ls
121      else if (action .eq. act_accept) then
122c
123c     Accept the current point ... next will be a precond
124c
125         mode = 'accept step'
126         action = act_precond
127      endif
128c
129      if (action .eq. act_ls) then
130c
131c     We have a new energy+gradient for a line search underway.  Determine
132c     if the LS has converged or what the next step should be
133c
134         dwork(ls_e+nstep-1) = e
135         dwork(ls_g+nstep-1) = ga_ddot(g_grad, g_s)
136         dwork(ls_s+nstep-1) = step
137         ls_print = oprint_ls .and. (ga_nodeid().eq.0)
138         if (line_search(nstep, ls_max, dwork(ls_e), dwork(ls_g),
139     $        dwork(ls_s), eacc, ls_tol, ls_print)) then
140c
141c     The LS has converged.  Return the current step as accepted.
142c     Next action will be to precondition the gradient.
143c
144            mode = 'accept step'
145            action = act_precond
146         else
147c
148c     Just continue with the line search
149c
150            step = dwork(ls_s+nstep-1)
151         endif
152      endif
153c
154      if (mode .eq. 'accept step') then
155         if (oprint) then
156            gnorm = sqrt(ga_ddot(g_grad, g_grad))
157            call ga_maxelt(g_grad, gmax)
158            if (ga_nodeid() .eq. 0) then
159c
160c     Don't write out the header unless are also printing LS info
161c     except on the first iteration ... keeps output compact
162c
163               if (ls_print .or. iter.eq.1) write(6,1)
164 1             format(/,
165     $              13x,' iter       energy          gnorm     gmax   ',
166     $              '    time'/
167     $              13x,'----- ------------------- --------- ---------',
168     $              ' --------')
169               write(6,2) iter, e, gnorm, gmax, util_cpusec()
170 2             format(13x,i5,f20.10,1p,2d10.2,0p,f9.1)
171               call util_flush(6)
172            endif
173         endif
174         iter = iter + 1
175      endif
176c
177      iwork(iwork_action) = action
178      iwork(iwork_ls_nstep) = nstep
179c
180      end
181
182