Setting up a generic solver program

There are many different methods for iterative least-square estimation some of which will be discussed later in this book. The conjugate-gradient (CG) family (including the first order conjugate-direction method described above) share the property that theoretically they achieve the solution in n iterations, where n is the number of unknowns. The various CG methods differ in their numerical errors, memory required, adaptability to non-linear optimization, and their requirements on accuracy of the adjoint. What we do in this section is to show you the generic interface.

None of us is an expert in both geophysics and in optimization theory (OT), yet we need to handle both. We would like to have each group write its own code with a relatively easy interface. The problem is that the OT codes must invoke the physical operators yet the OT codes should not need to deal with all the data and parameters needed by the physical operators.

In other words, if a practitioner decides to swap one solver for another, the only thing needed is the name of the new solver.

The operator entrance is for the geophysicist, who formulates the estimation problem. The solver entrance is for the specialist in numerical algebra, who designs a new optimization method.

The Fortran-90 programming language allows us to achieve this design goal by means of generic function interfaces.

A generic solver subroutine solver() is shown in module simple_solver. It is simplified substantially from the library version, which has a much longer list of optional arguments  

module simple_solver  {
  logical, parameter, private  :: T = .true., F = .false.
contains
  subroutine solver( oper, solv, x, dat, niter, x0, res)  {
    optional                                 :: x0, res   
    interface {
       function oper( adj, add, x, dat)  result (stat) {
            integer              :: stat
            logical, intent (in) :: adj, add
            real, dimension (:)  :: x, dat
            }
       function solv( forget, x, g, rr, gg)  result (stat) {
            integer             :: stat
            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
          }
    do i = 1, niter  {
       stat = oper( T, F, g, rr)                          # g  <- F' rr
       stat = oper( F, F, g, gg)                          # gg <- F  g
       stat = solv( F, x, g, rr, gg)                      # step in x and rr
       }
    if( present( res)) res = rr
  }
}

(The forget parameter is not needed by the solvers we discuss first.)

The two most important arguments in solver() are the operator function oper, which is defined by the interface from Chapter , and the solver function solv, which implements one step of an iterative estimation. For example, a practitioner who choses to use our new cgstep() for iterative solving the operator matmult would write the call

call solver ( matmult_lop, cgstep, ...

The other required parameters to solver() are dat (the data we want to fit), x (the model we want to estimate), and niter (the maximum number of iterations). There is also a couple of optional arguments. For example, x0 is the starting guess for the model. If this parameter is omitted, the model is initialized to zero. To output the final residual vector, we include a parameter called res, which is optional as well. We will watch how the list of optional parameters to the generic solver routine grows as we attack more and more complex problems in later chapters.