[tools] g_sans - add OpenMP based parallelism

This commit adds OpenMP based parallelism for
calc_radial_distribution_histogram function

Also it adds check for mcover parameter which should be -1 or (0,1]

Change-Id: Ic14267782c69c25364b69a91719cb451fd8e3a89
Signed-off-by: Alexey Shvetsov <alexxy@omrb.pnpi.spb.ru>

diff --git a/src/tools/gmx_sans.c b/src/tools/gmx_sans.c
index 311fa9c3a4faf43dc943bd57582cc0ba670c834d..ab3730a36b7a506176dcc11a2b072e816c0c9643 100644 (file)
--- a/src/tools/gmx_sans.c
+++ b/src/tools/gmx_sans.c
@@ -55,6 +55,10 @@
 #include "gmx_ana.h"
 #include "nsfactor.h"
 
+#ifdef GMX_OPENMP
+#include <omp.h>
+#endif
+
 
 int gmx_sans(int argc,char *argv[])
 {
@@ -80,6 +84,7 @@ int gmx_sans(int argc,char *argv[])
     static real start_q=0.0, end_q=2.0, q_step=0.01;
     static real mcover=-1;
     static unsigned int  seed=0;
+    static int           nthreads=-1;
 
     static const char *emode[]= { NULL, "direct", "mc", NULL };
     static const char *emethod[]={ NULL, "debye", "fft", NULL };
@@ -95,7 +100,7 @@ int gmx_sans(int argc,char *argv[])
         { "-mode", FALSE, etENUM, {emode},
           "Mode for sans spectra calculation" },
         { "-mcover", FALSE, etREAL, {&mcover},
-          "Monte-Carlo coverage"},
+          "Monte-Carlo coverage should be -1(default) or (0,1]"},
         { "-method", FALSE, etENUM, {emethod},
           "[HIDDEN]Method for sans spectra calculation" },
         { "-pbc", FALSE, etBOOL, {&bPBC},
@@ -110,6 +115,10 @@ int gmx_sans(int argc,char *argv[])
           "Stepping in q (1/nm)"},
         { "-seed",     FALSE, etINT,  {&seed},
           "Random seed for Monte-Carlo"},
+#ifdef GMX_OPENMP
+        { "-nt",  FALSE, etINT, {&nthreads},
+          "Number of threads to start"},
+#endif
     };
   FILE      *fp;
   const char *fnTPX,*fnNDX,*fnDAT=NULL;
@@ -144,13 +153,22 @@ int gmx_sans(int argc,char *argv[])
       { efXVG, "-pr",         "pr",   ffWRITE }
   };
 
+#ifdef GMX_OPENMP
+    nthreads = omp_get_max_threads();
+#endif
+
   CopyRight(stderr,argv[0]);
   parse_common_args(&argc,argv,PCA_BE_NICE,
                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 
   /* check that binwidth not smaller than smallers distance */
   check_binwidth(binwidth);
+  check_mcover(mcover);
 
+  /* setting number of omp threads globaly */
+#ifdef GMX_OPENMP
+  omp_set_num_threads(nthreads);
+#endif
   /* Now try to parse opts for modes */
   switch(emethod[0][0]) {
   case 'd':

diff --git a/src/tools/nsfactor.c b/src/tools/nsfactor.c
index 79ebdbcebcdbfe37a066ec07cb502c9b2e5076af..e7879e7ca0bf613b5e3b4d7a4fb90c75f9b0e235 100644 (file)
--- a/src/tools/nsfactor.c
+++ b/src/tools/nsfactor.c
@@ -46,12 +46,26 @@
 #include "vec.h"
 #include "nsfactor.h"
 
+#ifdef GMX_OPENMP
+#include <omp.h>
+#endif
+
 void check_binwidth(real binwidth) {
     real smallest_bin=0.1;
     if (binwidth<smallest_bin)
         gmx_fatal(FARGS,"Binwidth shouldnt be smaller then smallest bond length (H-H bond ~0.1nm) in a box");
 }
 
+void check_mcover(real mcover) {
+    if (mcover>1.0) {
+        gmx_fatal(FARGS,"mcover should be -1 or (0,1]");
+    } else if ((mcover<0)&(mcover != -1)) {
+        gmx_fatal(FARGS,"mcover should be -1 or (0,1]");
+    } else {
+        return;
+    }
+}
+
 void normalize_probability(int n,double *a){
     int i;
     double norm=0.0;
@@ -148,11 +162,17 @@ gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram (
                             real mcover,
                             unsigned int seed) {
     gmx_radial_distribution_histogram_t    *pr=NULL;
-    rvec        dist;
-    double      rmax;
-    int         i,j;
-    gmx_large_int_t   mc=0,max;
-    gmx_rng_t   rng=NULL;
+    rvec            dist;
+    double          rmax;
+    int             i,j;
+#ifdef GMX_OPENMP
+    double          **tgr;
+    int             tid;
+    int             nthreads;
+    gmx_rng_t       *trng=NULL;
+#endif
+    gmx_large_int_t mc=0,max;
+    gmx_rng_t       rng=NULL;
 
     /* allocate memory for pr */
     snew(pr,1);
@@ -181,18 +201,86 @@ gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram (
             max=(gmx_large_int_t)floor(0.5*mcover*isize*(isize-1));
         }
         rng=gmx_rng_init(seed);
-        while(mc<max) {
+#ifdef GMX_OPENMP
+        nthreads = omp_get_max_threads();
+        snew(tgr,nthreads);
+        snew(trng,nthreads);
+        for(i=0;i<nthreads;i++){
+            snew(tgr[i],pr->grn);
+            trng[i]=gmx_rng_init(gmx_rng_uniform_uint32(rng));
+        }
+        #pragma omp parallel shared(tgr,trng,mc) private(tid,i,j)
+        {
+            tid = omp_get_thread_num();
+            /* now starting parallel threads */
+            #pragma omp for
+            for(mc=0;mc<max;mc++) {
+                i=(int)floor(gmx_rng_uniform_real(trng[tid])*isize);
+                j=(int)floor(gmx_rng_uniform_real(trng[tid])*isize);
+                if(i!=j) {
+                    tgr[tid][(int)floor(sqrt(distance2(x[index[i]],x[index[j]]))/binwidth)]+=gsans->slength[index[i]]*gsans->slength[index[j]];
+                }
+            }
+        }
+        /* collecting data from threads */
+        for(i=0;i<pr->grn;i++) {
+            for(j=0;j<nthreads;j++) {
+                pr->gr[i] += tgr[j][i];
+            }
+        }
+        /* freeing memory for tgr and destroying trng */
+        for(i=0;i<nthreads;i++) {
+            sfree(tgr[i]);
+            gmx_rng_destroy(trng[i]);
+        }
+        sfree(tgr);
+        sfree(trng);
+#else
+        for(mc=0;mc<max;mc++) {
             i=(int)floor(gmx_rng_uniform_real(rng)*isize);
             j=(int)floor(gmx_rng_uniform_real(rng)*isize);
             if(i!=j)
                 pr->gr[(int)floor(sqrt(distance2(x[index[i]],x[index[j]]))/binwidth)]+=gsans->slength[index[i]]*gsans->slength[index[j]];
-            mc++;
         }
+#endif
         gmx_rng_destroy(rng);
     } else {
-        for(i=0;i<isize;i++)
-            for(j=0;j<i;j++)
+#ifdef GMX_OPENMP
+        nthreads = omp_get_max_threads();
+        /* Allocating memory for tgr arrays */
+        snew(tgr,nthreads);
+        for(i=0;i<nthreads;i++) {
+            snew(tgr[i],pr->grn);
+        }
+        #pragma omp parallel shared(tgr) private(tid,i,j)
+        {
+            tid = omp_get_thread_num();
+            /* starting parallel threads */
+            #pragma omp for
+            for(i=0;i<isize;i++) {
+                for(j=0;j<i;j++) {
+                    tgr[tid][(int)floor(sqrt(distance2(x[index[i]],x[index[j]]))/binwidth)]+=gsans->slength[index[i]]*gsans->slength[index[j]];
+                }
+            }
+        }
+        /* collecating data for pr->gr */
+        for(i=0;i<pr->grn;i++) {
+            for(j=0;j<nthreads;j++) {
+                pr->gr[i] += tgr[j][i];
+            }
+        }
+        /* freeing memory for tgr */
+        for(i=0;i<nthreads;i++) {
+            sfree(tgr[i]);
+        }
+        sfree(tgr);
+#else
+        for(i=0;i<isize;i++) {
+            for(j=0;j<i;j++) {
                 pr->gr[(int)floor(sqrt(distance2(x[index[i]],x[index[j]]))/binwidth)]+=gsans->slength[index[i]]*gsans->slength[index[j]];
+            }
+        }
+#endif
     }
 
     /* normalize */

diff --git a/src/tools/nsfactor.h b/src/tools/nsfactor.h
index 7f2d9cc59e78c051f3d7a5559bcbee5e8dbc1ef1..3e6e5fa562cab728d3732e15b23f7403ea420a93 100644 (file)
--- a/src/tools/nsfactor.h
+++ b/src/tools/nsfactor.h
@@ -74,13 +74,15 @@ typedef struct gmx_static_structurefator_t {
 
 void check_binwidth(real binwidth);
 
+void check_mcover(real mcover);
+
 void normalize_probability(int n, double *a);
 
 gmx_nentron_atomic_structurefactors_t *gmx_neutronstructurefactors_init(const char *datfn);
 
 gmx_sans_t *gmx_sans_init(t_topology *top, gmx_nentron_atomic_structurefactors_t *gnsf);
 
-gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram  ( gmx_sans_t *gsans,
+gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram  (gmx_sans_t *gsans,
                             rvec *x,
                             matrix box,
                             atom_id *index,