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.
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
module dottest
use random_mod
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.