Exercise 5: Writing a gauge field
Goals: Write files using QDP I/O
Concepts: I/O write, Subsets
Directory: $HOME/examples/05-writing-gauge
How to: $ cd $HOME/examples/05-writing-gauge
$ make
$ submit plaquette /scratch/data/propagator_8x8x8x16 new-gauge

Now let us write some new gauge field files. QDP provides different modes for writing data files. For this reason QDP_open_write() takes an extra argument, volfmt describing what kind of file is to be written. On Blue Gene we can safely use QDP_SINGLEFILE. Presently there are two file formats supported by QDP:

/* volume formats */
#define QDP_SINGLEFILE QIO_SINGLEFILE 
#define QDP_MULTIFILE QIO_MULTIFILE

The begining of write-gauge.c is similar to read-gauge.c from the previous exercise. We start with associating a writer with a file name on line 10, then write an array of four lattice matrices into it on line 16 and close the writer at line 17. 

  1   #include <qdp-config.h>
  2
  3   int
  4   write_gauge(char *name, QDP_ColorMatrix *link[])
  5   {
  6     QDP_String *id = QDP_string_create();
  7     QDP_Writer *writer;
  8
  9     QDP_string_set(id, "write_gauge file metadata");
 10     writer = QDP_open_write(id, name, QDP_SINGLEFILE);
 11     if (writer == 0) {
 12       printf0("open_write(%s) failed\n", name);
 13       return 1;
 14     }
 15     QDP_string_set(id, "write_gauge record metadata");
 16     QDP_vwrite_M(writer, id, link, 4);
 17     QDP_close_write(writer);
 18     QDP_string_destroy(id);
 19     return 0; 
 20   }

In main() we start with reading the gauge on line 40, computing our favorite plaquette on line 46 and printing it at line 49. 

 39   /* read gauge field */
 40   if (read_gauge(U, argv[1]) != 0) {
 41     printf0("ERROR: read_gauge(%s)\n", argv[1]);
 42     goto end;
 43   }
 44
 45   /* Compute plaquette */
 46   plaq = plaquette(U);
 47
 48   /* Display the value */
 49   printf0("initial plaquette = %g\n",
 50           plaq / (QDP_volume() * QDP_Nc * NDIM * (NDIM - 1) / 2 ));

Our next step is to make a random link on sites where x+y+z+t is even (that is what QDP_even stands for. Yes, there is also QDP_odd. The formerly magical QDP_all stands for all sites on the lattice. You can even define your own set of points on a lattice with QDP_create_subset(), but it is beyond the scope of our exercise). Notice that we are being frugal here: the random link is created only on even sites on line 53. On lines 56 and 57 we multiply U[0] by force and store the result back to U[0]. The new value of the plaquette is computed on line 60. 

 35     /* allocate and initialize random state */
 36     rs = QDP_create_S();
 37     init_random(rs, 32582657);
... 
 52     /* random link at even sites */
 53     random_link(force, rs, QDP_even);
 54
 55     /* update U[0] at even sites */
 56     QDP_M_eq_M_times_M(link, force, U[0], QDP_even);
 57     QDP_M_eq_M(U[0], link, QDP_even);
 58
 59     /* compute and print the new plaquette */
 60     plaq = plaquette(U);
 61     printf0("new plaquette = %g\n",
 62     plaq / (QDP_volume() * QDP_Nc * NDIM * (NDIM - 1) / 2));

Freshly updated gauge field is stored into file argv[2] on line 65. 

 64     /* write U to a new file */
 65     if (write_gauge(argv[2], U) != 0) {
 66       printf("ERROR: write_gauge(%s)\n", argv[2]);
 67       goto end;
 68     }

Problems