Because Fortran does not recognize the difference between upper- and
lower-case letters,
the conjugate vectors G and S in the program are denoted by
gg and ss.
The inner part of the conjugate-direction task is in
function cgstep().
module cgstep_mod {
real, dimension (:), allocatable, private :: s, ss
contains
function cgstep( forget, x, g, rr, gg) result (stat) {
integer :: stat
real, dimension (:) :: x, g, rr, gg
logical :: forget
double precision :: sds, gdg, gds, determ, gdr, sdr, alfa, beta
if( .not. allocated (s)) { forget = .true.
allocate ( s (size ( x)))
allocate (ss (size (rr)))
}
if( forget){ s = 0.; ss = 0.; beta = 0.d0 # steepest descent
if( dot_product(gg, gg) == 0 )
call erexit('cgstep: grad vanishes identically')
alfa = - sum( dprod(gg,rr)) / sum( dprod(gg, gg))
}
else{ gdg = sum( dprod(gg, gg)) # search plane by solving 2-by-2
sds = sum( dprod(ss, ss)) # G . (R - G*alfa - S*beta) = 0
gds = sum( dprod(gg, ss)) # S . (R - G*alfa - S*beta) = 0
if( gdg==0. .or. sds==0.) { stat = 1; return }
determ = gdg * sds * max (1.d0 - (gds/gdg)*(gds/sds), 1.d-12)
gdr = - sum( dprod(gg,rr))
sdr = - sum( dprod(ss,rr))
alfa = ( sds * gdr - gds * sdr ) / determ
beta = (-gds * gdr + gdg * sdr ) / determ
}
s = alfa * g + beta * s # update solution step
ss = alfa * gg + beta * ss # update residual step
x = x + s # update solution
rr = rr + ss # update residual
forget = .false.; stat = 0
}
subroutine cgstep_close ( ) {
if (allocated (s)) deallocate (s, ss)
}
}