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 $y_i = \sum_j b_{ij} x_j$ to multiply a vector $\bold x$ by a matrix $\bold B$ and the operation $\tilde x_j = \sum_i b_{ij} y_i$ to multiply the data vector $\bold y$ by the transpose matrix $\bold B'$. 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

subroutine adjnull (adj, add, x, y)

logical, intent (in) :: adj, add
    real, dimension (:)  :: x, y

if (.not. add) then
       if (adj) then 
          x = 0.
       else 		    
          y = 0.
       end if
    end if

end subroutine adjnull

end module adj_mod

module matmult
  use adj_mod

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

contains

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

function matmult_op (adj, add, x, y) result (stat)
    integer              :: stat
    logical, intent (in) :: adj, add
    real, dimension (:)  :: x, y

integer              :: i

call adjnull (adj, add, x, y)

do i = 1, size (x)
       if (adj) then
          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.