Next: About this document ... Up: Fomel: Regularization Previous: REFERENCES

# 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
```

Next: About this document ... Up: Fomel: Regularization Previous: REFERENCES
Stanford Exploration Project
11/11/1997