LINKING WITH FORTRAN77

We implemented the two-dimensional tcai operator in the HCL framework by calling Claerbout's Fortran77 routine tcai2. It was a simple exercise, and the relevant portions of code are shown below. (Notice that this example, being tcai alone and not combined with tcaf, does not emphasize the bilinearity of transient convolution. That is, the forward and only one of the two adjoints are coded together here.)

// if TCAIF77 is not defined, the C++ version is compiled
// if TCAIF77 is     defined, the F77 version is compiled
#ifdef TCAIF77
extern "C" void tcai2_(int &adj, int &add,
                       float *aa, int &na1, int &na2,
                       float *xx, int &nx1, int &nx2,
                       float *yy, int &ny1, int &ny2
                       );
#endif

void tcai2D::apply(int adj, RGF &xx, RGF &yy) const {
  int m1 = filter.GetSpace().GetNSample(0);
  int m2 = filter.GetSpace().GetNSample(1);
  int n1 =     xx.GetSpace().GetNSample(0);
  int n2 =     xx.GetSpace().GetNSample(1);
#ifdef TCAIF77

int l1 =     yy.GetSpace().GetNSample(0);
  int l2 =     yy.GetSpace().GetNSample(1);
  float *a = (float *) filter;
  float *x = (float *) xx;
  float *y = (float *) yy;
  int zeroAdd = 0;

tcai2_( adj, zeroAdd,
          a, m1, m2   ,
          x, n1, n2   ,
          y, l1, l2   );
#else

if (adj==0) { yy.Zero(); }
  else        { xx.Zero(); }

for (int a2=0; a2<m2; a2++){ for (int x2=0; x2<n2; x2++){ int y2 = x2+a2;
  for (int a1=0; a1<m1; a1++){ for (int x1=0; x1<n1; x1++){ int y1 = x1+a1;
    if (adj==0)
      yy(y1,y2) = yy(y1,y2) + xx(x1,x2) * filter(a1,a2);
    else
      xx(x1,x2) = xx(x1,x2) + yy(y1,y2) * filter(a1,a2);
  }}}}
#endif
}