Krylov space solver in Fortran 2009: Beta version

Introduction

A geophysicist understands and/or approximates how a given set of earth properties (model) would create a given set of measurements (data). Geophysics is often an attempt to do the inverse: from a set of recorded data, estimate a model. When the set of measurements and/or the number of model points is large, matrix-based approaches become impractical. Iterative approaches are often the method of choice for large-scale estimation problems. Iterative solvers can become quite complex and are generally more the domain of the mathematician than the geophysicist. Ideally we would like to leverage the mathematician's expertise without having to understand all of the details of the implementation. Gockenbach (1994); Nichols et al. (1993) all implemented model estimation through an object-oriented framework, allowing this separation using C++. Schwab (1998) described a java-based approach to this problem, and Clapp (2005) described a python-based approach for large, out-of-core solvers. SEP chose instead to use Fortran 90. Unfortunately, Fortran 90 is far from a complete object-oriented language, and as a result complicated inversion problems are cumbersome to describe given its limitations. The problems encountered in implementing the hybrid norm (Zhang and Claerbout, 2010; Claerbout, 2009) is but one example of the limitations of Fortran 90 for solving inverse problems. Recently, Fortran compilers have begun to support Fortran's latest incarnation, Fortran 2008, a more complete implementation of the needed object oriented constructs. This paper is a follow up to Clapp (2010), which described how an implementation could be done in Fortran 2003, but was hampered by the immaturity of the Fortran compiler.

In this paper, I show the implementation of an abstract solver class in Fortran 2008. I begin by describing one method to implement an abstract operator-based iterative solver. I describe an abstract vector, abstract operator, and solver class. I finish by showing how a small geophysical inverse problems can be solved using these classes and concrete classes inherited from them.

Krylov space solver in Fortran 2009: Beta version