Tutorial 2: Programming with QDP++ and Chroma

Introduction

This tutorial takes you through the basics of programming with Chroma and linking against the Chroma library. It also shows you how we deal with XML in code.

Skill Level

This tutorial is rated intermediate. No heavy knowledge of C++ is needed. You can just "go with the flow".

Getting the files

Download the files for this tutorial by using from this link

Unpack the files with the following command:

gunzip tut2_files.tar.gz tar xvf tut2_files.tar

Now enter the working directory by typing:

cd Tutorial2

You should find several files in this directory:

input.xml

The input XML file for this tutorial

propagator_0

A test propagator for this tutorial

test_purgaug.cfg1

Your favourite old V=4x4x4x8 test lattice

Makefile

A Makefile to build the tutorial main program

tut2.cc

A skeleton Chroma application we will use during this tutorial

The `Makefile`

First let us look at (and edit) the Makefile. Bring it up in emacs

The first line defines a Makefile variable called CHROMA, containing the path to the chroma installation. We need to fix this so that it points to the installation on tutorial machines. You need to substitute the your chroma installation directory here.
The second line sets the location of the chroma-config program. This is based below the path defined in the CHROMA variable so you don't need to change it
The next four lines define Makefile variables:
- CXX -- The C++ compiler used to compile chroma
- CXXFLAGS -- The compiler flags to compile chroma
- LDFLAGS -- Additional linker flags to link chroma executables
- LIBS -- external libraries that need to be linked with Chroma
In the case of the CXXFLAGS we put also the option -I. instructing the compiler to look in the local directories for header files too. In this tutorial we don't have extra header files so this is a little redundant.

Note that these variables are set using commands of the form:
$(shell $(CONFIG) --cxx)
This instructs make to execute the command line following the string shell. In this case it will run the chroma-config program with various flags to determine characteristics of the installation needed to work with Chroma. If you change the CHROMA variable, this Makefile should work even on the QCDOC as long as you use GNU Make.

You can use this command yourself without putting it in a Makefile. Try typing:
< your chroma installation directory >/bin/chroma-config
to get a list of options. Now try to use some of them:
< your chroma installation directory >/bin/chroma-config --cxx
We have an empty HDRS variables, since we don't have extra headers.
We have a variable called OBJS that defines the object files used in our exercise. Currently we have only one object file: tut3.o which will be created from the source file tut3.cc.
We have the first rule of our makefile:
tut2: $(OBJS) $(CXX) -o $@ $(CXXFLAGS) $(OBJS) $(LDFLAGS) $(LIBS)
Hidden at the start of the second line is a TAB character.

This rule tells the makefile that we make target tut2 which depends on the files listed in the OBJS variable, by calling the command in the CXX macro with the following arguments. The $@ character is expanded to the name of the target, and the other variables are also expanded to their earlier definitions

Basically, this rule tells make, how to build our test executable
We have a rule to make .o object files from .cc source files:
%.o: %.cc $(HDRS) $(CXX) $(CXXFLAGS) -c $<
Here the $< character specifies the input file, and as with all make rule, the processing instruction has a TAB character in front of it.
Finally we have a rule for cleaning up. This is the rule that gets executed when you type make clean

So this makefile, will first try to make tut2.o as it is a dependency of the tut2 target. It uses the implicit rule to make .o files from .cc files. Finally, it creates the executable using the rule for tut2. When you type make without arguments, the first target will be processed which in this case is to make tut2.

Compiling and Running the Tutorial

Once you have changed the CHROMA macro in line one, you can just type make to build the executable. You should see output similar to

g++ -I/home/balint/gnu/include -I/home/balint/install/chroma_scalar_single/include -O2 -finline-limit=50000 -I/home/balint/install/qdp_scalar_single/include -I/opt/local/include/libxml2 -I. -c tut3.cc In file included from /home/balint/install/chroma_scalar_single/include/chromaba se.h:15, from /home/balint/install/chroma_scalar_single/include/chroma.h:29, from tut3.cc:6: /home/balint/install/qdp_scalar_single/include/qdp.h:174:2: warning: #warning "Using scalar architecture" g++ -o tut3 -I/home/balint/gnu/include -I/home/balint/install/chroma_scalar_single/include -O2 -finline-limit=50000 -I/home/balint/install/qdp_scalar_single/include -I/opt/local/include/libxml2 -I. tut3.o -R/home/balint/gnu/lib -L/home/balint/gnu/lib -L/home/balint/install/chroma_scalar_single/lib -L/home/balint/install/qdp_scalar_single/lib -lchroma -lqdp -lXPathReader -lxmlWriter -lqio -llime -L/opt/local/lib -R/opt/local/lib -lxml2 -lz -lpthread -lm -lsocket -lnsl -lgmp

You can ignore the warning. It is simply (intentionally) informing you that you are doing a scalar build.

You can now run the executable like you do with the stock applications:

./tut2 -i input.xml -o output.xl

and you should see the following output:

Hello World Finished init of RNG Finished lattice layout Starting up unit gauge (free) config Single Precision Read QIO_read_finished

and an output.xml file should have appeared in your local directory which we will inspect later.

So far you should have (apart from learning about makefiles) learned how you can use the chroma-config program to get the compiler and linker flags to link an application. If you want to do your own little projects you now can. We have shown you a simple makefile which can compile chroma applications

Hacking Chroma

Bring the tut2.cc file up in emacs and let us look through it. Lines beginning with // or text between /* and */ are comments.

First we include the file chroma.h which contains many definitions. It's location is given to the compiler by the -I flags given by chroma-config --cxxflags.

We tell the compiler that we are using the Chroma namespace. This is there so that other code can use the same names we do but in their own namespace

At this point it becomes useful to follow the comments in the code, and I'll just give a cursory overview of what is happening in the code. If you have done any C/C++ work this should be easy.

On line 18 we define a structure called Param_t which will be used to hold parameters. Currently the only parameter is

multi1d<int> nrow;

The multi1d is a C++ templated type for 1d arrays. It can have its size reset. The template parameter int between the < and > signifies that this will be an array of integers. It will be called nrow which is the conventional name for an array that holds the lattice size.

On line 25, we define a struct called Prop_t which will hold in itself a single string, the name of the propagator file. This string will be called prop_file and has C++ type std::string (or string from the std namespace).

On line 32 we define the complete input structure for the code called App_input_t. This consists of the previous Param_t, Prop_t structs and a struct called Cfg_t which is defined already in the chroma library.

Reading XML

We need to define XML reading routines for the structs we have just defined. This is done in the functions on lines 41, 49 and 58. Note that the functions all have the same name: read. The first two inputs are the same in all three cases but the third argument is different. This is a technique called function overloading in C++ and is a component of polymorphic programming.

Note that in the read routines, we specify as the first argument an object of type XMLReader& and const string& path. The XMLReader& is a reference to an XML document container, which allows us to read from the document. The path is an XPath expression. Consider the following XML snippet:

<?xml version="1.0"?> <foo> <bar> <fred> 5 </fred> </bar> </foo>

then we can select the value 5 using the Xpath expression:

/foo/bar/fred

We could do this in Chroma with the commands:

int integer; read(xml_reader, "/foo/bar/fred", integer);

Where the xml_reader object would be of type XMLReader and should contain the XML document in question. Because the reader function is declared as:

void read(XMLReader& reader, const std::string& path, int& i);

the &s mean that the arguments are passed by reference and are potentially modifiable (with the exception of the path which is declared const ie immutable).

We used the Xpath expression /foo/bar/fred. This is an absolute path. We can also do the equivalend of cd-ing in to the bar tag, and perform relative queries from there:

XMLReader reader2(xml_reader, "/foo/bar"); int integer; read(reader2,"fred",integer);

This technique is useful if a particular tag has several subtags of interest, and also to process individual groups. We use this trick on lines 41, 49 and 58 of our example code, to select the groups of interest for reading into individual structs.

Our technique, in lieu of data binding, is to build up recursive read functions using relative queries. This can be seen in around line 56, where we have just 3 reads to read in each part of the App_input_t structure, using the readers for the parts that we have just defined and are already defined in the chroma library

Exceptions

Occasionally things we may try and do are not possible. At times like this we need to interrupt what we are doing to report the error condition. This is done in C++ with the mechanism of exceptions. For example if we cannot open a file, or cannot match an XPath expression, we "throw" an exception. This exception can be caught by the program and handled gracefully. If it is not caught, then generally the program aborts.

For example, if the read in line 50 fails an exception will occur. This function does not handle the exception, so it will propagate up to the read defined on line 56, which is the only one to call the original read. We can see beginning on line 61 the following construct:

try { // Do stuff } catch( const string& e) { // Do this stuff if a string exception occurs }

What happens here is that the routines in the try block are executed, and if any of them throw an exception of type string that exception will be caught, and an appropriate error will be displayed. The XMLReaders only throw exceptions of type string. The strings contain error messages.

There can be other types of exceptions eg: file opening exceptions, casting exceptions, memory allocation exceptions that are more than just strings. Generally for each different type of exception you need a separate catch block. However there is a catch all case in which will catch all exceptions. An example of this is given on line 102

catch(...) // Catch all types of exceptions { // Handle exceptions here }

Main Program, Initialization and Ending

Once we have defined our XML readers, we get on to the main program proper on line 83, declared as:

int main(int argc, char **argv)

which is a standard declaration that should be familiar to you if you have ever used C. The argc is the number of command line arguments, and the argv is an array of C strings (char*s) holding the actual arguments.

On line 88 we have

Chroma::initialize(&argc,&argv)

which initialises QDP++ for us, and also grabs the input and output filenames out of the command line arguments.

Paired with this is a

Chroma::finalize()

at the end of the program at line 221

On line 89, there is a function call START_CODE() paired with a function END_CODE() on line 217. These functions invoke the profiling functionality of QDP++ if that has been enabled. You are encouraged to put START_CODE() and END_CODE() calls in later functions that you may write yourself.

The rest of the code

Lines 94-109 open our XML reader to read the input XML for the application.

We instantiate the App_input_t struct on line 94

We instantiate the XMLReader on line 100

We try to open the user specified input file or DATA between 101-109 catching errors as we go. Note the use of Chroma::getXMLInputFileName() -- read the comments for explanation

On line 112 we write Hello World to the terminal. Note we don't use the standard cout but rather the one defined in the QDPIO namespace:

QDPIO::cout << "Hello World" << endl;

The endl asks for a newline character. Using this special QDPIO::cout on a parallel machine will result only in 1 processor writing to the terminal.

Line 116, reads the app input structure from root tag <tutorial2>.

Lines 120 and 123 initialise the lattice layout from the XML that has been parsed into the App_input_t struct input making use of the nrow parameters:

Layout::setLattSize(input.param.nrow); // Initialise Layout::create();

This is something that always needs to be done before you want to do anythig lattice-y.

Lines 133-143 read the gauge configuration specified in the <Cfg> tag of the input XML file. We now commonly use XML headers in our configuration files. Further, if using the SciDAC format, there are also file related XML headers too. The gaugeStartup function returns these in the readers gauge_file_xml and gauge_xml. The actual gauge field has type:

multi1d<LatticeColorMatrix> u(Nd);

In other words it is a 1D array (with Nd elements) of LatticeColorMatrix objects. Each array element represents the SU(3) link fields in one of the Nd directions.

Finally the unitarity of the field is checked in line 141.

lines 152-163 read a SciDAC LIME format propagator from the disk. This is similar to reading the gauge fields, but with different function calls.

There is some amount of XML writing following all this, we will come back to it in the next section

On line 203 we measure the plaquette and link trace, after which there is room for you to code

Writing XML

QDP++ provides Chroma with three ways of writing XML:

XMLFileWriter -- These objects write to files
XMLBufferWriter -- These objects write to memory - ie to strings
XMLArrayWriter -- These can be used to write sequences seperated by <elem> tags

The writers support two kinds of operations:

push() and pop() operations. These open up and close group XML tags respectively. So

XMLFileWriter foo("my_file.xml"); // Opens a new XML file called my_file.xml push(foo,"Root"); pop(foo);

will open an XML file. The push will open a tag: <Root> and the pop will close it by writing </Root>.

write() -- These operations will write data to the XML file. The instruction

write(foo, "Barf", 5);

will write into the file controlled by the writer foo the following XML:

<Barf>5</Barf>

The operations can be combined to produce fully fledged XML documents:

XMLFileWriter foo("my_file.xml"); // Opens a new XML file called my_file.xml push(foo,"Root"); write(foo, "Barf", 5); pop(foo);

should produce:

<?xml version="1.0" ?> <Root> <Barf>5</Barf> </Root>

in the file my_file.xml.

Additionally, you can use XML writers to spew back the contents of readers. This will strip the readers of their XML headers.

XMLReader in("file.xml"); XMLFileWriter out("outfile.xml"); push(out, "Root"); write(out, "Fred", in); pop(out);

should spew the contents of file.xml into outfile.xml with a root tag of Root and surrounded by tags <Fred> </Fred>.

Some functions write their own tags, and do not need surrounding push and pop functions examples of this are:

proginfo() on line 174
MesPlq() on line 199

Sometimes the underlying file stream of an XMLFileWriter is buffered, so although you write to it, it may not immediately appear in the resulting file. If at this point your code crashes, it may not get written out at all. You can explicitly flush the streams associated with a writer by calling their member function flush. This may help you track more easily where crash-es have occurred, or make output from long calculations more "regular". You can see examples of this on lines 195 and 200 of the example program tut2.cc

At this point we have covered everything in the tut2.cc file and it is time for a few exercises

Exercise 1: XML Processing

Look at the output file we produced earlier: output.xml. Correlate the tags therein with the push(), pop, write, proginfo, MesPlq functions.

Exercise 2: More XML Processing

Look at the input file input.xml. Correlate the tags there, with the reader functions and XPath expressions

Exercise 3: More XML Processing

Make the read function for the Param_t struct read in your name from tags <my_name> and </my_name>.

You will need to add a member to the Param_t struct near line 18. This should be of type std::string. You will need to add a line near line 50, to read the new tag.

You will need to add the tag to the input.xml file.

Now remake the program and run it.

Exercise 4: Computing the pion correlator

Lets do some physics

Find the part of the program marked out with comments:

// Do something exciting here yourself. // Suggested exercise: Compute the zero mom pion on the propagator you // have read.

We will write our code below this. You should add in the code, into the file as we go along.

The program will have read for you a propagator into an object called quark_propagator. One thing we can do with this is to use it to compute the pseudoscalar operator:

C(t) = sum_{x} Tr gamma_5 \bar{G}(0,x) gamma_5 G(0,x)

where G(0,x) is the quark propagator and \bar{G} is the antiquark propagator. While it is straightforward to optimise the pion to be just the norm of the propagator, lets work with this form, it leaves room for more generalisation.

Let us form the antiquark propagator:

\bar{G}(0,x) = gamma_5 G(0,x) gamma_5

we can do this in chroma as follows:

LatticePropagator anti_quark = Gamma(15)*quark_propagator*Gamma(15);

remembering the spin conventions from tutorial 1 (15 is all four bits set so: gamma_1 gamma_2 gamma_3 gamma_4 which is the definition of gamma_5).

Now we need to trace this with the original propagator inserting gamma_5 matrices. We can do this in Chroma with the line:

LatticeComplex traced_props = trace( Gamma(15)*adj(anti_quark)*Gamma(15)*quark_propagator);

Now we need to sum this over spatial sites only, to give us a correlation function. Usually we also need to convolve the correlation function with momentum phase factors. Chroma has special technology to do this an object called SftMom which stands for SlowFourier TransformMomenta. This object can precompute the momentum phase factors, and use them to Fourier transform correlation functions like traced_props. Let us create such an object with the phases for ZERO momentum. In Chroma this is accomplished by

SftMom phases(0, true, Nd-1);

This creates the momentum phases, with a maximum mom^2 of zero. The true argument asks for equivalent momenta to be averaged, which is not relevant in this case. The Nd-1 nominates direction Nd-1=3 as the time direction (remember these are indexed from 0).

Let us now fourier transform our propagator and sum it over space:

multi2d<DComplex> hsum; // A 2D array of double prec complex numbers // One dimension is indexed by momenta // The other is the timeslice hsum = phases.sft(traced_props); // Apply the fourier transform

The array hsum is indexed as:

hsum[ momentum_index ][ t ];

with the number of momenta determined from the SftMom class by its member function numMom(). Further we can convert the momentum index to the actual 3-vector momentum, using the member function numToMom(). Let us use these to write out our correlator into an XML file. We will use an XMLArrayWriter to write this out.

// Create an XMLArrayWriter to write out am array of numMom() elements XMLArrayWriter momenta(xml_out, phases.numMom()); push(momenta, "PseudoScalar"); // Array will be called PseudoScalar // loop over all the momenta for(int i=0; i < phases.numMom(); i++) { // Special XMLArrayWriter instruction to start a new array element push(momenta); // Write out the momentum index write(momenta, "mom_index", i); // Write out the 3 momentum corresponding to this index write(momenta, "mom", phases.numToMom(i)); // Write correlator write(momenta, "correlator", hsum[i]); // We are done with this element pop(momenta); } pop(momenta); // Pop the toplevel tag

If you have been conciensously adding all this code as we went along you should now save your program, and run make to remake it. Now run it with (with input and output XML files) and look at the output XML file. It should have at its end the following snippet:

<PseudoScalar> <elem> <mom_index>0</mom_index> <mom>0 0 0</mom> <correlator> <elem> <re>0.147564</re> <im>0</im> </elem> <elem> <re>0.00923146</re> <im>0</im> </elem> <elem> <re>0.0028872</re> <im>0</im> </elem> <elem> <re>0.00109399</re> <im>0</im> </elem> <elem> <re>0.000683482</re> <im>0</im> </elem> <elem> <re>0.00110199</re> <im>0</im> </elem> <elem> <re>0.00286503</re> <im>0</im> </elem> <elem> <re>0.00919643</re> <im>0</im> </elem> </correlator> </elem> </PseudoScalar>

Exercise 5: 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.

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 6: 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 7: 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.

Phew !!! You have completed Tutorial 3

You should now be able to read a lot of the measurement code, understand basic XML manipulation, write your own readers/writers, measurements and link against an installed library