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