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 ! scalingreal, 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 modreal, 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 containssubroutine 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-12contains
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