From afe2a11219723242db246c5dca7aea4e4739f2aa Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Sun, 2 Oct 2011 14:46:12 -0400
Subject: [PATCH] multi-thread one more neighbor list build type

---
 src/neigh_half_nsq.cpp | 89 ++++++++++++++++++++++++++----------------
 1 file changed, 55 insertions(+), 34 deletions(-)

diff --git a/src/neigh_half_nsq.cpp b/src/neigh_half_nsq.cpp
index 32792a678f..5ad36f756f 100644
--- a/src/neigh_half_nsq.cpp
+++ b/src/neigh_half_nsq.cpp
@@ -14,9 +14,14 @@
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "atom.h"
+#include "comm.h"
 #include "group.h"
 #include "error.h"
 
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
+
 using namespace LAMMPS_NS;
 
 /* ----------------------------------------------------------------------
@@ -27,7 +32,17 @@ using namespace LAMMPS_NS;
 
 void Neighbor::half_nsq_no_newton(NeighList *list)
 {
-  int i,j,n,itype,jtype,which,bitmask;
+  const int nthreads = comm->nthreads;
+  const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
+  const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0;
+
+  NEIGH_OMP_INIT;
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(list)
+#endif
+  NEIGH_OMP_SETUP(nlocal);
+
+  int i,j,n,itype,jtype,which;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *neighptr;
 
@@ -41,32 +56,30 @@ void Neighbor::half_nsq_no_newton(NeighList *list)
   int *type = atom->type;
   int *mask = atom->mask;
   int *molecule = atom->molecule;
-  int nlocal = atom->nlocal;
-  int nall = nlocal + atom->nghost;
+  int nall = atom->nlocal + atom->nghost;
   int molecular = atom->molecular;
-  if (includegroup) {
-    nlocal = atom->nfirst;
-    bitmask = group->bitmask[includegroup];
-  }
 
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
-  int **pages = list->pages;
 
-  int inum = 0;
-  int npage = 0;
+  int npage = tid;
   int npnt = 0;
 
-  for (i = 0; i < nlocal; i++) {
+  for (i = ifrom; i < ito; i++) {
 
     if (pgsize - npnt < oneatom) {
       npnt = 0;
-      npage++;
-      if (npage == list->maxpage) pages = list->add_pages();
+      npage += nthreads;
+      // only one thread at a time may check whether we
+      // need new neighbor list pages and then add to them.
+#if defined(_OPENMP)
+#pragma omp critical
+#endif
+      if (npage >= list->maxpage) list->add_pages(nthreads);
     }
 
-    neighptr = &pages[npage][npnt];
+    neighptr = &(list->pages[npage][npnt]);
     n = 0;
 
     itype = type[i];
@@ -94,15 +107,15 @@ void Neighbor::half_nsq_no_newton(NeighList *list)
       }
     }
 
-    ilist[inum++] = i;
+    ilist[i] = i;
     firstneigh[i] = neighptr;
     numneigh[i] = n;
     npnt += n;
     if (n > oneatom || npnt >= pgsize)
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
   }
-
-  list->inum = inum;
+  NEIGH_OMP_CLOSE;
+  list->inum = nlocal;
 }
 
 /* ----------------------------------------------------------------------
@@ -113,7 +126,17 @@ void Neighbor::half_nsq_no_newton(NeighList *list)
 
 void Neighbor::half_nsq_newton(NeighList *list)
 {
-  int i,j,n,itype,jtype,itag,jtag,which,bitmask;
+  const int nthreads = comm->nthreads;
+  const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
+  const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0;
+
+  NEIGH_OMP_INIT;
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(list)
+#endif
+  NEIGH_OMP_SETUP(nlocal);
+
+  int i,j,n,itype,jtype,itag,jtag,which;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *neighptr;
 
@@ -127,32 +150,30 @@ void Neighbor::half_nsq_newton(NeighList *list)
   int *type = atom->type;
   int *mask = atom->mask;
   int *molecule = atom->molecule;
-  int nlocal = atom->nlocal;
-  int nall = nlocal + atom->nghost;
+  int nall = atom->nlocal + atom->nghost;
   int molecular = atom->molecular;
-  if (includegroup) {
-    nlocal = atom->nfirst;
-    bitmask = group->bitmask[includegroup];
-  }
 
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
-  int **pages = list->pages;
 
-  int inum = 0;
-  int npage = 0;
+  int npage = tid;
   int npnt = 0;
 
-  for (i = 0; i < nlocal; i++) {
+  for (i = ifrom; i < ito; i++) {
 
     if (pgsize - npnt < oneatom) {
       npnt = 0;
-      npage++;
-      if (npage == list->maxpage) pages = list->add_pages();
+      npage += nthreads;
+      // only one thread at a time may check whether we
+      // need new neighbor list pages and then add to them.
+#if defined(_OPENMP)
+#pragma omp critical
+#endif
+      if (npage >= list->maxpage) list->add_pages(nthreads);
     }
 
-    neighptr = &pages[npage][npnt];
+    neighptr = &(list->pages[npage][npnt]);
     n = 0;
 
     itag = tag[i];
@@ -198,13 +219,13 @@ void Neighbor::half_nsq_newton(NeighList *list)
       }
     }
 
-    ilist[inum++] = i;
+    ilist[i] = i;
     firstneigh[i] = neighptr;
     numneigh[i] = n;
     npnt += n;
     if (n > oneatom || npnt >= pgsize)
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
   }
-
-  list->inum = inum;
+  NEIGH_OMP_CLOSE;
+  list->inum = nlocal;
 }