Privacy and Security Notice

Tutorial 2: Handles, linear operators and “factories”

Introduction

This tutorial exploits the versatility of Chroma to define objects, and “factories”

Exercise 1: Handles and LinearOperators

In this exercise we will create a basic Wilson linear operator and do things with it. Here we will work directly with the objects. In the next exercise we will play with the factories.  Work in the same directory as tutorial 1, and go back to our file tut1.cc.

Find the comment:

// Do something exciting here yourself
// Suggested exercise: Create a linear operator and apply it
// to a gaussian noise vector.

and we will work below this

Getting hold of a linear operator involves the following steps:

  • First we must specify our boundary conditions
  • Then we must make a CreateFermState object. This is an indirection so that we can use the same fermion actions for both thin link and smeared link states. Information about the boundaries is buried in the CreateFermState
  • We can then create a FermionAction using the CreateFermState. This object will allow us to create fermionic states and linear operators.
  • We use the FermionAction to create a gauge field state that can bhave a LinearOperator built over it (ie one that may be smeared/stouted and have the boundaries applied to it).
  • We then use the FermionAction to create the linear operator over the state which we just created.

This chain seems elaborate but really it is quite logical. The basic idea is that once we built the fermion action, it can take a gauge field and modify it appropriately to build the linear operator over. Once this is done, we can create the linear operator. So this is essentially the last two lines. We also need to specify the boundary conditions, clearly, so that is the first line. But why this messing around with CreateFermState... Well it didn't use to be this way until we realised that we may want to build fermion operators over stout links say, but also to keep the thin links for the gauge actions. So we'd have to have either

  • Stout the gauge fields first OR
  • Copy the whole fermion action changing only its createState function
  • Do what we chose to do: put in an extra indirection.

The extra indirection being the CreateFermState intermediary. This is a base class, whose derived classes can create either simple or stouted (and in the future perhaps other) states.

Secondly we should mention that FermBC, CreateFermState, FermAct and LinearOperator all take template parameters. We use the following 3 template parameters:

  • T the type of the FermionField (LatticeFermion)
  • P the type for the LieAlgebra (momenta) - even differentiable linear operators need to know this - the time derivatives (forces) live in the Lie Algebra. We typically use multi1d<LatticeColorMatrix>
  • Q the type of the Group (gauge fields), again we typically use multi1d<LatticeColorMatrix>

Typing out the full types (can be done and is done in many places throughout chroma) is unwieldy, with long lines and lots of angle brackets. We can save ourselves some tyding and make the code readable by defining aliases to these types. We can do this anywhere in the code but only once within one scoping unit. Let's do it now. Add the following typedef lines into the file:

typedef QDP::LatticeFermion T;
typedef QDP::multi1d<LatticeColorMatrix> P;
typedef QDP::multi1d<LatticeColorMatrix> Q;

Now we need to specify the boundary conditions and we need a FermBC<T,P,Q> object for this. Recall that this is a virtual base class, so we cannot instantiate it directly. We have to create its subclass and keep a reference (or in this case a reference counted handle) to it. To keep things simple, we'll use a SimpleFermBC<T,P,Q> which allows periodic, antiperiodic and Dirichlet boundaries. We'll use periodic for the first 3 dims and antiperiodic for the 4th (time) dim.

// Parameters for boundary conditions
SimpleFermBCParams bc_params;

bc_params.boundary.resize(Nd);
for(int i=0; i < Nd-1; i++) {
   bc_params.boundary[i] = BC_TYPE_PERIODIC;
}
bc_params.boundary[Nd-1] = BC_TYPE_ANTIPERIODIC;

// The boundary conditions themselves - we will get one from the
// heap using new and drop it into the handle, since that is what
// the CreateFermState constructor wants.
Handle< FermBC<T,P,Q> > fbc_handle( new SimpleFermBC< T,P,Q >(bc_params) );

So here we got an instance of the derived class (SimpleFermBC) and returned a pointer to it. The pointer got down-cast to be a pointer to the base class (FermBC) which we used in the constructor of the Handle. We don't need to worry about freeing this. When the Handle goes out of scope and the reference count drops to 0, it will be disposed of automatically.

Now we need to create a CreateFermState< T,P,Q > object which can be used to apply the Boundaries we just created to a simple thin field. Again to keep things simple, we'll use what is called a CreateSimpleFermState class which is a derived class of CreateFermState that only applies boundaries, and does not do things to the links otherwise. We follow a similar pattern to creating the SimpleFermBC. Create the base class from the heap using new, down cast the pointer to the base class, and drop it into the a handle. All this can be done by the line:

// Create a FermState Creator with boundaries
Handle<CreateFermState<T,P,Q> > cfs( new CreateSimpleFermState<T,P,Q>(fbc_handle));

Now we have enough information (objects) to create a FermionAction. We will start off with an EvenOddPreconditioned WilsonFermion action:

// Now create the FermAct using the FermState Creator
EvenOddPrecWilsonFermAct S( cfs, Real(0.3));

Here the Real(0.3) is the desired quark mass which is the only real parameter for this action.

Let us start working with the FermionAction. We need to take our gauge field and create a state from it which has the boundaries applied:

// Use the FermAct to create a suitable FermState
Handle< FermState<T,P,Q> > fs( S.createState(u) );

Now we can create a linear operator over this state:

// Use the FermAct to create a linear Operator over the FermState
Handle< LinearOperator<T> > M( S.linOp(fs) );

Now that we have a Linear Operator. We need to apply it to somethig. We will use a LatticeFermion filled with gaussian noise. The code for initialising the vectors and applying the operator is below.

// Declare the vectors
LatticeFermion source_vector;
LatticeFermion target_vector;

// fill the source with gaussian noise
// We use the subset method of the linear operator
// so we only fill the part of the vector the operator cares about

gaussian(source_vector, M->subset());
target_vector[ M->subset() ] = zero;

// Apply the linear operator
(*M)(target_vector, source_vector, PLUS);

Finally we should do something we can check. One thing we can check is that applying the linear operator and its dagger in succession gives the same result as applying M^\dagger M. We can get the "squared operator" by using the lMdagM() function in the fermion action. Since this is a 'test' we will make it first fail by leaving a deliberate bug in the code. The code is as follows. The comments explain the intention. Imagine it continuing on from before:

// We already have M source.
// Now apply M^{\dagger} onto the target
// to get (M^\dagger M) source_vector

LatticeFermion mdagm_target_1;

// DELIBERATE MISTAKE TO TEST THE TEST: Use PLUS instead of MINUS
(*M)(mdagm_target_1, target_vector, PLUS); // MINUS means apply M^\dagger

// Now check against M^\dagger M method.
// Create the linear operator
Handle< const LinearOperator<LatticeFermion> > MdM( S.lMdagM(fs) );

// A target for the MdM operator
// I will not initialize it as I would just overwrite the initialisation
LatticeFermion mdagm_target_2;

// Apply the "squared" operator
(*MdM)(mdagm_target_2, source_vector, PLUS);

// A vector for the difference of the 2 results
LatticeFermion diff;

// Take the difference on the subset of the linear operator
// Note: Subset index is only on the target
diff[M->subset()] = mdagm_target_2 - mdagm_target_1;

// Compute the norm of the difference on the subset
// norm2() computes the squared norm.
Double norm_diff = sqrt(norm2(diff, M->subset()));

// Print the difference.
QDPIO::cout << "LinearOperator test: norm2(diff) = " << norm_diff << endl;

// Check that the difference is small
// Note: I use the toBool() function because I am comparing QDP++
// types rather than C++ primitive types and I have not overloaded

if( toBool( norm_diff > Double(1.0e-5) ) ) {
  QDPIO::cout << "LinearOperator test FAILED" << endl << flush;
  QDP_abort(1);
}
else {
  QDPIO::cout << "LinearOperator test PASSED" << endl << flush;
}

When you compile and run the program you should see telling you that the test has FAILED. This of course is unsurprising since we put in a deliberate bug to check that the test does its job. Correct the mistake and make the test pass, by changing the incorrect PLUS into a MINUS so that the second application of M applies the daggered operator.

What you have written is essentially a unit test. This is different from a regression test as it checks expected behaviour rather than that previous code has not been broken.

Exercise 2: Factories

In the previous exercise we explicitly created the boundary condition parameters, boundary conditions and the fermion action. In this exercise we will use a factory to make this process more generic. This will

  • Allow us to switch fermion actions using an input XML file
  • Allow us to use stout smearing
  • Allow us to use different boundaries using an input XML file
  • See how to use the factories

Find the line in the file after the

// Do something exciting here yourself
// Suggested exercise: Create a linear operator and apply it
// to a gaussian noise vector.
typedef QDP::LatticeFermion T;
typedef QDP::multi1d<LatticeColorMatrix> P;
typedef QDP::multi1d<LatticeColorMatrix> Q;

and delete the code from below this all the way until you reach the line

// Use the FermAct to create a suitable FermState

We will insert the code here, essentially creating the same fermion action S by a different means. We can add in the following piece of code

// I define this up front because I can't do it within the
// try{} catch{} block.
// I initialise it to 0.

WilsonTypeFermAct<T,P,Q>* ferm_act_pointer = 0;
// Try and read a fermion action definition
try {
  // Get the name of the fermion action
  std::string fermact_name;
  read(xml_in, "/tutorial3/FermionAction/FermAct", fermact_name);

  // Get the factory to create the FermionAction
  ferm_act_pointer =
    TheWilsonTypeFermActFactory::Instance().createObject(fermact_name,
             xml_in,
             "/tutorial3/FermionAction");
}
// Handle any errors
catch( const std::string& e) {
  QDPIO::cerr << "Caught Exception: " << e << endl << flush;
  QDP_abort(1);
}

// Drop the pointer into a handle - so that it gets cleaned up later
Handle< WilsonTypeFermAct<T,P,Q> S_handle(ferm_act_pointer);

// We could just use the handle, but in the previous code we used
// the object. We can emulate this behaviour by creating a reference
// to the dereferenced handle.

WilsonTypeFermAct<T,P,Q>& S = *S_handle;

What has happened? We have replaced the concrete invocation of the Wilson Fermion action with a more generic class: WilsonTypeFermAct. So this code should work now for all available WilsonTypeFermActs that are registered with our factory.

Secondly, we are reading the FermionAction tag from the XML file. This contains the FermState tag and the FermionBC tag which specify the wanted subclasses of CreateFermState and FermBC respectively. So we can now also run this program with stout states and other FermionBCs.

Before running the program we need to add the FermionAction we want to the input XML file: Within the first group (ie between the <tutorial3> </tutorial3> tags, below the <Prop> tag) add in the following snippet:

<FermionAction>
   <FermAct>WILSON</FermAct>
   <Mass>0.3</Mass>
   <FermState>
      <Name>SIMPLE_FERM_STATE</Name>
      <FermionBC>
         <FermBC>SIMPLE_FERMBC</FermBC>
         <boundary>1 1 1 -1</boundary>
      </FermionBC>
   </FermState>
</FermionAction>

Note that the fermion action says WILSON but in principle you now switch this to CLOVER or PARWILSON (twisted mass) with their extra parameters. So in principle your code can now unit test the MdagM method of all Wilson-like fermions. Likewise you can change the FermState as well (eg to STOUT_FERM_STATE) and of course you can also vary teh boundary types.

If you run this program you may find that it ends in dismal FAILURE, with the message:

Couldnt find key WILSON in the map:
Available Keys are :
Caught Exception: Factory error: unknown identifier: id = WILSON

This is the problem of linkage. I mentioned in the talk on Design Patterns in Chroma. We need to generate linkage for the fermion actions in the factory. If we don't do this, all our nice Fermion Actions, Boundary Conditions etc do not get registered in the respective factories.

To register all the actions and measurements we use in chroma, we need to add the line

WilsonTypeFermActs4DEnv::registerAll();

somewhere before the line:

  ferm_act_pointer =
    TheWilsonTypeFermActFactory::Instance().createObject(fermact_name,
             xml_in,
             "/tutorial3/FermionAction");

This will bring in linkage for all the WilsonType FermionActions and their dependencies, and should make the test pass. In the main programs this kind of linkage hack is moved to an explicit linkageHack() function.

Exercise 3: Play Away

You can now do lots of things and should play:

  • Insert Different Gamma Matrices instead of gamma_5, and generate other mesons
  • Increase the maximum allowed momenta in SftMom
  • Turn off averaging over momenta when you create SftMom
  • Try different linear operators. Look in chroma/tests/chroma/hadron/propagator for input files.
  • Try using a stout smeared state. Here is one with periodic boundaries in space and antiperiodic in time:

<FermState>
   <Name>STOUT_FERM_STATE</Name>
   <rho>0.22</rho>
   <n_smear>3</n_smear>
   <orthog_dir>3</orthog_dir>
   <FermionBC>
         <FermBC>SIMPLE_FERMBC</FermBC>
         <boundary>1 1 1 -1</boundary>
   </FermionBC>
</FermState>

This should do 3 levels of stout smearing with an isotropic coefficient rho = 0.22 in space only (orthog_dir=3 excludes smearing in the time direction.