Next: Transient convolution Up: EXAMPLES OF LINEAR OPERATORS Previous: EXAMPLES OF LINEAR OPERATORS

Matrix multiply

If both the input and the output of a linear operator are vectors, composed of real numbers, the operator is equivalent to multiplication by a real matrix. In data processing problems, we hardly ever see this matrix. We prefer instead to represent it in a compact non-matrix form and to write a subroutine that performs the forward operator and its adjoint. Nevertheless, multiplication by a matrix is a very important example of a linear operator, serving as the basic prototype for many other operators.

Module matmult implements both the operation to multiply a vector by a matrix and the operation to multiply the data vector by the transpose matrix . The matrix multiplication is preceded by a call to subroutine adjnull, defined in module adj, which optionally erases the output according to the value of the logical parameter add. When add=.true., the result of applying the operator is accumulated in the output, e.g. y=y+B*x.

module adj_mod
contains

real, dimension (:)  :: x, y

x = 0.
else
y = 0.
end if
end if



module matmult

real, dimension (:,:), pointer, private :: b

contains

subroutine matmult_init (matrix)
real, dimension(:,:), pointer :: matrix
b => matrix
end subroutine matmult_init

integer              :: stat
real, dimension (:)  :: x, y

integer              :: i

do i = 1, size (x)
x (i) = x (i) + sum (b (:, i) * y (:))
else
y (:) = y (:) +      b (:, i) * x (i)
end if
end do

stat = 0
end function matmult_op

end module matmult


Subroutine matmult_init serves as the object constructor, copying a pointer to the matrix in a private variable b, local with respect to the module, but accessible by the actual operator function matmult_op. The operator function returns an integer stat. In our convention, the value stat=0 corresponds to a normal execution of the operator routine. Otherwise, it can return the diagnostic of an error.

Fortran 90 has intrinsic functions matmul and dot_product for matrix-vector manipulations, which we could have used in this simple example. We prefer more explicit style to show how the algebraic symmetry of mathematical constructions translate into actual code. We will see more examples of this symmetry in the following examples.

Next: Transient convolution Up: EXAMPLES OF LINEAR OPERATORS Previous: EXAMPLES OF LINEAR OPERATORS
Stanford Exploration Project
11/11/1997