The next step was to write a few simple operators. We coded the four
basic convolution operators used at SEP: `tcai`, `tcaf`, `
icai`, and `icaf`. These acronyms indicate convolution in which
the boundary conditions are either transient or internal and in which
the adjoint gives either the filter or the input.

We further realized that `tcai` and `tcaf` can be combined
into a single bilinear operator. A bilinear operator takes two vector
arguments and returns a vector value which is linear in each of the
arguments considered by itself. We wrote a transient convolution
class which is a subclass of HCL's bilinear operator class. In so
doing, we unified the code for `tcai` and `tcaf`. (Likewise,
we unified `icai` and `icaf` into a single internal
convolution bilinear operator). Jon Claerbout stresses that the key
code segments of an operator's forward and adjoint should be as close
together as possible. However, by writing both `tcai` and `
tcaf` as separate routines, he duplicates the forward operation and
isolates the two adjoints in separate routines. Our transient
convolution operator encapsulates the forward and both adjoints in a
single routine. The following is pseudo-code for the key part of our
transient convolution bilinear operator in one dimension:

% Pseudo-code for transient convolution as a bilinear operator (in 1D)do fi = 0,m-1 % Loop over the filter components (index begins at 0) do xi = 0,n-1 % Loop over the input components yi = xi + fi % Calculate the corresponding output components if (result = output) then y(yi) += f(fi) * x(xi) % forward if (result = input ) then x(xi) += y(yi) * f(fi) % adjoint of tcai if (result = filter) then f(fi) += x(xi) * y(yi) % adjoint of tcaf end do end do

Another feature of our convolutions is that they are written recursively. This allows convolution in any dimension without having to write separate routines. Recursion, however, is usually a speed bottleneck. An alternative would be to write a single iterative (i.e. non-recursive) version in a sufficiently large dimension. That is, because of our indexing mechanism of RGFs in C++, one could use (iterative) convolution code written for ten dimensions on a filter, input, and output of only two dimensions.

Furthermore, we have envisioned a more general convolution operator in which the boundary conditions (transient or internal) are specified for each dimension individually. We believe that this would be easy to implement with a recursive algorithm, although we are not sure if there is a demand for such an operator. Even more general would be a convolution in which the boundary conditions are specified at each boundary; however, making the code elegant and efficient might be difficult. Nonetheless, we could write it in C++ whereas in Fortran, because of restrictions on recursion, it would be quite difficult.

In addition to convolutions, we wrote masking operators, which zero specified components of an RGF. These are important in some of Claerbout's missing data problems.

11/11/1997