Applying a block operator to composite IGF90 objects

Blocks of IGF90 operators can also be easily implemented within the HCL mainframe. The following example implements a basic block operator that gets applied to a composite IGF90 object:

#include "HCL_BlockOp.h"
#include "HCL_ProductVector.h"
char HCL_ProductVectorSpace::Type;
  ...

  // Create data vector and its space
  IGF90 data( "data" );
  IGF90Space *Dspace = &( (IGF90Space&) data.Space() );

  // Build axes
  Axis a0( 0. , 1.,  "trace", "number" , 1 );
  Axis a1( -112.85 , 0.001,  "sx", "meters" , 50 );
  Axis a2( -16.012501 , .00025,  "sy", "meters" , 150 );

  const int NN=3; AxisArray aa(NN);
  aa[0] = a0; aa[1] = a1; aa[2] = a2;

  IGF90 y( NN, aa );
  y.Random();   // model some noise

  // Create model vector and its space
  IGF90 x( NN, aa );
  IGF90Space *Sspace = &( (IGF90Space&) x.Space() );
  x.Zero(); 

  // Create a 2x1 left-hand-side (lhs) block vector
  HCL_ProductVector lhs(&data,0,  &y,0);
  HCL_ProductVectorSpace *Lhsspace = &((HCL_ProductVectorSpace &) lhs.Space());

  // Create the 2x1 block operator
  HCL_BlockLinearOpAdj N( Sspace,  Lhsspace );
  N.SetNext( &L, 0 );                   //      N = | L |
  N.SetNext( &C, 0 );                   //          | C |

  //                                               |data|
  // Try the adjoint block operator  x =  |L'  C'| |    |
  //                                               | y  |
  N.Adjoint().Image( lhs, x );
  x.Write();

The lhs vector of the HCL_ProductVector class implements a $2\times 1$ matrix of HCL_vectors. The N object of the HCL_BlockLinearOpAdj class implements a $2\times 1$ block operator. The shape of the domain and the range vector spaces determines the shape of N operator. The SetNext() member function fills the element of the block operator by rows, from top to bottom.

Adding the next couple of lines to the previous example implements the dot product test for the block operator:

  IGF90 xp( *Sspace );  IGF90 dp ( *Dspace );
  IGF90 yp( *Sspace );

  HCL_ProductVector lhsp(&dp,0,  &yp,0);

  data.Random(); y.Random(); x.Random();
  lhsp.Zero();   xp.Zero();

  // Applied the Forward mapping 
  N.Image( x, lhsp );
  // Applied the Adjoint mapping 
  N.Adjoint().Image( lhs, xp );

  // Dot product for block operator
  cerr << "  " << lhsp.Inner(lhs)<< " " << xp.Inner(x) << endl;

Figure 3 shows the result of the block operation. The block operator interpolates the Seabeam data and adds to it some background filter noise. The filter used is a Laplacian filter.

 

block Figure 3 Row operation of interpolation of Seabeam dataset and convolution of random noise.