Hybrid-norm and Fortran 2003: Separating the physics from the solver

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 estimating problems. Iterative solvers can become quite complex, quite quickly, 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 viewed 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 recent introduction of 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 2003, which is more complete object-oriented language.

In this paper, I show how to implement an abstract solver class in Fortran 2003. I begin by describing one method to implement an abstract operator-based iterative solver. I describe how the definition of a vector needs to be extended to support the hybrid norm. I then show how the abstract operator and vector classes can be implemented in Fortran2003.

Hybrid-norm and Fortran 2003: Separating the physics from the solver