module gmres_mod use ddot_mod use llist_mod implicit none type(list), private :: grads real, dimension(:), allocatable, private :: s, ss contains subroutine gmres_init () grads = llist_init () end subroutine gmres_init subroutine gmres_close () call llist_close (grads) end subroutine gmres_close function gmres(forget, x, grad, rr, cgrad) result(stat) integer :: stat real, dimension(:):: x, grad, rr, cgrad logical :: forget real, dimension(:), pointer :: g, gg, gi double precision :: alpha, beta, gn double precision, parameter :: eps = 1.d-30 integer :: i, n, nm, nd nm = size (x) nd = size (rr) allocate(gg (nm + nd)) ; g => gg(: nm) g = grad ; gg(nm + 1:) = cgrad if (.not. allocated(s)) allocate(s(nm), ss(nd)) n = llist_depth (grads) call llist_rewind (grads) do i = 1, n call llist_down (grads, gi, beta) alpha = ddot (grad, gi (:nm)) / beta gg = gg - alpha * gi end do gn = ddot (g, g) if (n > 0) then alpha = gn / beta else alpha = 0.d0 end if s = g + alpha * s ss = gg (nm+1 :) + alpha * ss beta = ddot (ss, ss) if (beta > eps) then alpha = gn / beta x = x + alpha * s rr = rr - alpha * ss end if if (gn > eps) then call llist_add (grads, gg, gn) if (forget) call llist_chop (grads) end if stat = 0 end function gmres end module gmres_mod