Next: EXAMPLES OF LINEAR OPERATORS Up: Fomel & Claerbout: F90 Previous: Introduction

# UNIVERSAL DOT-PRODUCT ROUTINE

Dot-product routine is a numerical implementation of the definition of adjoint operator. We use it to confirm that the dot product of vectors and is equal (up to the machine precision) to the dot product of and for random and .

The dot-product test is crucially important for applications, where the composite operators or need to be positive definite. An important example is the conjugate gradient optimization. As a rule, an unusually high misfit in the dot-product test indicates either a conceptual error or a programming bug.

The dot-product test is an example of a mathematical abstraction, which requires object-oriented programming style to be expressed adequately. With Fortran 77 tools, we have to write a separate dot-product test program for each operator. Each time the dot-product test fails, we need to check both the operator implementation, and the implementation of the test itself. Ideally, the dot-product test routine should be isolated from the operator and precompiled in a library. The concept of function interface, introduced in Fortran 90, allows us to achieve this goal in a straightforward way.

Module dottest is the actual implementation of a universal routine. An operator, passed to the dot_test subroutine is defined by its interface

function oper (adj, add, mod, dat) result (stat)
integer              :: stat
real, dimension (:)  :: mod, dat
end function oper

When the logical argument adj is .false., the operator maps the model vector mod (defined as a one-dimensional real array) to the data vector dat. When adj=.true., dat is mapped into mod. When add is .true., the output is added to the previous value of the corresponding vector.

module dottest
use random_mod

contains

subroutine dot_test (oper, n_mod, n_dat, dotprod, dotadd)
integer,             intent (in)  :: n_mod, n_dat
real, dimension (2), intent (out) :: dotprod, dotadd
interface
integer              :: stat
real, dimension (:)  :: mod, dat
end function oper
end interface

real, dimension (n_mod)           :: mod1, mod2
real, dimension (n_dat)           :: dat1, dat2
integer                           :: stat
logical, parameter                :: T = .true., F = .false.

call random (mod1) ; call random (dat2)

stat = oper (F, F, mod1, dat1) ; dotprod (1) = dot_product (dat1, dat2)
stat = oper (T, F, mod2, dat2) ; dotprod (2) = dot_product (mod1, mod2)

write (0,*) dotprod ; if (.not. present (dotadd)) return

stat = oper (F, T, mod1, dat1) ; dotadd  (1) = dot_product (dat1, dat2)
stat = oper (T, T, mod2, dat2) ; dotadd  (2) = dot_product (mod1, mod2)