DOT PRODUCT TEST OF ADJOINT OPERATORS

The adjoint, A, of an operator, F, is generally defined as satisfying for all vectors and (in the appropriate vector spaces). (Here, the dot denotes the dot product.) SEP has traditionally used this definition to test if its operators are coded correctly. Thus, if the implementation of an operator's forward and adjoint satisfy the equality when and are random, we trust that they are indeed adjoints.

We wrote the dot product test in HCL. It can test any HCL_LinearOpAdj. Unfortunately, the user must supply the random vectors, since there is no provision for a random vector in HCL. There is no general procedure mathematically to choose a random vector from a vector space in general. However, we provided RGF with the method Random(), which fills the RGF with random numbers. Therefore, we could write a dot product test which requires only the operator, as long as it acts on and returns RGFs.

/*
  DotProductTest checks if a HCL_LinearOpAdj is coded correctly in the
  sense that the implementations of the forward and adjoint are actually
  adjoints of each other.

Since there is no general random vector concept in HCL_Vector,
  the user should provide random vectors w and v
  in the domain and range of FWD, respectively.
*/

void DotProductTest(       ostream         & ostr, // stream to send output
                           HCL_LinearOpAdj & FWD , // operator to test
                     const HCL_Vector      & w   , // random vector in domain
                     const HCL_Vector      & v   , // random vector in range
                     int   verbose               ) // flags output verbosity
{
  if (verbose) {                           // if verbose: display FWD, w, v
    ostr << "DOT PRODUCT TEST"
         << "--------"
         << "Operator (A):" << endl;
    FWD.Inspect(ostr);
    ostr << "Vector 1 (w):" << endl;
    w.Inspect(ostr);
    ostr << "Vector 2 (v):" << endl;
    v.Inspect(ostr); }

HCL_Vector * FWDw = v.Space().Member();  // allocate space for (FWD )w
  HCL_Vector * ADJv = w.Space().Member();  // allocate space for (FWD')v

FWD.          Image(w, *FWDw);           // calculate (FWD )w
  FWD.Adjoint().Image(v, *ADJv);           // calculate (FWD')v

float vFWDw = v.Inner(*FWDw);            // calculate v'(FWD )w
  float wADJv = w.Inner(*ADJv);            // calculate w'(FWD')v
                                           // display results
                                           // ideally, v'(FWD)w = w'(FWD')v

ostr << "RESULTS:"                       << endl;
  ostr << "<v, A (w)> = " << vFWDw         << endl;
  ostr << "<w, A'(v)> = " << wADJv         << endl;
  ostr << "Difference = " << vFWDw - wADJv << endl;  // should be small
}