Programming linear operators

The operation $y_i = \sum_j b_{ij} x_j$ is the multiplication of a matrix $\bold B$ by a vector $\bold x$.The adjoint operation is $\tilde x_j = \sum_i b_{ij} y_i$.The operation adjoint to multiplication by a matrix is multiplication 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 $\tilde \bold x = \bold B' \bold y$:

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

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 adjoint and the operator itself.

Next we restate the matrix-multiply pseudo code in real code, in a language called Loptran, a language designed for exposition and research in model fitting and optimization in physical sciences.

The module matmult for matrix multiply and its adjoint exhibits the style that we will use repeatedly. At last count there were 53 such routines (operator with adjoint) in this book alone.  

#$
#$=head1 NAME
#$
#$matmult - matrix multiplication
#$
#$=head1 SYNOPSIS
#$
#$C<call matmult_init(bb)>
#$
#$C<ierr=matmult_lop(adj,add,model,data)>
#$  
#$=head1 PARAMETERS
#$  
#$=over 4
#$  
#$=item  bb  C<real(:,:)>   
#$  
#$       Matrix 
#$  
#$=item adj,add,model,data  
#$  
#$      Standard operator interface
#$  
#$=back
#$  
#$=head1 DESCRIPTION
#$  
#$Matrix multiplication
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut
#$
module matmult {   # matrix multiply and its adjoint
real, dimension (:,:), pointer :: bb 
#%  _init( bb)
#%  _lop(  x, y) 
integer ix, iy
do ix= 1, size(x) {
do iy= 1, size(y) {
        if( adj)
		x(ix) = x(ix) + bb(iy,ix) * y(iy)
        else	
		y(iy) = y(iy) + bb(iy,ix) * x(ix)
        }}
}

Notice that the module matmult does not explicitly erase its output before it begins, as does the psuedo code. That is because Loptran will always erase for you the space required for the operator's output. Loptran also defines a logical variable adj for you to to distinguish your computation of the adjoint x=x+B'*y from the forward operation y=y+B*x. In computerese, the two lines beginning #% are macro expansions that take compact bits of information which expand into the verbose boilerplate that Fortran requires. Loptran is Fortran with these macro expansions. You can always see how they expand by looking at http://sepwww.stanford.edu/sep/prof/gee/Lib/.

What is new in Fortran 90, and will be a big help to us, is that instead of a subroutine with a single entry, we now have a module with two entries, one named _init for the physical scientist who defines the physical problem by defining the matrix, and another named _lop for the least-squares problem solver, the computer scientist who will not be interested in how we specify $\bold B$, but who will be iteratively computing $\bold B\bold x$ and $\bold B' \bold y$ to optimize the model fitting. The lines beginning with #% are expanded by Loptran into more verbose and distracting Fortran 90 code. The second line in the module matmult, however, is pure Fortran syntax saying that bb is a pointer to a real-valued matrix.

To use matmult, two calls must be made, the first one

            call matmult_init( bb)

is done by the physical scientist after he or she has prepared the matrix. Most later calls will be done by numerical analysts in solving code like in Chapter . These calls look like

            stat = matmult_lop( adj, add, x, y)

where adj is the logical variable saying whether we desire the adjoint or the operator itself, and where add is a logical variable saying whether we want to accumulate like $\bold y \leftarrow \bold y+\bold B\bold x$or whether we want to erase first and thus do $\bold y \leftarrow \bold B\bold x$.The return value stat is an integer parameter, mostly useless (unless you want to use it for error codes).

Operator initialization often allocates memory. To release this memory, you can call matmult_close() although in this case nothing really happens.