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 $\bold y$ and ${\bf A x}$ is equal (up to the machine precision) to the dot product of $\bold x$ and ${\bf A' y}$ for random $\bold x$ and $\bold y$.

The dot-product test is crucially important for applications, where the composite operators ${\bf A' A}$ or ${\bf A A'}$ 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

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
       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.