From fd70a495b78f56181d6a16d61788268c2c91daf0 Mon Sep 17 00:00:00 2001
From: Andreas Nube <nubi@freenet.de>
Date: Thu, 3 May 2012 17:13:17 +0200
Subject: [PATCH 1/2] added sopport for openmp

---
 qed/2MN_integrator.c    | 10 ++++++++++
 qed/CMakeLists.txt      | 36 ++++++++++++++++++++++++++++++++++++
 qed/dirac.c             |  2 ++
 qed/leapfrog.c          |  2 ++
 qed/linalg.c            | 11 ++++++++++-
 qed/qed.c               | 15 ++++++++++++++-
 qed/rand/CMakeLists.txt |  3 +++
 7 files changed, 77 insertions(+), 2 deletions(-)
 create mode 100644 qed/CMakeLists.txt
 create mode 100644 qed/rand/CMakeLists.txt

diff --git a/qed/2MN_integrator.c b/qed/2MN_integrator.c
index 420c0ea..180c018 100644
--- a/qed/2MN_integrator.c
+++ b/qed/2MN_integrator.c
@@ -92,10 +92,13 @@ 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);
+#pragma omp parallel for private(f1,f2)
   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 +115,13 @@ 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);
+#pragma omp parallel for private(f1,f2)
   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 +139,13 @@ 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);
+#pragma omp parallel for private(f1,f2)
   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 +160,7 @@ void update_momenta_PF1(const double dtau) {
 void update_momenta_gauge(const double dtau) {
   int i;
   double f1=0., f2=0., sqsum = 0.;
+#pragma omp parallel for private (f1,f2)
   for(i = 0; i < GRIDPOINTS; i++) {
     f1 = DS_G1(i);
     f2 = DS_G2(i);
diff --git a/qed/CMakeLists.txt b/qed/CMakeLists.txt
new file mode 100644
index 0000000..9ce253f
--- /dev/null
+++ b/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} )
+
+
+
+
diff --git a/qed/dirac.c b/qed/dirac.c
index 4aa9f82..3916f2b 100644
--- a/qed/dirac.c
+++ b/qed/dirac.c
@@ -39,6 +39,8 @@
 void gam5D_wilson(spinor *out, spinor *in) {
   int i;
   double factor = (2*g_R + g_mass);
+
+#pragma omp parallel for 
   for(i=0; i<GRIDPOINTS; i++) {
     complex double link1_i = link1[i];
     spinor in_i = in[i];
diff --git a/qed/leapfrog.c b/qed/leapfrog.c
index 4600670..f7ce30a 100644
--- a/qed/leapfrog.c
+++ b/qed/leapfrog.c
@@ -29,6 +29,7 @@ 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);
+#pragma omp parallel for
   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 +39,7 @@ void update_momenta(const double dtau)
 
 void update_gauge(const double dtau) {
   int i;
+#pragma omp parallel for
   for(i = 0; i < GRIDPOINTS; i++) {
     gauge1[i] = gauge1[i] + dtau*gp1[i];
     gauge2[i] = gauge2[i] + dtau*gp2[i];
diff --git a/qed/linalg.c b/qed/linalg.c
index 4e2bb04..65bb21b 100644
--- a/qed/linalg.c
+++ b/qed/linalg.c
@@ -10,6 +10,7 @@ void add(spinor *Q, spinor *R, spinor *S)
   int ix;
   spinor *q,*r,*s;
   
+#pragma omp parallel for private(q,r,s)
   for (ix=0;ix<GRIDPOINTS;ix++)
   {
     q = (spinor *) Q + ix;
@@ -43,6 +44,7 @@ void assign(spinor *R, spinor *S)
   int ix;
   spinor *r,*s;
   
+#pragma omp parallel for private(r,s)
   for (ix=0;ix<GRIDPOINTS;ix++){
     r=(spinor *) R + ix;
     s=(spinor *) S + ix;
@@ -59,6 +61,7 @@ void assign_add_mul(spinor *P, spinor *Q, complex double c)
   spinor *r,*s;
 
   
+#pragma omp parallel for private(r,s)
   for (ix=0;ix<GRIDPOINTS;ix++){
     r=(spinor *) P + ix;
     s=(spinor *) Q + ix;
@@ -74,7 +77,7 @@ void assign_add_mul_r(spinor *P, spinor *Q, double c)
   static double fact;
 
   fact=c;
-   
+#pragma omp parallel for
   for (ix=0;ix<GRIDPOINTS;ix++){
 
     P[ix].s1+=fact*Q[ix].s1;
@@ -86,6 +89,7 @@ void assign_diff_mul(spinor *R, spinor *S, complex double c){
   int ix;
   spinor *r, *s;
 
+#pragma omp parallel for private(r,s)
   for (ix=0;ix<GRIDPOINTS;ix++)
   {
     s = (spinor *) S + ix;
@@ -104,6 +108,7 @@ void assign_mul_add_r(spinor *R, spinor *S, double c)
   
   fact=c;
   
+#pragma omp parallel for private(r,s)
   for (ix=0;ix<GRIDPOINTS;ix++){
     r=(spinor *) R + ix;
     s=(spinor *) S + ix;
@@ -118,6 +123,7 @@ void diff(spinor *Q, spinor *R, spinor *S)
   int ix;
   spinor *q,*r,*s;
   
+#pragma omp parallel for private(q,r,s)
   for (ix=0;ix<GRIDPOINTS;ix++){
     q=(spinor *) Q + ix;
     r=(spinor *) R + ix;
@@ -133,6 +139,7 @@ void mul_r(spinor *R, double c, spinor *S)
   int ix;
   spinor *r,*s;
 
+#pragma omp parallel for private(r,s)
   for (ix=0;ix<GRIDPOINTS;ix++){
     r=(spinor *) R + ix;
     s=(spinor *) S + ix;
@@ -147,6 +154,7 @@ void mul_c(spinor *R, complex double c, spinor *S)
   int ix;
   spinor *r,*s;
 
+#pragma omp parallel for private(r,s)
   for (ix=0;ix<GRIDPOINTS;ix++){
     r=(spinor *) R + ix;
     s=(spinor *) S + ix;
@@ -166,6 +174,7 @@ complex double scalar_prod(spinor *S, spinor *R)
 
   ds=0.0 + I*0.0;
   
+#pragma omp parallel for private(r,s)
   for (ix=0;ix<GRIDPOINTS;ix++){
     s=(spinor *) S + ix;
     r=(spinor *) R + ix;
diff --git a/qed/qed.c b/qed/qed.c
index 863d62b..d3434a6 100644
--- a/qed/qed.c
+++ b/qed/qed.c
@@ -5,6 +5,7 @@
 #include <getopt.h>
 #include <string.h>
 #include <unistd.h>
+#include <omp.h>
 
 #include "statistics.h"
 #include "lattice.h"
@@ -27,6 +28,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 +65,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 +84,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 +111,8 @@ 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)
+      g_omp_nthreads = optionIntValue;
   }
   
   // setup integration parameters
@@ -148,6 +155,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 +361,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 ++)
diff --git a/qed/rand/CMakeLists.txt b/qed/rand/CMakeLists.txt
new file mode 100644
index 0000000..c302a8e
--- /dev/null
+++ b/qed/rand/CMakeLists.txt
@@ -0,0 +1,3 @@
+
+
+set( RAND_SRC gauss.c ranlxd.c ranlxs.c )

From 500df254a1b1a917a8dee8241417ee911b221ad2 Mon Sep 17 00:00:00 2001
From: Andreas Nube <annube@ifh.de>
Date: Mon, 7 May 2012 13:28:59 +0200
Subject: [PATCH 2/2] The "#pramga omp ..." statements are now compiled only if
 "OMP" macro is defined. A warning will be displayed if OMP was not defined
 but the a value was given with the --omp_nthreads command line option.

---
 qed/2MN_integrator.c |  8 ++++++++
 qed/dirac.c          |  2 ++
 qed/leapfrog.c       |  4 ++++
 qed/linalg.c         | 20 ++++++++++++++++++++
 qed/qed.c            | 11 ++++++++++-
 5 files changed, 44 insertions(+), 1 deletion(-)

diff --git a/qed/2MN_integrator.c b/qed/2MN_integrator.c
index 180c018..f7e5afe 100644
--- a/qed/2MN_integrator.c
+++ b/qed/2MN_integrator.c
@@ -92,7 +92,9 @@ 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);
@@ -115,7 +117,9 @@ 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);
@@ -139,7 +143,9 @@ 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);
@@ -160,7 +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);
diff --git a/qed/dirac.c b/qed/dirac.c
index 3916f2b..3e52115 100644
--- a/qed/dirac.c
+++ b/qed/dirac.c
@@ -40,7 +40,9 @@ 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];
diff --git a/qed/leapfrog.c b/qed/leapfrog.c
index f7ce30a..63031e4 100644
--- a/qed/leapfrog.c
+++ b/qed/leapfrog.c
@@ -29,7 +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));
@@ -39,7 +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];
diff --git a/qed/linalg.c b/qed/linalg.c
index 65bb21b..6280f6d 100644
--- a/qed/linalg.c
+++ b/qed/linalg.c
@@ -10,7 +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;
@@ -44,7 +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;
@@ -61,7 +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;
@@ -77,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;
@@ -89,7 +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;
@@ -108,7 +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;
@@ -123,7 +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;
@@ -139,7 +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;
@@ -154,7 +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;
@@ -174,7 +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;
diff --git a/qed/qed.c b/qed/qed.c
index d3434a6..649e295 100644
--- a/qed/qed.c
+++ b/qed/qed.c
@@ -5,7 +5,10 @@
 #include <getopt.h>
 #include <string.h>
 #include <unistd.h>
+
+#ifdef OMP
 #include <omp.h>
+#endif
 
 #include "statistics.h"
 #include "lattice.h"
@@ -111,8 +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)
+    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