module weight_solver { logical, parameter, private :: T = .true., F = .false. contains subroutine solver( oper, solv, x, dat, niter, nmem, nfreq, wght) { interface { integer function wght( res, w) { real, dimension (:) :: res, w } integer function oper( adj, add, x, dat) { logical, intent (in) :: adj, add real, dimension (:) :: x, dat } integer function solv( forget, x, g, rr, gg) { logical :: forget real, dimension (:) :: x, g, rr, gg } } real, dimension (:), intent (in) :: dat # data real, dimension (:), intent (out) :: x # solution integer, intent (in) :: niter, nmem, nfreq # iterations real, dimension (size (x)) :: g # gradient real, dimension (size (dat)) :: rr, gg, wt # res, CG, weight integer :: i, stat logical :: forget rr = -dat; x = 0.; wt = 1.; forget = F # initialize do i = 1, niter { forget = (i > nmem) # restart if( forget) stat = wght( rr, wt) # compute weighting rr = rr * wt # rr = W (Fx - d) stat = oper( T, F, g, wt*rr) # g = F' W' rr stat = oper( F, F, g, gg) # gg = Fg gg = gg * wt # gg = W F g if( forget) forget = ( mod( i, nfreq) == 0) # periodic restart stat = solv( forget, x, g, rr, gg) # step in x and rr rr = - dat stat = oper( F, T, x, rr) # rr = Fx - d } } }