module conjgrad_mod {
  real, dimension (:), allocatable, private 	:: s, ss
contains
  subroutine conjgrad_close ()  {
        if( allocated( s))  deallocate( s, ss)
        }
  function conjgrad( forget, x, g, rr, gg)  result( stat) {
    integer              :: stat
    real, dimension (:)  :: x, g, rr, gg
    logical              :: forget
    real, save           :: rnp
    double precision     :: rn, alpha, beta
    rn = sum( dprod( g, g))
    if( .not. allocated( s)) {   forget = .true.
           allocate( s  (size (x )));  s  = 0. 
           allocate( ss (size (rr)));  ss = 0.
           }
    if( forget .or. rnp < epsilon (rnp))
           alpha = 0.d0
    else
           alpha = rn / rnp
    s  =  g + alpha *  s
    ss = gg + alpha * ss
    beta = sum( dprod( ss, ss))
    if( beta > epsilon( beta)) {
           alpha = - rn / beta
           x =  x  + alpha *  s
           rr = rr + alpha * ss
           }
    rnp = rn;    forget = .false.;   stat = 0
    }
}