module smallsolver  {
  logical, parameter, private  :: T = .true., F = .false.
  logical,            private  :: forget
contains
  subroutine solver( oper, solv, x, dat, niter, x0, res)  {
    optional                                 :: x0, res   
    interface {
       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, x0     # data, initial
    real,    dimension (:),   intent (out) :: x, res      # solution, residual
    integer,                  intent (in)  :: niter       # iterations
    real, dimension (size (x))             :: g           # gradient
    real, dimension (size (dat))           :: rr, gg      # residual, conj grad
    integer                                :: i, stat

    rr = - dat
    if( present( x0)) {  
	stat = oper( F, T, x0, rr)                        # rr <- F x0 - dat  
        x = x0                                            # start with x0
	}
    else {                              
	x = 0.                                            # start with zero
	}
    forget = F
    do i = 1, niter  {
       stat = oper( T, F, g, rr)                          # g  <- F' rr
       stat = oper( F, F, g, gg)                          # gg <- F  g
       stat = solv( forget, x, g, rr, gg)                      # step in x and rr
       }
    if( present( res)) res = rr
  }
}