Summary

This document describes the Fortran version of the code KronLinInv.

If you use this code for research or else, please cite the related paper:

Andrea Zunino, Klaus Mosegaard, An efficient method to solve large linearizable inverse problems under Gaussian and separability assumptions, Computers & Geosciences, 2018 ISSN 0098-3004, https://doi.org/10.1016/j.cageo.2018.09.005.

See the paper mentioned above for a detailed description.

Two versions of the code are provided:

kronlininv.f08, a serial or parallel version for shared memory (OpenMP), depending on chosen compiler options;
ompi_kronlininv.f08, a parallel version for distributed memory (OpenMPI).

The latter should be used for solving large size problems e.g. on a cluster computer, while the former is simpler and requires less dependencies. Shell scripts to compile the two versions are provided.

KronLinInv is thought as a set of subroutines, such that it can be embedded in other code. However, two standalone test programs are provided, test_kronlininv.f08 and test_ompi_kronlininv.f08, one for each version, which read input data from a HDF5 file (HDF5 library for Fortran must be installed) and write the results to another HDF5 file.

KronLinInv solves the linear inverse problem with Gaussian uncertainties represented by the following objective function:

\[ S( \mathbf{m}) = \frac{1}{2} ( \mathbf{G} \mathbf{m} - \mathbf{d}_{\sf{obs}} )^{\sf{T}} \mathbf{C}^{-1}_{\rm{D}} ( \mathbf{G} \mathbf{m} - \mathbf{d}_{\sf{obs}} ) + \frac{1}{2} ( \mathbf{m} - \mathbf{m}_{\sf{prior}} )^{\sf{T}} \mathbf{C}^{-1}_{\rm{M}} ( \mathbf{m} - \mathbf{m}_{\sf{prior}} ) \]

The posterior covariance matrix is given by

\[ \mathbf{\widetilde{C}}_{\rm{M}} = \left( \mathbf{G}^{\sf{T}} \, \mathbf{C}^{-1}_{\rm{D}} \, \mathbf{G} + \mathbf{C}^{-1}_{\rm{M}} \right)^{-1} \]

and the center of posterior gaussian is

\[ \mathbf{\widetilde{m}} = \mathbf{m}_{\rm{prior}}+ \mathbf{\widetilde{C}}_{\rm{M}} \, \mathbf{G}^{\sf{T}} \, \mathbf{C}^{-1}_{\rm{D}} \left(\mathbf{d}_{\rm{obs}} - \mathbf{G} \mathbf{m}_{\rm{prior}} \right) . \]

The paper describes how to transform the previous problem in order to obtain what follows.

\[ \mathbf{U}_1 \mathbf{\Lambda}_1 \mathbf{U}_1^{-1} = \mathbf{C}_{\rm{M}}^{\rm{x}} (\mathbf{G}^{\rm{x}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{x}})^{-1} \mathbf{G}^{\rm{x}} \]

\[ \mathbf{U}_2 \mathbf{\Lambda}_2 \mathbf{U}_2^{-1} = \mathbf{C}_{\rm{M}}^{\rm{y}} (\mathbf{G}^{\rm{y}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{y}})^{-1} \mathbf{G}^{\rm{y}} \]

\[ \mathbf{U}_3 \mathbf{\Lambda}_3 \mathbf{U}_3^{-1} = \mathbf{C}_{\rm{M}}^{\rm{z}} (\mathbf{G}^{\rm{z}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{z}})^{-1} \mathbf{G}^{\rm{z}} \]

The posterior covariance is expressed as

\[ \mathbf{\widetilde{C}}_{\rm{M}} = \left( \mathbf{U}_1 \otimes \mathbf{U}_2 \otimes \mathbf{U}_3 \right) \big( \mathbf{I} + \mathbf{\Lambda}_1 \! \otimes \! \mathbf{\Lambda}_2 \! \otimes \! \mathbf{\Lambda}_3 \big)^{-1} \big( \mathbf{U}_1^{-1} \mathbf{C}_{\rm{M}}^{\rm{x}} \otimes \mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{y}} \otimes \mathbf{U}_3^{-1} \mathbf{C}_{\rm{M}}^{\rm{z}} \big) \]

and the posterior mean model as

\[ \mathbf{\widetilde{m}} = \mathbf{m}_{\rm{prior}} + \Big[ \! \left( \mathbf{U}_1 \otimes \mathbf{U}_2 \otimes \mathbf{U}_3 \right) \big( \mathbf{I} + \mathbf{\Lambda}_1\! \otimes \! \mathbf{\Lambda}_2 \! \otimes\! \mathbf{\Lambda}_3 \big)^{-1} \\ \times \Big( \left( \mathbf{U}_1^{-1} \mathbf{C}_{\rm{M}}^{\rm{x}} (\mathbf{G}^{\rm{x}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{x}})^{-1} \right) \! \otimes \left( \mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{y}} (\mathbf{G}^{\rm{y}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{y}})^{-1} \right) \! \\ \otimes \left( \mathbf{U}_3^{-1} \mathbf{C}_{\rm{M}}^{\rm{z}} (\mathbf{G}^{\rm{z}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{z}})^{-1} \right) \Big) \Big] \\ \times \Big( \mathbf{d}_{\rm{obs}} - \big( \mathbf{G}^{\rm{x}} \otimes \mathbf{G}^{\rm{y}} \otimes \mathbf{G}^{\rm{z}} \big) \, \mathbf{m}_{\rm{prior}} \Big) \]

These last two formulae are those used by the KronLinInv.

Several routines are provided:

calcfactors(): Computes the factors necessary to solve the inverse problem
blockpostcov(): Computes a block (or all) of the posterior covariance
bandpostcov(): Computes a band of the posterior covariance
posteriormean(): Computes the posterior mean model

Author: Andrea Zunino
Niels Bohr Institute, University of Copenhagen