A basic 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 basic solver is solver_tiny(). solver_tinytiny solver

The two most important arguments in solver_tiny() are the operator function Fop, which is defined by the interface from Chapter , and the stepper function stepper, 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_tiny ( matmult_lop, cgstep, ...

so while you are reading the solver_tiny module, you should visualize the Fop() function as being matmult_lop, and you should visualize the stepper() function as being cgstep.

The other required parameters to solver_tiny() are d (the data we want to fit), m (the model we want to estimate), and niter (the maximum number of iterations). There are also a couple of optional arguments. For example, m0 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 resd, 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.