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
logical, intent (in) :: adj, add
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
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 function oper (adj, add, mod, dat) result (stat) integer :: stat logical, intent (in) :: adj, add real, dimension (:) :: mod, dat end function oper end interface optional :: dotadd
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)
write (0,*) dotadd end subroutine dot_test
end module dottest
Before the dot_test program calls a particular operator, we need to supply the operator with the additional parameters it needs for execution. In order to do that, we encapsulate the operator function in one module with the initialization routine. The following section provides examples of modules for simple basic linear operators.