Basic operators and adjoints

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}\T \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 preceding example is a prototype of many to follow; therefore observe carefully the similarities and differences between the operator and its adjoint.

Next, we restate the matrix-multiply pseudo code in real code.

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

user/pwd/matmult.c

void matmult_lop (bool adj, bool add, 
		  int nx, int ny, float* x, float*y) 
/*< linear operator >*/
{
    int ix, iy;
    sf_adjnull (adj, add, nx, ny, x, y);
    for (ix = 0; ix < nx; ix++) {
		for (iy = 0; iy < ny; iy++) {
			if (adj) x[ix] += B[iy][ix] * y[iy];
			else     y[iy] += B[iy][ix] * x[ix];
		}
    }
}

We now have a module with two entries/ One is named _init for the physical scientist who defines the physical problem by defining the operator. The other is named _lop for the least-squares problem solvers, computer scientists not interested in how we specify $\bold B$. They will be iteratively computing $\bold B\bold x$ and $\bold B\T \bold y$ to optimize the model fitting.

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

           matmult_init( bb);

is done by physical scientists to set up the operation. Most later calls are done by numerical analysts in solving code like in Chapter . These calls look like

            matmult_lop( adj, add, nx, ny, 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$.

We split operators into two independent processes; the first is used for geophysical set up, while the second is invoked by mathematical library code (introduced in the next chapter) to find the model that best fits the data. Here is why we do so. It is important that the math code contain nothing about the geophysical particulars. This independence enables us to use the same math code on many different geophysical applications. This concept of ``information hiding'' arrived late in human understanding of what is desirable in a computer language. Subroutines and functions are the way new programs use old ones. Object modules are the way old programs (math solvers) are able to use new ones (geophysical operators).

Basic operators and adjoints