Regularized Reusable Linear Solver

Experimenting with different forms of regularization, different kinds of interpolation and regularization operators, and different optimization algorithms was made easy with an improved design of the generic solver subroutine (module.)

The Fortran-90 subprograms in this appendix are simplified from the actual version, which is more heavily loaded with optional parameters. For example, the actual solver routine can implement an iterative reweighting that I used in the deburst example, work with non-linear operators, etc. Two solver subroutines together with some others (e.g. LSQR solver) are packed in a MODULE, with overloading of the subroutine names, so that the user can easily switch from one solver to another without changing the name of the program (and without even knowing the details of its implementation.)

The data-space regularization routine takes three functions as its arguments. Functions oper and reg correspond to the linear operators L and P. They comply with the generic linear operator interface, defined by Fomel and Claerbout (1996). Function solv implements one step of an optimization descent. Its arguments are a logical parameter forget, which controls a conditional restarting of the optimization, the current effective model mod, the gradient vector g, the data residual vector rr, and the conjugate gradient vector gg. An example of a solv function is conjgrad , which implements a classic version of the conjugate-gradient method. Subroutine lin_solver_dat constructs the effective model vector x, which consists of the model-space part xm and the data-space part xd. Similarly, the effective gradient vector g is split into the the model-space part gm and the data-space part gd. Subroutine chain from module chainmod is called to compute a chain of the operators P and L (reg and oper).

subroutine lin_solver_dat (oper, solv, reg, nreg, mod, dat, niter, eps, mod0) 
  interface
     integer function oper (adj, add, mod, dat)           
       logical, intent (in) :: adj, add
       real, dimension (:)  :: mod, dat
     end function oper
     integer function solv (forget, mod, g, rr, gg)       
       logical              :: forget
       real, dimension (:)  :: mod, g, rr, gg
     end function solv
     integer function reg (adj, add, mod, dat)            
       logical, intent (in) :: adj, add
       real, dimension (:)  :: mod, dat
     end function reg
  end interface
  real, dimension (:),   intent (in)           :: dat         ! data 
  real, dimension (:),   intent (in), optional :: mod0        ! initial model
  real, dimension (:),   intent (out)          :: mod         ! model 
  integer,               intent (in)           :: niter, nreg ! size of x
  real,                  intent (in)           :: eps         ! scaling

real, dimension (size (dat) + nreg), target  :: x, g
  real, dimension (size (dat))                 :: rr, gg
  real, dimension (:), pointer                 :: xm, xd, gm, gd
  integer                                      :: iter, stat
  logical                                      :: forget

xm => x (1 : nreg) ; xd => x (1 + nreg:) ; xd = 0.
  gm => g (1 : nreg) ; gd => g (1 + nreg:)
  if (present (mod0)) then
     xm = mod0 ; call chain (oper, reg, .false., .false., xm, rr, mod) 
     rr = dat - rr
  else
     xm = 0. ; rr = dat
  end if

forget = .false.
  do iter = 1, niter
     call chain (oper, reg, .true. , .false., gm, rr, mod)  ; gd =   eps*rr
     call chain (oper, reg, .false., .false., gm, gg, mod)  ; gg =gg+eps*gd
     stat = solv (forget, x, g, rr, gg)
  end do
  stat = reg (.false., .false., xm, mod)
end subroutine lin_solver_dat

The model-space regularization routine has an analogous design. In this case, function reg corresponds to the model regularization operator D. We construct the effective residual vector rr, which consists of the model-space part rm and the data-space part rd. Similarly, the effective conjugate gradient vector gg is split into the the model-space part gm and the data-space part gd. Subroutine array from module chainmod is called to compute an array of the operators L and D (oper and reg).

subroutine lin_solver_mod (oper, solv, reg, nreg, mod, dat, niter, mod0)
  interface
     integer function oper (adj, add, mod, dat)
       logical, intent (in) :: adj, add
       real, dimension (:)  :: mod, dat
     end function oper
     integer function solv (forget, mod, g, rr, gg)
       logical             ::forget
       real, dimension (:) :: mod, g, rr, gg
     end function solv
     integer function reg (adj, add, mod, dat)
       logical, intent (in) :: adj, add
       real, dimension (:)  :: mod, dat
     end function reg
  end interface
  real, dimension (:), intent (in)            :: dat         ! data
  real, dimension (:), intent (in), optional  :: mod0        ! initial model  
  real, dimension (:), intent (out)           :: mod         ! model
  integer,             intent (in)            :: niter, nreg ! size of D mod

real, dimension (size (mod))                :: g
  real, dimension (size (dat) + nreg), target :: rr, gg
  real, dimension (:), pointer                :: rd, rm, gd, gm
  integer                                     :: iter, stat
  logical                                     :: forget

rm => rr (1 : nreg) ; rd => rr (1 + nreg :)
  gm => gg (1 : nreg) ; gd => gg (1 + nreg :)

if (present (mod0)) then
     mod = mod0 
     stat = oper (.false., .false., mod, rr) ; rr = dat - rr
     stat = reg  (.false., .false., mod, rm) ; rm =     - rm
  else
     mod = 0. ; rd = dat ; rm = 0.
  end if

forget = .false.
  do iter = 1, niter
     call array (oper, reg, .true. , .false., g, rd, rm) 
     call array (oper, reg, .false., .false., g, gd, gm)
     stat = solv (forget, mod, g, rr, gg)
  end do
end subroutine lin_solver_mod

module chainmod
contains

subroutine chain (oper1, oper2, adj, add, mod, dat, tmp)
    interface
       integer function oper1 (adj, add, mod, dat)
         logical, intent (in)  :: adj, add
         real, dimension (:)   :: mod, dat
       end function oper1
       integer function oper2 (adj, add, mod, dat)
         logical, intent (in)  :: adj, add
         real, dimension (:)   :: mod, dat
       end function oper2
    end interface
    logical, intent (in)  :: adj, add
    real, dimension (:)   :: mod, dat, tmp

integer               :: stat1, stat2
    if (adj) then
       stat1 = oper1 ( .true., .false., tmp, dat)
       stat2 = oper2 ( .true.,     add, mod, tmp)
    else
       stat2 = oper2 (.false., .false., mod, tmp)
       stat1 = oper1 (.false.,     add, tmp, dat)
    end if
  end subroutine chain

subroutine array (oper1, oper2, adj, add, mod, dat1, dat2)
    interface
       integer function oper1 (adj, add, mod, dat)
         logical, intent (in) :: adj, add
         real, dimension (:)  :: mod, dat
       end function oper1
       integer function oper2 (adj, add, mod, dat)
         logical, intent (in) :: adj, add
         real, dimension (:)  :: mod, dat
       end function oper2
    end interface
    logical, intent (in)  :: adj, add
    real, dimension (:)   :: mod, dat1, dat2

integer               :: stat1, stat2
    if (adj) then
       stat1 = oper1 (.true.,    add, mod, dat1)
       stat2 = oper2 (.true., .true., mod, dat2)
    else
       stat1 = oper1 (.false.,   add, mod, dat1)
       stat2 = oper2 (.false.,   add, mod, dat2)
    end if
  end subroutine array

end module chainmod

module conjgrad_mod  
  real, dimension (:), allocatable, private  :: s, ss
  real, parameter,                  private  :: eps = 1.e-12

contains

subroutine conjgrad_close ()
    deallocate (s, ss)
  end subroutine conjgrad_close

function conjgrad (forget, x, g, rr, gg) result (stat)
    integer             :: stat
    real, dimension (:) :: x, g, rr, gg
    logical             :: forget

real, save          :: rnp
    real                :: rn, alpha, beta

rn = dot_product (g, g)
    if (.not. allocated (s)) then
       allocate (s  (size (x ))) 
       allocate (ss (size (rr)))
       forget = .true.
    end if

if (forget .or. rnp < eps) then
       alpha = 0.d0
    else
       alpha = rn / rnp
    end if

s =   g + alpha *  s    ! model step
    ss = gg + alpha * ss    ! data  step

beta = dot_product (ss, ss)
    if (beta > eps) then
       alpha = rn / beta
       x  =  x + alpha *  s ! update model
       rr = rr - alpha * ss ! update residual
    end if

rnp = rn ; forget = .false. 
    stat = 0
  end function conjgrad

end module conjgrad_mod