first tests with mpi4py

bin/beyn.py

@@ -7,6 +7,7 @@ Technische Universitaet Darmstadt.
 See the LICENSE.txt and README.md for more license and copyright information.
 """
 
+from   mpi4py import MPI
 import solver as Solver
 import sys    as sys
 import numpy  as np
@@ -42,7 +43,9 @@ class Beyn:
         self.__lStart   = lStart
         self.__lStep    = lStep
         self.__rankTol  = rankTol
-        self.__verbose  = verbose
+        self.__myProc   = MPI.COMM_WORLD.Get_rank()
+        self.__nProc    = MPI.COMM_WORLD.Get_size()
+        self.__verbose  = verbose and (self.__myProc == 0)
 
     def simple(self):
         """Solves the eigenvalue problem using Beyn's algorithm (simple version)
@@ -137,10 +140,14 @@ class Beyn:
 
         If beQuite is True, this solver becomes quite; it is verbose otherwise
         """
+        if(self.__myProc == 0):
             self.__verbose = not beQuite
 
     def display(self):
         """Prints this Beyn solver parameters"""
+        if(self.__myProc != 0):
+            return
+
         print "Beyn's contour integral method (simple)"
         print "---------------------------------------"
         print " # Nodes used for the trapezoidal rule:"+ " " +str(self.__nNode)
@@ -177,18 +184,38 @@ class Beyn:
         vHat -- the RHS matrix defing Beyn's integrals
         """
 
-        # Initialise integrals
-        A0 = np.matlib.zeros(vHat.shape, dtype=complex)
-        A1 = np.matlib.zeros(vHat.shape, dtype=complex)
+        # Share integration points among processes
+        baseShare     = self.__nNode // self.__nProc      #Base share for proc
+        residualShare = self.__nNode %  self.__nProc      #Residual share
+        allStart      = [0]                               #Inital start
+
+        for i in range(residualShare):                    #Proc with extra share
+            allStart.append(allStart[-1] + baseShare + 1) # -> previous+base+1
+        for i in range(residualShare, self.__nProc):      #Proc with base share
+            allStart.append(allStart[-1] + baseShare)     # -> previous+base+0
+
+        myStart = allStart[self.__myProc]                 #This proc start
+        myStop  = allStart[self.__myProc+1]               #This proc end
+
+        # Initialise integrals (partial results)
+        A0_partial = np.matlib.zeros(vHat.shape, dtype=complex)
+        A1_partial = np.matlib.zeros(vHat.shape, dtype=complex)
 
         # Loop over integation points and integrate
-        for i in range(self.__nNode):
+        for i in range(myStart, myStop):
             t   = 2 * np.pi * i / self.__nNode;
             phi = self.__origin + self.__radius * np.exp(1j * t)
             tmp = self.__multiSolve(vHat, phi)
 
-            A0 += tmp * np.power(phi, 0) * np.exp(1j * t)
-            A1 += tmp * np.power(phi, 1) * np.exp(1j * t)
+            A0_partial += tmp * np.power(phi, 0) * np.exp(1j * t)
+            A1_partial += tmp * np.power(phi, 1) * np.exp(1j * t)
+        MPI.COMM_WORLD.Barrier()
+
+        # Synchronize parital results
+        A0 = np.matlib.zeros(vHat.shape, dtype=complex)
+        A1 = np.matlib.zeros(vHat.shape, dtype=complex)
+        MPI.COMM_WORLD.Allreduce(A0_partial, A0, op=MPI.SUM)
+        MPI.COMM_WORLD.Allreduce(A1_partial, A1, op=MPI.SUM)
 
         # Final scale
         A0 *= self.__radius / self.__nNode
@@ -204,7 +231,6 @@ class Beyn:
         B -- the matrix of RHS
         w -- the complex frequency to use
         """
-
         # Solve
         self.__operator.solve(B, w)
 
@@ -218,7 +244,13 @@ class Beyn:
         # Done
         return X
 
-    @staticmethod
-    def __randomMatrix(n, m):
+    def __randomMatrix(self, n, m):
         """Returns a random complex matrix of the given size (n x m)"""
-        return np.matlib.rand(n, m) + np.matlib.rand(n, m) * 1j
+        if(self.__myProc == 0):
+            A = np.matlib.rand(n, m) + np.matlib.rand(n, m) * 1j
+        else:
+            A = None
+
+        A = MPI.COMM_WORLD.bcast(A, root=0)
+        MPI.COMM_WORLD.Barrier()
+        return A

bin/cim.py

@@ -11,6 +11,7 @@ See the LICENSE.txt and README.md for more license and copyright information.
 from solver       import GetDPWave
 from beyn         import Beyn
 from fixed_point  import FixedPoint
+from mpi4py       import MPI
 import args       as args
 import numpy      as np
 import time       as timer
@@ -21,6 +22,10 @@ tStart = timer.time()
 # Parse arguments
 arg = args.parse()
 
+# Who is master process and verbosity
+master  = (MPI.COMM_WORLD.Get_rank() == 0)
+verbose = (not arg.quiet) and (master)
+
 # Prepare variables of -setnumber for GetDP
 setnumber = []
 for i in arg.setnumber:
@@ -44,44 +49,46 @@ if(not arg.non_holomorphic):
 else:
     fixed   = FixedPoint(beyn, arg.nhMaxIt, arg.nhMinRratio, not arg.quiet)
                         # This variable ('fixed') is meant to be called
-    if(not arg.quiet):  # elsewhere (i.e. Python[] call in GetDP) thanks
+    if(verbose):        # elsewhere (i.e. Python[] call in GetDP) thanks
         fixed.display() # to the python interpreter.
         beyn.display()  # Please don't change its name!
 
     l, V, r = fixed.solve()
 
 # Display the computed eigenvalues
-if(not arg.quiet):
+if(verbose):
     print
     print "Eigenvalues:"
     print "----------- "
 
-for i in range(l.shape[0]):
-    if(not arg.quiet): print " # " + str(i) + ": ",
-    print "(" + format(np.real(l[i]).tolist(), '+.12e') + ")",
-    print "+",
-    print "(" + format(np.imag(l[i]).tolist(), '+.12e') + ")*j",
-    print "| " + format(r[i], 'e'),
-    print "| " + format(r[i]/np.abs(l[i]), 'e')
-
-if(not arg.quiet):
+#if(master):
+#    for i in range(l.shape[0]):
+#        if(verbose): print " # " + str(i) + ": ",
+#        print "(" + format(np.real(l[i]).tolist(), '+.12e') + ")",
+#        print "+",
+#        print "(" + format(np.imag(l[i]).tolist(), '+.12e') + ")*j",
+#        print "| " + format(r[i], 'e'),
+#        print "| " + format(r[i]/np.abs(l[i]), 'e')
+#
+if(verbose):
     print "----------- "
     print "Found: " + str(l.shape[0])
 
 # Generate GetDP views for eigenvectors
-if(not arg.quiet):
+if(verbose):
     print
     print "Generating eigenvector views..."
 
-for i in range(l.shape[0]):
-    operator.view(V[:, i], "cim" + str(i) + ".pos")
-    if(not arg.quiet): print " # View: " + str(i) + " done!"
-
+#if(master):
+#    for i in range(l.shape[0]):
+#        operator.view(V[:, i], "cim" + str(i) + ".pos")
+#        if(not arg.quiet): print " # View: " + str(i) + " done!"
+#
 # Stop timer
 tStop = timer.time()
 
 # Done
-if(not arg.quiet):
+if(verbose):
     print
     print "Elapsed time: " + str(round(tStop - tStart)) + "s"
     print "Bye!"

bin/fixed_point.py

@@ -11,6 +11,7 @@ See the LICENSE.txt and README.md for more license and copyright information.
 TODO: cluster identification / can a spline be holomorphic?
 """
 
+from   mpi4py import MPI
 import sys    as sys
 import numpy  as np
 
@@ -56,7 +57,8 @@ class FixedPoint:
         self.__solver      = solver
         self.__maxIt       = maxIt
         self.__minNHHRatio = minNHHRatio
-        self.__verbose     = verbose
+        self.__myProc      = MPI.COMM_WORLD.Get_rank()
+        self.__verbose     = verbose and (self.__myProc == 0)
 
         # State
         self.__cluster     = []

bin/solver.py

@@ -7,6 +7,7 @@ Technische Universitaet Darmstadt.
 See the LICENSE.txt and README.md for more license and copyright information.
 """
 
+from   mpi4py import MPI # Initializes correctly MPI if solver.py is standalone
 import getdp  as GetDP
 import numpy  as np
 

example/maxwell_linear/Makefile

@@ -5,7 +5,7 @@ test:
 	@echo "make ref: compute the first eigenvalue directly with GetDP"
 
 ref: init
-	getdp ref.pro -solve Eig -pos Eig -msh square.msh
+	getdp ref.pro -solve Ref -pos Ref -msh square.msh
 
 cim: init
 	python ${BIN}/cim.py maxwell.pro square.msh Maxwell 9e8 1e8 -lStart 1 -lStep 1

example/maxwell_linear/maxwell.pro

@@ -93,6 +93,7 @@ Resolution{
       { Name A; NameOfFormulation Maxwell; Type Complex; }
     }
     Operation{
+      MPI_SetCommSelf;
       Generate[A];
 
       If(doInit)

example/maxwell_linear/ref.pro

@@ -56,7 +56,7 @@ FunctionSpace{
 }
 
 Formulation{
-  { Name Eig; Type FemEquation;
+  { Name Ref; Type FemEquation;
     Quantity{
       { Name e; Type Local; NameOfSpace HCurl; }
     }
@@ -70,9 +70,9 @@ Formulation{
 }
 
 Resolution{
-  { Name Eig;
+  { Name Ref;
     System{
-      { Name A; NameOfFormulation Eig; Type Complex; }
+      { Name A; NameOfFormulation Ref; Type Complex; }
     }
     Operation{
       GenerateSeparate[A];
@@ -82,7 +82,7 @@ Resolution{
 }
 
 PostProcessing{
-  { Name Eig; NameOfFormulation Eig;
+  { Name Ref; NameOfFormulation Ref;
     Quantity{
       { Name e; Value{ Local{ [{e}]; In Omega; Jacobian JVol; } } }
     }
@@ -90,7 +90,7 @@ PostProcessing{
 }
 
 PostOperation{
-  { Name Eig; NameOfPostProcessing Eig;
+  { Name Ref; NameOfPostProcessing Ref;
     Operation{
       Print[e, OnElementsOf Omega, File "eig.pos", EigenvalueLegend];
     }

example/maxwell_sibc/cimResolution.pro

@@ -7,6 +7,7 @@ Resolution{
   { Name Solve;
     System{ { Name A; NameOfFormulation FormulationCIM; Type ComplexValue; } }
     Operation{
+      MPI_SetCommSelf;
       Generate[A];
 
       If(doInit)

example/maxwell_sibc_lagrange/cimResolution.pro

@@ -7,6 +7,7 @@ Resolution{
   { Name Solve;
     System{ { Name A; NameOfFormulation FormulationCIM; Type ComplexValue; } }
     Operation{
+      MPI_SetCommSelf;
       If(!doInit)
         Evaluate[Python[]{"real_corrector_init.py"}];
       EndIf