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.