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 operWhen the logical argument

module dottest use random_modcontains

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.

11/11/1997