Introduction

Geophysical data analysis often calls for reconstruction of a model of the Earth's subsurface (the model, $\bold m$) from measurements of a physical quantity recorded at some distance (the data, $\bold D$). The relationship between the model and the data are often non-linear. However, in practice, we linearize this relationship, so that we can formulate and attempt to solve a problem in linear form.

Mathematically, we can represent our optimization problem through an equation like

$$\mathcal L\bold m\approx \bold D,$$

where $\mathcal L$ is the linearized operator relating the model $\bold m$and the data $\bold D$ Claerbout (1999).

Common industrial practice is to recover the model by applying the adjoint of the operator $\mathcal L$ to the recorded data

$$\bold m\approx \mathcal L^* \bold D,$$

where $\mathcal L^*$, the adjoint operator, is an approximation to the operator we would really like to apply, $\mathcal L^{-1}$.

The adjoint approximation is not only very convenient to use, but it also appears to produce good results in practice. However, there are situations when, given that $\mathcal L$ is not at all unitary, the model we produce is only a distant approximation of the true one.

In the particular case of seismic imaging, the model (reflectivity map) and the data (seismic data) are extremely large. It becomes therefore, completely infeasible to compute $\mathcal L^{-1}$ other than by an iterative application of the operator $\mathcal L$ and its adjoint. Furthermore, the model and the data are far too large to be kept in the RAM memory of current computers.

A solution for the large size problems is to implement inversion in an out-of-core fashion, where only limited chunks of the model and data are kept in memory at any given time. This is the purpose of the optimization library I am introducing in this paper.

Generally speaking, the types of problems that can be solved using this library are regularized inversion in standard form  

$$\left [\mathcal W\mathcal L\atop\epsilon\mathcal A\right ] \bold m\approx \left [\mathcal W\bold D\atop0\right ],$$ (1)

or in its preconditioned form  

$$\left [\mathcal W\mathcal L\P \atop\epsilon I\right ] \bold p\approx \left [\mathcal W\bold D\atop0\right ],$$ (2)

where $\mathcal A$ is a regularization operator, $\mathcal W$ is a weighting operator, $\P$ is a preconditioning operator, and $\bold p$ is the preconditioned form of $\bold m$.

The operators $\mathcal L$, $\mathcal A$, $\mathcal W$, $\P$ are application-dependent and therefore have to be externally implemented, likely in an out-of-core manner, by the user of the library. All other operations needed to solve the inversion problem are build into oclib.

Although other languages allow for more creative implementation, this entire library is implemented in Fortran90 mainly because this is still the programming language of choice in scientific computing and also the language most commonly used at SEP Claerbout (1999).