trying to implement mpi in zobov

Makefile

@@ -11,7 +11,7 @@ qhull:
 	make -C qhull/src
 
 voboz:
-	make -C src install
+	make -C src all
 
 fastmodules:
 	python2.7 python_tools/setup.py build_ext --inplace

parameters/params.py

@@ -84,9 +84,8 @@
 
 # -- Tessellation options -- #
 do_tessellation = True    # if True, does tessellation; if False, only post-processes a previous run
-# vozisol code performs entire tessellation in one shot: more memory-intensive, but handles survey data better
-use_vozisol = True  # set True for survey data or small simulation boxes
-# if not using vozisol, tessellation code divides data into chunks to save memory and requires following two options
+use_mpi = True
+# if using MPI, the number of CPUs is taken from nthreads above; following two options control how tasks are split
 zobov_box_div = 2   # no. of subdivisions per box side
 zobov_buffer = 0.1  # fraction of box length overlap between subdivisions
 # -------------------------- #

python_tools/zobov.py

@@ -636,7 +636,7 @@ def read_config(self):
         self.num_part_total = parms.num_mocks + parms.num_tracers
         self.num_tracers = parms.num_tracers
 
-    def zobov_wrapper(self, use_vozisol=False, zobov_box_div=2, zobov_buffer=0.1):
+    def zobov_wrapper(self, use_mpi=False, zobov_box_div=2, zobov_buffer=0.1, nthreads=2):
         """Wrapper function to call C-based ZOBOV codes
 
         Arguments:
@@ -648,7 +648,7 @@ def zobov_wrapper(self, use_vozisol=False, zobov_box_div=2, zobov_buffer=0.1):
 
         begin = time.time()
         # ---run the tessellation--- #
-        if use_vozisol:
+        if not use_mpi:
             print("Calling vozisol to do the tessellation...")
             sys.stdout.flush()
             logfolder = self.output_folder + 'log/'
@@ -665,36 +665,34 @@ def zobov_wrapper(self, use_vozisol=False, zobov_box_div=2, zobov_buffer=0.1):
             if not os.access("%s.vol" % self.handle, os.F_OK):
                 sys.exit("Something went wrong with the tessellation. Aborting ...")
         else:
-            print("Calling vozinit, voz1b1 and voztie to do the tessellation...")
+            print("MPI run: calling voz1b1 and voztie to do the tessellation...")
             sys.stdout.flush()
 
             # ---Step 1: call vozinit to write the script used to call voz1b1 and voztie--- #
             logfolder = self.output_folder + 'log/'
             if not os.access(logfolder, os.F_OK):
                 os.makedirs(logfolder)
-            logfile = logfolder + self.handle + '.out'
+            logfile = logfolder + self.handle + '-mpirun.out'
             log = open(logfile, "w")
-            cmd = ["./bin/vozinit", self.posn_file, str(zobov_buffer), str(self.box_length),
-                   str(zobov_box_div), self.handle]
+            cmd = 'parts="%d"\nparts2=`expr $parts \* $parts`\n' % zobov_box_div
+            cmd += 'parts3=`expr $parts \* $parts \* $parts`\npartslist="$(seq 0 $((parts3-1)))"\n'
+            cmd += 'ncpu="%d"\n' % nthreads
+            cmd += 'mpirun -np $ncpu ./bin/mpifor '
+            cmd += '"./bin/voz1b1 %s %f %f %s %d \$((i/$parts2)) \$((i/$parts % $parts)) \$((i % $parts))" ' + \
+                    '$partslist' % (self.posn_file, zobov_buffer, self.box_length, zobov_box_div)
             subprocess.call(cmd, stdout=log, stderr=log)
             log.close()
 
-            # ---Step 2: call this script to do the tessellation--- #
-            voz_script = "scr" + self.handle
-            cmd = ["./%s" % voz_script]
-            log = open(logfile, 'a')
+            log = open(logfile, "a")
+            cmd = ["./bin/voztie", str(zobov_box_div), self.handle]
             subprocess.call(cmd, stdout=log, stderr=log)
             log.close()
 
             # ---Step 3: check the tessellation was successful--- #
             if not os.access("%s.vol" % self.handle, os.F_OK):
                 sys.exit("Something went wrong with the tessellation. Aborting ...")
 
-            # ---Step 4: remove the script file--- #
-            if os.access(voz_script, os.F_OK):
-                os.unlink(voz_script)
-
-            # ---Step 5: copy the .vol files to .trvol--- #
+             # ---Step 5: copy the .vol files to .trvol--- #
             cmd = ["cp", "%s.vol" % self.handle, "%s.trvol" % self.handle]
             subprocess.call(cmd)
 
@@ -706,6 +704,47 @@ def zobov_wrapper(self, use_vozisol=False, zobov_box_div=2, zobov_buffer=0.1):
                 subprocess.call(cmd, stdout=log, stderr=log)
                 log.close()
 
+            # print("Calling vozinit, voz1b1 and voztie to do the tessellation...")
+            # sys.stdout.flush()
+            #
+            # # ---Step 1: call vozinit to write the script used to call voz1b1 and voztie--- #
+            # logfolder = self.output_folder + 'log/'
+            # if not os.access(logfolder, os.F_OK):
+            #     os.makedirs(logfolder)
+            # logfile = logfolder + self.handle + '.out'
+            # log = open(logfile, "w")
+            # cmd = ["./bin/vozinit", self.posn_file, str(zobov_buffer), str(self.box_length),
+            #        str(zobov_box_div), self.handle]
+            # subprocess.call(cmd, stdout=log, stderr=log)
+            # log.close()
+            #
+            # # ---Step 2: call this script to do the tessellation--- #
+            # voz_script = "scr" + self.handle
+            # cmd = ["./%s" % voz_script]
+            # log = open(logfile, 'a')
+            # subprocess.call(cmd, stdout=log, stderr=log)
+            # log.close()
+            #
+            # # ---Step 3: check the tessellation was successful--- #
+            # if not os.access("%s.vol" % self.handle, os.F_OK):
+            #     sys.exit("Something went wrong with the tessellation. Aborting ...")
+            #
+            # # # ---Step 4: remove the script file--- #
+            # if os.access(voz_script, os.F_OK):
+            #     os.unlink(voz_script)
+            #
+            # # ---Step 5: copy the .vol files to .trvol--- #
+            # cmd = ["cp", "%s.vol" % self.handle, "%s.trvol" % self.handle]
+            # subprocess.call(cmd)
+            #
+            # # ---Step 6: if buffer mocks were used, remove them and flag edge galaxies--- #
+            # # (necessary because voz1b1 and voztie do not do this automatically)
+            # if self.num_mocks > 0:
+            #     cmd = ["./bin/checkedges", self.handle, str(self.num_tracers), str(0.9e30)]
+            #     log = open(logfile, 'a')
+            #     subprocess.call(cmd, stdout=log, stderr=log)
+            #     log.close()
+
         print("Tessellation done.\n")
         sys.stdout.flush()
         finish = time.time()

revolver.py

@@ -202,8 +202,8 @@
         # write a config file
         voidcat.write_config()
         # run ZOBOV
-        voidcat.zobov_wrapper(use_vozisol=parms.use_vozisol, zobov_box_div=parms.zobov_box_div,
-                              zobov_buffer=parms.zobov_buffer)
+        voidcat.zobov_wrapper(use_mpi=parms.use_mpi, zobov_box_div=parms.zobov_box_div,
+                              zobov_buffer=parms.zobov_buffer, nthreads=parms.nthreads)
     else:
         # read the config file from a previous run
         voidcat.read_config()

src/Makefile

@@ -1,20 +1,19 @@
 # Compiler choice:
 CC    = gcc
+# CC    = /usr/local/Cellar/clang-omp/2015-04-01/bin/clang-omp
 MPICC = mpicc
+MPICXX = mpicxx
 CFLAGS = -O2 -g -ansi -Wno-unused-result
+MPICFLAGS = -O2 -g -ansi
 
 MLIBS	=   -lm
 
 QLIB = -L../qhull/src/ -lqhull
 QINC = -I../qhull/src/
 
-EXES = vozinit voz1b1 voztie jovoz boz jozov jozovtrvol jozov-grid vozisol voz1b1_mpi checkedges
+EXES = vozinit voz1b1 voztie jovoz boz jozov jozovtrvol jozov-grid vozisol voz1b1_mpi mpifor checkedges
 EXES_NOMPI = vozinit voz1b1 voztie jovoz boz jozov jozovtrvol jozov-grid vozisol checkedges
 
-all: $(EXES)
-
-all_nompi: $(EXES_NOMPI)
-
 jovoz: jovoz.o findrtop.o voz.h
 	$(CC) $(CFLAGS) -o jovoz jovoz.o findrtop.o $(MLIBS)
 
@@ -45,9 +44,20 @@ voztie: voztie.o readfiles.o voz.h
 vozisol: vozisol.o readfiles.o vozutil.o  voz.h
 	$(CC)  -o vozisol $(CFLAGS) vozisol.o readfiles.o vozutil.o -L. $(QLIB) $(MLIBS)
 
+voz1b1_openmp: voz1b1_openmp.o voz1b1.o readfiles.o vozutil.o  voz.h
+	$(CC)  -o voz1b1_openmp $(CFLAGS) voz1b1_openmp.o voz1b1.o readfiles.o vozutil.o -L. $(QLIB) $(MLIBS)
+
 voz1b1_mpi: voz1b1_mpi.o voz1b1.o readfiles.o vozutil.o voz.h
 	$(MPICC) -o voz1b1_mpi.o $(CFLAGS) -c voz1b1_mpi.c $(QINC)
-	$(MPICC) $(CFLAGS) -o voz1b1_mpi voz1b1_mpi.o voz1b1.o readfiles.o vozutil.o -L. $(QLIB) $(MLIBS)
+	$(MPICC) $(MPICFLAGS) -o voz1b1_mpi voz1b1_mpi.o voz1b1.o readfiles.o vozutil.o -L. $(QLIB) $(MLIBS)
+
+mpifor:
+	$(MPICXX) mpifor.cpp -o mpifor
+
+all: $(EXES)
+	mv $(EXES) ../bin
+
+all_nompi: $(EXES_NOMPI)
 
 install: $(EXES_NOMPI)
 	mv $(EXES_NOMPI) ../bin

src/checkedges.c

@@ -87,9 +87,9 @@ int main(int argc, char *argv[]) {
 	    p[i].adj[p[i].nadj_count] = j;
 	    p[j].adj[p[j].nadj_count] = i;
 	    p[i].nadj_count++; p[j].nadj_count++;
-	  } else{
+	  } /*else{
 	    printf("weird#adj,p %d or %d\t",i,j);
-	  }
+	  }*/
 	} else {
 	  printf("adj(%d)=%d>np\n",i,j); FF;
 	}

src/mpifor.cpp

@@ -0,0 +1,53 @@
+#include <cstdlib>
+#include <cstdio>
+#include <glob.h>
+#include <vector>
+#include <signal.h>
+#include <mpi.h>
+
+using namespace std;
+
+int mod(int a, int b) { int c = a%b; return c < 0 ? c+b : c; }
+
+int main(int argc, char ** argv)
+{
+	vector<char*> args;
+	for(char **i = argv+1; *i; i++)
+		if(false);
+		else args.push_back(*i);
+
+	if(args.empty()) {
+		fprintf(stderr,"Syntax: mpifor [options] command [args]\n");
+		return 1;
+	}
+
+	int nprocs, procid;
+
+	MPI_Init (&argc, &argv);
+	MPI_Comm_rank (MPI_COMM_WORLD, &procid);
+	MPI_Comm_size (MPI_COMM_WORLD, &nprocs);
+
+	// And parse and run
+	char* command = args[0];
+	for(int i = 1; i < args.size(); i++)
+	{
+		int k = i-1;
+		glob_t p;
+		p.gl_offs = 1;
+		glob(args[i], GLOB_BRACE | GLOB_TILDE | GLOB_NOMAGIC, 0, &p);
+		for(int j = mod(procid-k,nprocs); j < p.gl_pathc; j += nprocs)
+		{
+			setenv("i", p.gl_pathv[j], 1);
+			//setenv("id", getenv("OMPI_COMM_WORLD_RANK"), 1);
+			//setenv("nproc", getenv("OMPI_COMM_WORLD_SIZE"), 1);
+			int status = system(command);
+			if (WIFSIGNALED(status) &&
+				(WTERMSIG(status) == SIGINT || WTERMSIG(status) == SIGQUIT))
+				return 1;
+		}
+		globfree(&p);
+	}
+
+	MPI_Finalize();
+	return 0;
+}