Merge branch 'master' into remove_readers

Makefile

@@ -1,5 +1,4 @@
 MPI?=1
-VALGRIND?=1
 DEBUG?=0
 UNDERSCORE?=1
 CC?=mpicc
@@ -15,7 +14,7 @@ TESTDIR =$(SRCROOT)/tests
 
 TARGET=parRSB
 TESTS=$(TESTDIR)/con/con-test
-LIB=src/lib$(TARGET).so
+LIB=src/lib$(TARGET).a
 
 INCFLAGS=-I$(SRCDIR) -I$(GSLIBDIR)/include
 

README.md

@@ -22,27 +22,22 @@ mpirun -np 4 ./gmsh-test twistedrod.msh
 
 We provide a simple C interface to use parRSB as a library.
 
-```sh
-int parRSB_partMesh(long long *egl, long long *vl, int *negl,
-                    long long *eglin, long long *vlin, int neglin,
-                    int nve, int *opt, MPI_Comm comm)
+```C
+int parRSB_partMesh(int *part, long long *vtx, int nel, int nve,
+                    int *options, MPI_Comm comm);
 ```
 
 See example here, see `tests/con/con-test.c`.
 
 ### Parameters
 
 ```sh
-egl     (out)    ... local list of global element IDs
-vl      (out)    ... local list of vertex IDs for all elements in egl
-negl    (in/out) ... length of egl on input
-                     local partition size on output
-eglin   (in)     ... local list of global element IDs
-vlin    (in)     ... local list of verticies making up each element in eglin (adjacency structure) 
-neglin  (in)     ... length of eglin
-opt     (in)     ... additional parameters (to use defaults set opt[0] = 0)
-nve     (in)     ... number of vertices of a single element (has to be the same for all)
-comm    (in)     ... MPI Communicator (size determines number of partitions)
+part    (out)   ... Destination MPI rank for each element.
+vtx     (in)    ... Vertices of all the elements (size = nel *nve)
+nel     (in)    ... Total number of local elements to MPI rank
+opt     (in)    ... Additional parameters (to use defaults set opt[0] = 0)
+nve     (in)    ... Number of vertices of a single element (has to be the same for all)
+comm    (in)    ... MPI Communicator (size determines number of partitions)
 ```
 
-Note, any initial distribution of mesh elements (eglin) is valid. 
+Note, any initial distribution of mesh elements is valid. 

src/genmap-comm.c

@@ -82,3 +82,22 @@ void GenmapSplitComm(GenmapHandle h, GenmapComm *c, int bin) {
   crystal_init(&(h->cr), &((*c)->gsComm));
 }
 
+int GenmapCrystalInit(GenmapHandle h, GenmapComm c) {
+  crystal_init(&(h->cr), &(c->gsComm));
+  return 0;
+}
+
+int GenmapCrystalTransfer(GenmapHandle h, int field) {
+  if(field == GENMAP_ORIGIN)
+    sarray_transfer(struct GenmapElement_private, &(h->elementArray), origin, 0,
+                    &(h->cr));
+  else if(field == GENMAP_PROC)
+    sarray_transfer(struct GenmapElement_private, &(h->elementArray), proc, 0,
+                    &(h->cr));
+  return 0;
+}
+
+int GenmapCrystalFinalize(GenmapHandle h) {
+  crystal_free(&(h->cr));
+  return 0;
+}

src/genmap-impl.h

@@ -12,6 +12,13 @@
 
 #include "genmap.h"
 //
+// Fiedler fields
+//
+#define GENMAP_FIEDLER 0
+#define GENMAP_GLOBALID 1
+#define GENMAP_PROC 2
+#define GENMAP_ORIGIN 3
+//
 // GenmapComm
 //
 struct GenmapComm_private {
@@ -28,6 +35,7 @@ struct GenmapElement_private {
   GenmapLong globalId;
   GenmapLong vertices[8];
   GenmapInt proc;
+  GenmapInt origin;
 };
 //
 // GenmapElements: Create, Destroy

src/genmap-rsb.c

@@ -220,18 +220,40 @@ int GenmapFiedler(GenmapHandle h, GenmapComm c, int maxIter, int global) {
   return iter;
 }
 
-void GenmapSplitByMedian(GenmapHandle h, int *bin) {
+void GenmapSplitByGlobalId(GenmapHandle h) {
   GenmapLong start = GenmapGetLocalStartIndex(h);
   GenmapLong nel = GenmapGetNGlobalElements(h);
-  GenmapLong id = GenmapCommRank(GenmapGetLocalComm(h));
-  GenmapLong np = GenmapCommSize(GenmapGetLocalComm(h));
+  GenmapInt id = GenmapCommRank(GenmapGetLocalComm(h));
+  GenmapInt np = GenmapCommSize(GenmapGetLocalComm(h));
   GenmapInt lelt = GenmapGetNLocalElements(h);
   GenmapElements elements = GenmapGetElements(h);
 
-  if(id < (np + 1) / 2)
-    *bin = 0;
-  else
-    *bin = 1;
+  GenmapInt pNel = (GenmapInt)(nel / np);
+  GenmapInt nrem = (GenmapInt)(nel - pNel * np);
+  GenmapInt idCount = 0;
+  while(idCount * pNel + ((idCount < nrem) ? idCount : nrem) < start)
+    idCount++;
+
+  GenmapLong upLimit = idCount * pNel + ((idCount < nrem) ? idCount : nrem);
+  GenmapLong downLimit = start;
+  do {
+    GenmapInt end = upLimit - start < lelt ? (GenmapInt)(upLimit - start) : lelt;
+    GenmapInt i;
+    for(i = (GenmapInt)(downLimit - start); i < end; i++)
+      elements[i].proc = idCount - 1;
+    downLimit = upLimit;
+    idCount++;
+    upLimit = idCount * pNel + ((idCount < nrem) ? idCount : nrem);
+  } while(downLimit - start < lelt);
+}
+
+void GenmapSplitByMedian(GenmapHandle h) {
+  GenmapLong start = GenmapGetLocalStartIndex(h);
+  GenmapLong nel = GenmapGetNGlobalElements(h);
+  GenmapInt id = GenmapCommRank(GenmapGetLocalComm(h));
+  GenmapInt np = GenmapCommSize(GenmapGetLocalComm(h));
+  GenmapInt lelt = GenmapGetNLocalElements(h);
+  GenmapElements elements = GenmapGetElements(h);
 
   GenmapInt pNel = (GenmapInt)(nel / np);
   GenmapInt nrem = (GenmapInt)(nel - pNel * np);
@@ -256,13 +278,13 @@ void GenmapAssignBins(GenmapHandle h, int field, buffer *buf0) {
   GenmapElements elements = GenmapGetElements(h);
   GenmapInt lelt = GenmapGetNLocalElements(h);
 
-  if(field == 0) { // Fiedler
+  if(field == GENMAP_FIEDLER) { // Fiedler
     // sort locally according to Fiedler vector
     sarray_sort_2(struct GenmapElement_private, elements, (GenmapUInt)lelt, fiedler,
                   TYPE_DOUBLE, globalId, TYPE_LONG, buf0);
     // set the bin based on Fiedler vector
     GenmapSetFiedlerBin(h);
-  } else {
+  } else if(GENMAP_GLOBALID) {
     // sort locally according to globalId
     sarray_sort_2(struct GenmapElement_private, elements, (GenmapUInt)lelt,
                   globalId, TYPE_LONG, globalId, TYPE_LONG, buf0);
@@ -275,14 +297,14 @@ void GenmapTransferToBins(GenmapHandle h, int field, buffer *buf0) {
   GenmapElements elements = GenmapGetElements(h);
   GenmapInt lelt = GenmapGetNLocalElements(h);
 
-  if(field == 0) { // Fiedler
+  if(field == GENMAP_FIEDLER) { // Fiedler
     sarray_transfer(struct GenmapElement_private, &(h->elementArray), proc, 0,
                     &(h->cr));
     elements = GenmapGetElements(h);
     lelt = GenmapGetNLocalElements(h);
     sarray_sort_2(struct GenmapElement_private, elements, (GenmapUInt)lelt, fiedler,
                   TYPE_DOUBLE, globalId, TYPE_LONG, buf0);
-  } else {
+  } else if(field == GENMAP_GLOBALID) {
     sarray_transfer(struct GenmapElement_private, &(h->elementArray), proc, 0,
                     &(h->cr));
     elements = GenmapGetElements(h);
@@ -299,6 +321,12 @@ void GenmapBinSort(GenmapHandle h, int field, buffer *buf0) {
   GenmapAssignBins(h, field, buf0);
   GenmapTransferToBins(h, field, buf0);
   GenmapScan(h, GenmapGetLocalComm(h));
+  if(field == GENMAP_FIEDLER) {
+    GenmapSplitByMedian(h);
+  } else if(field == GENMAP_GLOBALID) {
+    GenmapSplitByGlobalId(h);
+  }
+  GenmapTransferToBins(h, field, buf0);
 }
 
 void GenmapRSB(GenmapHandle h) {
@@ -331,10 +359,11 @@ void GenmapRSB(GenmapHandle h) {
     } while(ipass < npass && iter == maxIter);
 
     GenmapBinSort(h, 0, &buf0);
-
     int bin;
-    GenmapSplitByMedian(h, &bin);
-    GenmapTransferToBins(h, 0, &buf0);
+    GenmapInt np = GenmapCommSize(GenmapGetLocalComm(h));
+    GenmapInt id = GenmapCommRank(GenmapGetLocalComm(h));
+    if(id < (np + 1) / 2) bin = 0;
+    else bin = 1;
 
     GenmapComm c = GenmapGetLocalComm(h);
     GenmapSplitComm(h, &c, bin);

src/genmap.h

@@ -100,6 +100,9 @@ int GenmapGop(GenmapComm c, void *v, GenmapInt size, GenmapDataType type,
               GenmapInt op);
 int GenmapDestroyComm(GenmapComm c);
 void GenmapSplitComm(GenmapHandle h, GenmapComm *c, int bin);
+int GenmapCrystalInit(GenmapHandle h, GenmapComm c);
+int GenmapCrystalTransfer(GenmapHandle h, int field);
+int GenmapCrystalFinalize(GenmapHandle h);
 //
 // Function to read/write from/to FILE
 //

src/parRSB.c

@@ -8,9 +8,8 @@
 #include "genmap-io.h"
 #include "parRSB.h"
 
-void fparRSB_partMesh(long long *egl, long long *vl, int *negl,
-                      long long *eglcon, long long *vlcon, int *neglcon,
-                      int *nve, int *options, int *comm, int *err) {
+void fparRSB_partMesh(int *part, long long *vtx, int *nel, int *nve,
+                      int *options, int *comm, int *err) {
   *err = 1;
 
   GenmapCommExternal c;
@@ -20,14 +19,11 @@ void fparRSB_partMesh(long long *egl, long long *vl, int *negl,
   c = 0;
 #endif
 
-  *err = parRSB_partMesh(egl, vl, negl,
-                         eglcon, vlcon, *neglcon,
-                         *nve, options, c);
+  *err = parRSB_partMesh(part, vtx, *nel, *nve, options, c);
 }
 
-int parRSB_partMesh(long long *egl, long long *vl, int *negl,
-                    long long *eglcon, long long *vlcon, int neglcon,
-                    int nve, int *options, MPI_Comm comm) {
+int parRSB_partMesh(int *part, long long *vtx, int nel, int nve, int *options,
+                    MPI_Comm comm) {
   GenmapHandle h;
   GenmapInit(&h, comm);
 
@@ -36,43 +32,38 @@ int parRSB_partMesh(long long *egl, long long *vl, int *negl,
     h->printStat = options[2];
   }
 
-  // Check if negl is large enough
-  GenmapLong neglcon_ = (GenmapLong) neglcon;
-  GenmapGop(GenmapGetGlobalComm(h), &neglcon_, 1, GENMAP_LONG, GENMAP_SUM);
-  GenmapInt negl_max = (GenmapInt)(neglcon_ / GenmapCommSize(
-                                     GenmapGetGlobalComm(h))) + 1;
-  if(negl_max > *negl) {
-    printf("ERROR: negl to small to hold resulting partition!\n");
-    return 1;
-  }
+  // Assert that nel is greater then zero. Will be remove in future.
+  assert(nel > 0);
 
-  GenmapSetNLocalElements(h, (GenmapInt)neglcon);
+  GenmapSetNLocalElements(h, (GenmapInt)nel);
   GenmapScan(h, GenmapGetGlobalComm(h));
   GenmapSetNVertices(h, nve);
 
+  GenmapInt id = GenmapCommRank(GenmapGetGlobalComm(h));
   GenmapElements e = GenmapGetElements(h);
+  GenmapLong start = GenmapGetLocalStartIndex(h);
   GenmapInt i, j;
 
-  for(i = 0; i < neglcon; i++) {
-    e[i].globalId = eglcon[i];
+  for(i = 0; i < nel; i++) {
+    e[i].globalId = start + i;
+    e[i].origin = id;
     for(j = 0; j < nve; j++) {
-      e[i].vertices[j] = vlcon[i * nve + j];
+      e[i].vertices[j] = vtx[i * nve + j];
     }
   }
 
   GenmapRSB(h);
 
-  GenmapElements elements = GenmapGetElements(h);
-  *negl = GenmapGetNLocalElements(h);
+  GenmapCrystalInit(h, GenmapGetGlobalComm(h));
+  GenmapCrystalTransfer(h, GENMAP_ORIGIN);
+  GenmapCrystalFinalize(h);
 
-  for(i = 0; i < *negl; i++) {
-    egl[i] = elements[i].globalId;
-    for(j = 0; j < nve; j++) {
-      vl[nve * i + j] = elements[i].vertices[j];
-    }
-  }
+  // This should hold true
+  assert(GenmapGetNLocalElements(h) == nel);
 
-  if(h->printStat > 0) GenmapPartitionQuality(h);
+  for(i = 0; i < nel; i++) {
+    part[i] = e[i].proc;
+  }
 
   GenmapFinalize(h);
 

src/parRSB.h

@@ -1,12 +1,9 @@
 // API
-
 #include "genmap-gslib.h"
 
 #define fparRSB_partMesh FORTRAN_UNPREFIXED(fparrsb_partmesh,FPARRSB_PARTMESH)
-void fparRSB_partMesh(long long *egl, long long *vl, int *negl,
-                      long long *eglcon, long long *vlcon, int *neglcon,
+void fparRSB_partMesh(int *part, long long *vtx, int *nel,
                       int *nve, int *options, int *comm, int *err);
 
-int parRSB_partMesh(long long *egl, long long *vl, int *negl,
-                    long long *eglcon, long long *vlcon, int neglcon,
-                    int nve, int *options, MPI_Comm comm);
+int parRSB_partMesh(int *part, long long *vtx, int nel, int nve,
+                    int *options, MPI_Comm comm);

tests/con/con-test.c

@@ -22,22 +22,16 @@ int main(int argc, char *argv[]) {
   int ierr = conRead(argv[1], &c, comm);
   if(ierr) goto quit;
 
-  int nel_max = c.nelg / np + 1;
-  long long *el = (long long*) malloc(nel_max * sizeof(long long));
-  long long *vl = (long long*) malloc(nel_max * c.nv * sizeof(long long));
+  int *el = (int*) malloc(c.nel * sizeof(int));
 
-  int nelo = nel_max;
   options[0] = 1; // use custom options
   options[1] = 5; // debug level
   options[2] = 1; // print statistics
-  ierr = parRSB_partMesh(el, vl, &nelo, c.el, c.vl, c.nel, c.nv, options,
-                         comm);
+  ierr = parRSB_partMesh(el, c.vl, c.nel, c.nv, options, comm);
   if(ierr) goto quit;
 
-
   conFree(&c);
   free(el);
-  free(vl);
 
 quit:
   MPI_Finalize();