module conjgrad_mod implicit none real, dimension (:), allocatable, private :: s, ss contains subroutine conjgrad_close () if ( allocated( s)) then deallocate( s, ss) end if end subroutine 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)) then forget = .true. allocate( s (size (x ))) s = 0. allocate( ss (size (rr))) ss = 0. end if if ( forget .or. rnp < epsilon (rnp)) then alpha = 0.d0 else alpha = rn / rnp end if s = g + alpha * s ss = gg + alpha * ss beta = sum( dprod( ss, ss)) if ( beta > epsilon( beta)) then alpha = - rn / beta x = x + alpha * s rr = rr + alpha * ss end if rnp = rn forget = .false. stat = 0 end function end module