FAMILIAR OPERATORS

The operation $y_i = \sum_j b_{ij} x_j$ is multiplication of a matrix $\bold B$ times a vector $\bold x$.The adjoint operation is $\tilde x_j = \sum_i b_{ij} y_i$.The operation adjoint to multiplying by a matrix is multiplying by the transposed matrix (unless the matrix has complex elements, in which case we need the complex-conjugated transpose). The following pseudocode does matrix multiplication $\bold y=\bold B\bold x$ and multiplication by the transpose matrix $\tilde \bold x = \bold B' \bold y$:

		if operator itself
		 		then erase y
		if adjoint
		 		then erase x
		do iy = 1, ny {
		do ix = 1, nx {
		 		if operator itself
		 		 		y(iy) = y(iy) + b(iy,ix) $\times$ x(ix)
		 		if adjoint
		 		 		x(ix) = x(ix) + b(iy,ix) $\times$ y(iy)
		      }}

Notice that the ``bottom line'' in the program is that x and y are simply interchanged. The above example is a prototype of many to follow, so observe carefully the similarities and differences between the operation and its adjoint.

A formal program for matrix multiply and its adjoint is found below. The first step is erasing the output. That may seem like too trivial a function to put in a separate library routine, but, at last count, 15 other routines in this book use the output-erasing subroutine adjnull() below.  

subroutine adjnull( adj, add, x, nx,  y, ny )
integer ix, iy,     adj, add,    nx,     ny
real                          x( nx), y( ny )
if( add == 0 )
        if( adj == 0 )
                do iy= 1, ny 
                        y(iy) = 0.
        else
                do ix= 1, nx 
                        x(ix) = 0. 
return; end

The subroutine matmult() for matrix multiply and its adjoint exhibits a style that we will use repeatedly.  

# matrix multiply and its adjoint
#
subroutine matmult( adj, bb,        nx,x,  ny,y)
integer ix, iy,    adj,            nx,    ny
real                    bb(ny,nx), x(nx), y(ny)
call adjnull(      adj, 0,         x,nx,  y,ny)
do ix= 1, nx {
do iy= 1, ny {
        if( adj == 0 )
                        y(iy) = y(iy) + bb(iy,ix) * x(ix)
        else
                        x(ix) = x(ix) + bb(iy,ix) * y(iy)
        }}
return; end