PRECONDITIONING AND INTERVAL VELOCITY

The basic idea of preconditioning as described above is that we start from fitting goals   

$$\begin{array} {lll} \bold 0 &\approx& \bold F \bold m - \bold d \\ \bold 0 &\approx& \bold A \bold m\end{array}$$ (8)

and change variables with $\bold m = \bold S \bold p$ getting

$$\begin{array} {lll} \bold 0 &\approx & \bold F \bold S \bo... ...- \bold d \\ \bold 0 &\approx & \epsilon\ \bold p \end{array}$$ (9)

where we have established $\bold A\bold S=\bold I$(to take advantage of the second CG miracle), either by finding $\bold S = \bold A^{-1}$or by revising the problem formulation $\bold A$ to be $\bold S^{-1}$.

It would seem that the next step is to invoke a regularized solver such as module reg_solver with the new operator $\bold F \bold S$and damping operator $\bold I$.Alternately, a new reusable preconditioned solver is the module prec_solver .  

module prec_solver {
  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 solv (forget, x, g, rr, gg) {
         logical             :: forget
         real, dimension (:) :: x, g, rr, gg
       }
       integer function prec (adj, add, x, dat) {
         logical, intent (in) :: adj, add
         real, dimension (:)  :: x, dat
       }
    }
    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 (nprec)           :: g, p              # gradient, precond
    real, dimension (size (dat) + nprec), target :: rr, gg # residual, conj grad
    real, dimension (:), pointer      :: rd, rm, gd, gm
    integer                           :: i, stat1, stat2, stat
    rm => rr (1 : nprec) ; rd => rr (1 + nprec :)          # model and data resids
    gm => gg (1 : nprec) ; gd => gg (1 + nprec :)          # model and data grads
                           rd = -dat                       # initialize r_d
    if (present (x0)) {
       stat1 = prec( F, F, x0, x)       # x = S p
       stat2 = oper( F, T, x, rd)       # r_d = L x - dat 
                  p = x0 ; rm = x0*eps} # start with x0      
    else {        p = 0. ; rm = 0.}     # start with zero
    do i = 1, niter {
       stat2 = oper( T, F, x, rd)    
       stat1 = prec( T, F, g, x)        # g = S' L' r_d
       g = g + eps*rm                   # g = S' L' r_d + eps I r_m
       stat1 = prec( F, F, g, x)     
       stat2 = oper( F, T, x, gd)       # g_d = L S g
       gm = eps*g                       # g_m = eps I g
       stat = solv (F, p, g, rr, gg)    # step in p and rr
    }
    stat1 = prec( F, F, p, x)           # x = S p
  }
}