This is my first approach to an openmp version of the code. Still has to be crosschecked against purely serial version.

qed/2MN_integrator.c

@@ -92,10 +92,15 @@ void update_momenta_fermion(const double dtau) {
   double f1=0., f2=0., sqsum = 0.;
   g_cgiterations1 += cg(g_X, g_fermion, ITER_MAX, DELTACG, &gam5D_SQR_wilson);
   gam5D_wilson(g_gam5DX, g_X);
+#ifdef OMP
+#pragma omp parallel for private(f1,f2)
+#endif
   for(i = 0; i < GRIDPOINTS; i++) {
     f1 = trX_dQ_wilson_dalpha1_X(i);
     f2 = trX_dQ_wilson_dalpha2_X(i);
+#ifdef _DEBUG_
     sqsum = f1*f1 + f2*f2;
+#endif
     gp1[i] = gp1[i] - dtau*(- f1);
     gp2[i] = gp2[i] - dtau*(- f2);
   }
@@ -112,10 +117,15 @@ void update_momenta_PF2(const double dtau) {
   double sqsum = 0.;
   g_cgiterations2 += cg(g_X, g_fermion2, ITER_MAX, DELTACG, &gam5D_SQR_wilson);
   gam5D_wilson(g_gam5DX, g_X);
+#ifdef OMP
+#pragma omp parallel for private(f1,f2)
+#endif
   for(i = 0; i < GRIDPOINTS; i++) {
     f1 = g_musqr*trX_dQ_wilson_dalpha1_X(i);
     f2 = g_musqr*trX_dQ_wilson_dalpha2_X(i);
+#ifdef _DEBUG_
     sqsum = f1*f1 + f2*f2;
+#endif
     gp1[i] = gp1[i] - dtau*(- f1);
     gp2[i] = gp2[i] - dtau*(- f2);
   }
@@ -133,10 +143,15 @@ void update_momenta_PF1(const double dtau) {
   double sqsum = 0.;
   g_cgiterations1 += cg(g_X, g_fermion, ITER_MAX, DELTACG, &gam5D_SQR_musqr_wilson);
   gam5D_wilson(g_gam5DX, g_X);
+#ifdef OMP
+#pragma omp parallel for private(f1,f2)
+#endif
   for(i = 0; i < GRIDPOINTS; i++) {
     f1 = trX_dQ_wilson_dalpha1_X(i);
     f2 = trX_dQ_wilson_dalpha2_X(i);
+#ifdef _DEBUG_
     sqsum = f1*f1 + f2*f2;
+#endif
     gp1[i] = gp1[i] - dtau*(- f1);
     gp2[i] = gp2[i] - dtau*(- f2);
   }
@@ -151,6 +166,9 @@ void update_momenta_PF1(const double dtau) {
 void update_momenta_gauge(const double dtau) {
   int i;
   double f1=0., f2=0., sqsum = 0.;
+#ifdef OMP
+#pragma omp parallel for private (f1,f2)
+#endif
   for(i = 0; i < GRIDPOINTS; i++) {
     f1 = DS_G1(i);
     f2 = DS_G2(i);

qed/CMakeLists.txt

@@ -0,0 +1,36 @@
+cmake_minimum_required(VERSION 2.4)
+
+# The name of our project is "HELLO".  CMakeLists files in this project can
+# refer to the root source directory of the project as ${HELLO_SOURCE_DIR} and
+# to the root binary directory of the project as ${HELLO_BINARY_DIR}.
+project (maxmize_tr)
+
+
+# intel compiling
+# SET ( CMAKE_C_COMPILER "icc" )
+
+SET ( CMAKE_C_FLAGS "-std=c99 -pg -fopenmp -DOMP" )
+## SET ( CMAKE_C_FLAGS "-std=c99 -pg" )
+
+INCLUDE_DIRECTORIES(./)
+
+add_subdirectory( rand )
+
+
+
+
+
+SET (BASE_SRC rand/gauss.c  rand/ranlxd.c  rand/ranlxs.c 2MN_integrator.c  fields.c   leapfrog.c   measurements.c  statistics.c  hmc.c   leapfrog2.c  dirac.c  lattice.c  linalg.c  rec_lf_integrator.c )
+
+
+
+ADD_EXECUTABLE(qed qed.c ${BASE_SRC})
+
+
+set( LIBS ${LIBS} "-lm" )
+
+target_link_libraries( qed ${LIBS} )
+
+
+
+

qed/dirac.c

@@ -39,6 +39,10 @@
 void gam5D_wilson(spinor *out, spinor *in) {
   int i;
   double factor = (2*g_R + g_mass);
+
+#ifdef OMP
+#pragma omp parallel for 
+#endif
   for(i=0; i<GRIDPOINTS; i++) {
     complex double link1_i = link1[i];
     spinor in_i = in[i];

qed/leapfrog.c

@@ -29,6 +29,9 @@ void update_momenta(const double dtau)
   int i;
   g_cgiterations1 += cg(g_X, g_fermion, ITER_MAX, DELTACG, &gam5D_SQR_wilson);
   gam5D_wilson(g_gam5DX, g_X);
+#ifdef OMP
+#pragma omp parallel for
+#endif
   for(i = 0; i < GRIDPOINTS; i++) {
     gp1[i] = gp1[i] - dtau*(DS_G1(i) - trX_dQ_wilson_dalpha1_X(i));
     gp2[i] = gp2[i] - dtau*(DS_G2(i) - trX_dQ_wilson_dalpha2_X(i));
@@ -38,6 +41,9 @@ void update_momenta(const double dtau)
 
 void update_gauge(const double dtau) {
   int i;
+#ifdef OMP
+#pragma omp parallel for
+#endif
   for(i = 0; i < GRIDPOINTS; i++) {
     gauge1[i] = gauge1[i] + dtau*gp1[i];
     gauge2[i] = gauge2[i] + dtau*gp2[i];

qed/linalg.c

@@ -10,6 +10,9 @@ void add(spinor *Q, spinor *R, spinor *S)
   int ix;
   spinor *q,*r,*s;
 
+#ifdef OMP
+#pragma omp parallel for private(q,r,s)
+#endif
   for (ix=0;ix<GRIDPOINTS;ix++)
   {
     q = (spinor *) Q + ix;
@@ -43,6 +46,9 @@ void assign(spinor *R, spinor *S)
   int ix;
   spinor *r,*s;
 
+#ifdef OMP
+#pragma omp parallel for private(r,s)
+#endif
   for (ix=0;ix<GRIDPOINTS;ix++){
     r=(spinor *) R + ix;
     s=(spinor *) S + ix;
@@ -59,6 +65,9 @@ void assign_add_mul(spinor *P, spinor *Q, complex double c)
   spinor *r,*s;
 
 
+#ifdef OMP
+#pragma omp parallel for private(r,s)
+#endif
   for (ix=0;ix<GRIDPOINTS;ix++){
     r=(spinor *) P + ix;
     s=(spinor *) Q + ix;
@@ -74,7 +83,9 @@ void assign_add_mul_r(spinor *P, spinor *Q, double c)
   static double fact;
 
   fact=c;
-
+#ifdef OMP
+#pragma omp parallel for
+#endif
   for (ix=0;ix<GRIDPOINTS;ix++){
 
     P[ix].s1+=fact*Q[ix].s1;
@@ -86,6 +97,9 @@ void assign_diff_mul(spinor *R, spinor *S, complex double c){
   int ix;
   spinor *r, *s;
 
+#ifdef OMP
+#pragma omp parallel for private(r,s)
+#endif
   for (ix=0;ix<GRIDPOINTS;ix++)
   {
     s = (spinor *) S + ix;
@@ -104,6 +118,9 @@ void assign_mul_add_r(spinor *R, spinor *S, double c)
 
   fact=c;
 
+#ifdef OMP
+#pragma omp parallel for private(r,s)
+#endif
   for (ix=0;ix<GRIDPOINTS;ix++){
     r=(spinor *) R + ix;
     s=(spinor *) S + ix;
@@ -118,6 +135,9 @@ void diff(spinor *Q, spinor *R, spinor *S)
   int ix;
   spinor *q,*r,*s;
 
+#ifdef OMP
+#pragma omp parallel for private(q,r,s)
+#endif
   for (ix=0;ix<GRIDPOINTS;ix++){
     q=(spinor *) Q + ix;
     r=(spinor *) R + ix;
@@ -133,6 +153,9 @@ void mul_r(spinor *R, double c, spinor *S)
   int ix;
   spinor *r,*s;
 
+#ifdef OMP
+#pragma omp parallel for private(r,s)
+#endif
   for (ix=0;ix<GRIDPOINTS;ix++){
     r=(spinor *) R + ix;
     s=(spinor *) S + ix;
@@ -147,6 +170,9 @@ void mul_c(spinor *R, complex double c, spinor *S)
   int ix;
   spinor *r,*s;
 
+#ifdef OMP
+#pragma omp parallel for private(r,s)
+#endif
   for (ix=0;ix<GRIDPOINTS;ix++){
     r=(spinor *) R + ix;
     s=(spinor *) S + ix;
@@ -166,6 +192,9 @@ complex double scalar_prod(spinor *S, spinor *R)
 
   ds=0.0 + I*0.0;
 
+#ifdef OMP
+#pragma omp parallel for private(r,s)
+#endif
   for (ix=0;ix<GRIDPOINTS;ix++){
     s=(spinor *) S + ix;
     r=(spinor *) R + ix;

qed/qed.c

@@ -6,6 +6,10 @@
 #include <string.h>
 #include <unistd.h>
 
+#ifdef OMP
+#include <omp.h>
+#endif
+
 #include "statistics.h"
 #include "lattice.h"
 #include "hmc.h"
@@ -27,6 +31,8 @@ double beta = 1.0;       //Coupling constant for the gauge field
 
 double thermalize_min_acc = 0.7; // minimum ratio of accepted thermalization steps
 
+int g_omp_nthreads = 1;
+
 void echo_sim_params();
 void save_gauge(const char *filename);
 void load_gauge(const char *filename);
@@ -62,6 +68,7 @@ int main(int argc, char **argv)
     {"tau", required_argument, NULL, 0},
     {"thermalize_min_acc", required_argument, NULL, 0},
     {"no_timescales", required_argument, NULL, 0},
+    {"omp_nthreads", required_argument, NULL, 0},
     {0, 0, 0, 0}
   };
 
@@ -80,6 +87,7 @@ int main(int argc, char **argv)
 
     const char *optionName = long_options[option_index].name;
     double optionDoubleValue = strtod(optarg, NULL);
+    int optionIntValue = atoi( optarg );
     printf("%s = %f\n", optionName, optionDoubleValue);
 
     if (strcmp(optionName, "thermalize") == 0)
@@ -106,6 +114,14 @@ int main(int argc, char **argv)
       thermalize_min_acc = optionDoubleValue;
     else if (strcmp(optionName, "no_timescales") == 0)
       no_timescales = optionDoubleValue;
+    else if (strcmp(optionName, "omp_nthreads") == 0){
+#ifdef OMP
+      g_omp_nthreads = optionIntValue;
+#else
+      fprintf(stderr," !!!  Warning code compiled without openMP capability.\
+ Option \"--omp_nthreads\" without effect. !!! \n");
+#endif
+    }
   }
 
   // setup integration parameters
@@ -148,6 +164,11 @@ int main(int argc, char **argv)
   coldstart();
   /* Print out the run parameters */
   echo_sim_params(); 
+
+
+#ifdef OMP
+  omp_set_num_threads(g_omp_nthreads);
+#endif
 
   /* thermalization */
   printf("\n Thermalization: \n\n");
@@ -349,7 +370,8 @@ void echo_sim_params()
 void save_gauge(const char *filename)
 {
   FILE *file;
-
+
+  fprintf(  stderr , "printing to file \"%s\" \n" , filename );
   file = fopen(filename, "w");
   fprintf(file, "%.16lg\n\n", s_g_old);
   for (int i = 0; i < GRIDPOINTS; i ++)

qed/rand/CMakeLists.txt

@@ -0,0 +1,3 @@
+
+
+set( RAND_SRC gauss.c ranlxd.c ranlxs.c )