module vr_solver {			            # virtual-residual fitter
  logical, parameter, private  :: T = .true., F = .false.
contains
  subroutine solver_prec( oper, solv, prec, nprec, x, dat, niter, eps, x0) {
    optional                                                        :: x0
    interface {
       integer function oper( adj, add, x, dat) {
         logical, intent (in) :: adj, add
         real, dimension (:)  :: x, dat
       }
       integer function prec( 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, x0         # data, initial
    real, dimension (:), intent (out) :: x               # solution
    real,                intent (in)  :: eps             # scaling
    integer,             intent (in)  :: niter, nprec    # iterations, size
    real, dimension (size (dat) + nprec), target :: p, g # new, gradient
    real, dimension (size (dat))      :: rr, gg          # residual, conj grad
    real, dimension (:), pointer      :: pm, pd, gm, gd
    integer                           :: i, stat1, stat2, stat

    pm => p(1:nprec);   pd => p(1+nprec:);  pd = 0.      # compound model
    gm => g(1:nprec);   gd => g(1+nprec:)                # data and model grads
                                     rr = - dat
    if( present( x0)) {              pm = x0             # start with x0
       stat2 = prec( F, F, pm, x)                        # x0 = S p0 
       stat1 = oper( F, T, x, rr)                        # r  = F x0 - dat
    } else {                         pm = 0.}            # start with zero
    do i = 1, niter {                                    
       stat1 = oper( T, F, x, rr)                
       stat2 = prec( T, F, gm, x)                    # g_m = S' F' rr
                           gd =     eps*rr           # g_d = eps rr
       stat2 = prec( F, F, gm, x)                
       stat1 = oper( F, F,  x, gg)                   # gg = F S g_m
                               gg += eps*gd          # gg = F S g_m + eps g_d
       stat = solv( F, p, g, rr, gg)                 # step in p and rr
    }
    stat2 =    prec( F, F, pm, x)                    # x = S p 
  }
}