Lanczos with restart + make astyle

Makefile

@@ -100,7 +100,7 @@ astyle:
 	astyle --style=google --indent=spaces=2 --max-code-length=72 \
 	    --keep-one-line-statements --keep-one-line-blocks --lineend=linux \
             --suffix=none --preserve-date --formatted --pad-oper \
-	    --unpad-paren tests/*.[ch] src/*.[ch] inc/*.[ch] io/*/*.[ch]
+	    --unpad-paren tests/*/*.[ch] src/*.[ch] inc/*.[ch] io/*/*.[ch]
 
 print-%:
 	$(info [ variable name]: $*)

src/genmap-algo.c

@@ -173,6 +173,12 @@ int GenmapLanczos(GenmapHandle h, GenmapComm c, GenmapVector init,
   assert(alpha->size == beta->size + 1);
   assert(init->size == h->header->lelt);
 
+  if(h->header->nel < iter) {
+    iter = h->header->nel;
+    alpha->size = iter;
+    beta->size = iter - 1;
+  }
+
   GenmapVector q0, q1, u;
   GenmapScalar normq1 = 0., b = 0., sum;
 
@@ -238,8 +244,7 @@ int GenmapLanczos(GenmapHandle h, GenmapComm c, GenmapVector init,
       if(fabs(b) < normq1 * GENMAP_TOL) {
         beta->size = k;
         alpha->size = k + 1;
-        iter = k + 1;
-        break;
+        return k + 1;
       }
 
       GenmapCopyVector(q0, q1);
@@ -320,7 +325,7 @@ int GenmapFiedler(GenmapHandle h, GenmapComm c, int maxIter,
   GenmapElements elements = GenmapGetElements(h);
 
   GenmapInt i;
-  if(global) {
+  if(global > 0) {
     for(i = 0;  i < lelt; i++) {
       initVec->data[i] = (GenmapScalar) elements[i].globalId;
     }
@@ -428,7 +433,13 @@ void GenmapRSB(GenmapHandle h) {
   GenmapLong nel = h->header->nel;
   GenmapLong start = h->header->start;
   GenmapElements elements = GenmapGetElements(h);
-  int maxIter = 50;
+
+  int maxIter = 5;
+  maxIter = (nel / 5 > maxIter) ? nel / 5 : maxIter;
+  maxIter = (maxIter > 50) ? 50 : maxIter;
+
+  int iter = maxIter;
+  int npass = 100, ipass = 0;
 
   if(h->Id(h->global) == 0) printf("Running RSB ... ");
 #if defined(GENMAP_MPI)
@@ -448,13 +459,15 @@ void GenmapRSB(GenmapHandle h) {
   // must be initialized using the global communicator, we never
   // touch global communicator
   while(h->Np(h->local) > 1) {
-    if(h->Np(h->local) == h->Np(h->global))
-      maxIter = GenmapFiedler(h, h->local, maxIter, 1);
-    else
-      maxIter = GenmapFiedler(h, h->local, maxIter, 0);
-
+    int global = (h->Np(h->local) == h->Np(h->global));
+    do {
+      maxIter = GenmapFiedler(h, h->local, maxIter, global);
+      ipass++;
+      global = 0;
+    } while(ipass < npass && iter < maxIter);
     // sort locally according to Fiedler vector
-    sarray_sort_2(struct GenmapElement_private, elements, (GenmapUInt)lelt, fiedler,
+    sarray_sort_2(struct GenmapElement_private, elements, (GenmapUInt)lelt,
+                  fiedler,
                   TYPE_DOUBLE, globalId, TYPE_INT, &buf0);
 
 #if defined(GENMAP_DEBUG) && defined(GENMAP_MPI)
@@ -479,11 +492,13 @@ void GenmapRSB(GenmapHandle h) {
     elements = GenmapGetElements(h);
     lelt = h->header->lelt = (GenmapInt)h->elementArray.n;
     // sort locally again -- now we have everything sorted
-    sarray_sort_2(struct GenmapElement_private, elements, (GenmapUInt)lelt, fiedler,
+    sarray_sort_2(struct GenmapElement_private, elements, (GenmapUInt)lelt,
+                  fiedler,
                   TYPE_DOUBLE, globalId, TYPE_INT, &buf0);
 
     GenmapLong lelt_ = (GenmapLong)lelt;
-    comm_scan(out, &(h->local->gsComm), genmap_gs_long, gs_add, &lelt_, 1, buf);
+    comm_scan(out, &(h->local->gsComm), genmap_gs_long, gs_add, &lelt_, 1,
+              buf);
     start = h->header->start = out[0][0];
     nel = h->header->nel = out[1][0];
     id = h->Id(h->local);
@@ -496,19 +511,20 @@ void GenmapRSB(GenmapHandle h) {
     else
       bin = 1;
 
-    GenmapInt pNel = (GenmapInt) (nel / np);
-    GenmapInt nrem = (GenmapInt) (nel - pNel * np);
+    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);
+                                           nrem);
     GenmapLong downLimit = start;
     do {
-      GenmapInt end = upLimit - start < lelt ? (GenmapInt) (upLimit - start) : lelt;
+      GenmapInt end = upLimit - start < lelt ? (GenmapInt)(
+                        upLimit - start) : lelt;
       GenmapInt i;
-      for(i = (GenmapInt) (downLimit - start); i < end;
+      for(i = (GenmapInt)(downLimit - start); i < end;
           i++) elements[i].proc = idCount - 1;
       downLimit = upLimit;
       idCount++;
@@ -518,10 +534,11 @@ void GenmapRSB(GenmapHandle h) {
     sarray_transfer(struct GenmapElement_private, &(h->elementArray), proc,
                     0, &cr);
     elements = GenmapGetElements(h);
-    lelt = h->header->lelt = (GenmapInt) (h->elementArray.n);
+    lelt = h->header->lelt = (GenmapInt)(h->elementArray.n);
 
     // sort locally again -- now we have everything sorted
-    sarray_sort_2(struct GenmapElement_private, elements, (GenmapUInt)lelt, fiedler,
+    sarray_sort_2(struct GenmapElement_private, elements, (GenmapUInt)lelt,
+                  fiedler,
                   TYPE_DOUBLE, globalId, TYPE_INT, &buf0);
 
 
@@ -542,7 +559,8 @@ void GenmapRSB(GenmapHandle h) {
     crystal_init(&cr, &(h->local->gsComm));
 
     lelt_ = (GenmapLong)lelt;
-    comm_scan(out, &(h->local->gsComm), genmap_gs_long, gs_add, &lelt_, 1, buf);
+    comm_scan(out, &(h->local->gsComm), genmap_gs_long, gs_add, &lelt_, 1,
+              buf);
     start = h->header->start = out[0][0];
     nel = h->header->nel = out[1][0];
     id = h->Id(h->local);

src/genmap-comm.c

@@ -5,7 +5,7 @@ int GenmapCreateComm(GenmapComm *c, GenmapCommExternal ce) {
   comm_init(&(*c)->gsComm, ce);
   (*c)->verticesHandle = NULL;
   (*c)->laplacianWeights = NULL;
-  buffer_init(&(*c)->buf,1024);
+  buffer_init(&(*c)->buf, 1024);
   return 0;
 }
 

src/genmap-quality.c

@@ -23,7 +23,8 @@ GenmapInt GenmapPartitionQuality(GenmapHandle h) {
 
   GenmapComm c;
   GenmapCreateComm(&c, h->global->gsComm.c);
-  c->verticesHandle = gs_setup(data, numPoints, &c->gsComm, 0, gs_pairwise,
+  c->verticesHandle = gs_setup(data, numPoints, &c->gsComm, 0,
+                               gs_pairwise,
                                0);
 
   GenmapInt neighborsCount = 0;
@@ -57,13 +58,14 @@ GenmapInt GenmapPartitionQuality(GenmapHandle h) {
   GenmapGop(c, &ncSum, 1, GENMAP_INT, GENMAP_SUM);
 
   if(GenmapId(h->global) == 0) {
-    printf("Max neighbors: "GenmapIntFormat,ncMax);
-    printf(" | Min neighbors: "GenmapIntFormat,ncMin);
-    printf(" | Avg neighbors: "GenmapScalarFormat"\n",(1.0*ncSum)/GenmapNp(h->global));
+    printf("Max neighbors: "GenmapIntFormat, ncMax);
+    printf(" | Min neighbors: "GenmapIntFormat, ncMin);
+    printf(" | Avg neighbors: "GenmapScalarFormat"\n",
+           (1.0 * ncSum) / GenmapNp(h->global));
   }
 
   GenmapFree(data);
-  GenmapDestroyComm(c); 
+  GenmapDestroyComm(c);
 
   return neighborsCount;
 }

src/genmap-types.h

@@ -6,35 +6,35 @@
 // Genmap types
 //
 #if defined(GENMAP_LONG_LONG)
-  typedef long long GenmapLong;
-  typedef unsigned long long GenmapULong;
-  #define GenmapLongFormat "%lld"
-  #define GenmapULongFormat "%llu"
-  
-  #if defined(GENMAP_MPI)
-  #define GENMAP_LONG MPI_LONG_LONG
-  #define GENMAP_UNSIGNED_LONG MPI_UNSIGNED_LONG_LONG
-  #else
-  #define GENMAP_LONG 0
-  #define GENMAP_UNSIGNED_LONG 0
-  #endif
-  
-  #define genmap_gs_long gs_long_long
+typedef long long GenmapLong;
+typedef unsigned long long GenmapULong;
+#define GenmapLongFormat "%lld"
+#define GenmapULongFormat "%llu"
+
+#if defined(GENMAP_MPI)
+#define GENMAP_LONG MPI_LONG_LONG
+#define GENMAP_UNSIGNED_LONG MPI_UNSIGNED_LONG_LONG
+#else
+#define GENMAP_LONG 0
+#define GENMAP_UNSIGNED_LONG 0
+#endif
+
+#define genmap_gs_long gs_long_long
 #else
-  typedef long GenmapLong;
-  typedef unsigned long GenmapULong;
-  #define GenmapLongFormat "%ld"
-  #define GenmapULongFormat "%lu"
-  
-  #if defined(GENMAP_MPI)
-  #define GENMAP_LONG MPI_LONG
-  #define GENMAP_UNSIGNED_LONG MPI_UNSIGNED_LONG
-  #else
-  #define GENMAP_LONG 0
-  #define GENMAP_UNSIGNED_LONG 0
-  #endif
+typedef long GenmapLong;
+typedef unsigned long GenmapULong;
+#define GenmapLongFormat "%ld"
+#define GenmapULongFormat "%lu"
+
+#if defined(GENMAP_MPI)
+#define GENMAP_LONG MPI_LONG
+#define GENMAP_UNSIGNED_LONG MPI_UNSIGNED_LONG
+#else
+#define GENMAP_LONG 0
+#define GENMAP_UNSIGNED_LONG 0
+#endif
 
-  #define genmap_gs_long gs_long
+#define genmap_gs_long gs_long
 #endif
 
 typedef int GenmapInt;

src/genmap-vector.c

@@ -127,7 +127,7 @@ int GenmapScaleVector(GenmapVector y, GenmapVector x,
 
 int GenmapCreateOnesVector(GenmapVector *x, GenmapInt size) {
   GenmapCreateVector(x, size);
-  
+
   GenmapInt i;
   for(i = 0; i < size; i++) {
     (*x)->data[i] = 1.;

src/parRSB.c

@@ -8,10 +8,9 @@
 #include "genmap-io.h"
 #include "parRSB.h"
 
-void fparRSB_partMesh(long long *egl, long long *vl, int *negl, 
+void fparRSB_partMesh(long long *egl, long long *vl, int *negl,
                       long long *eglcon, long long *vlcon, int *neglcon,
-                      int *nve, int *comm, int *err) 
-{
+                      int *nve, int *comm, int *err) {
   *err = 1;
   setbuf(stdout, NULL);
 
@@ -29,8 +28,7 @@ void fparRSB_partMesh(long long *egl, long long *vl, int *negl,
 
 int parRSB_partMesh(long long *egl, long long *vl, int *negl,
                     long long *eglcon, long long *vlcon, int neglcon,
-                    int nve, MPI_Comm comm) 
-{
+                    int nve, MPI_Comm comm) {
   GenmapHandle h;
   GenmapInit(&h, comm, "interface");
 
@@ -59,12 +57,12 @@ int parRSB_partMesh(long long *egl, long long *vl, int *negl,
   h->elementArray.n = neglcon;
 
   GenmapElements e = GenmapGetElements(h);
-  GenmapInt i,j;
+  GenmapInt i, j;
 
   for(i = 0; i < neglcon; i++) {
     e[i].globalId = eglcon[i];
-    for(j = 0; j<nve; j++) {
-      e[i].vertices[j] = vlcon[i*nve+j];
+    for(j = 0; j < nve; j++) {
+      e[i].vertices[j] = vlcon[i * nve + j];
     }
   }
 

tests/con/con-test.c

@@ -10,8 +10,7 @@ Parition mesh using Nek5000's vertex connectivity (con) file.
 
 MPI_Comm comm = MPI_COMM_WORLD;
 
-int main(int argc, char *argv[])
-{
+int main(int argc, char *argv[]) {
   MPI_Init(&argc, &argv);
   int myid, np;
   MPI_Comm_rank(comm, &myid);
@@ -21,9 +20,9 @@ int main(int argc, char *argv[])
   int ierr = conRead(argv[1], &c, comm);
   if(ierr) goto quit;
 
-  int nel_max = c.nelg/np + 1;
+  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));
+  long long *vl = (long long*) malloc(nel_max * c.nv * sizeof(long long));
 
   int nelo = nel_max;
   ierr = parRSB_partMesh(el, vl, &nelo, c.el, c.vl, c.nel, c.nv, comm);

tests/con/conReader.h

@@ -1,14 +1,14 @@
-struct con{
+struct con {
   int nv, nelg, nel;
   long long *el, *vl;
 };
 
-void conFree(struct con *c){
+void conFree(struct con *c) {
   free(c->el);
   free(c->vl);
 }
 
-int conRead(char *fname, struct con *c, MPI_Comm comm){
+int conRead(char *fname, struct con *c, MPI_Comm comm) {
   int i, j, ierr;
   int myid, np;
   MPI_Comm_rank(comm, &myid);
@@ -30,21 +30,22 @@ int conRead(char *fname, struct con *c, MPI_Comm comm){
   float byte_test;
   MPI_File_read_all(fh, &byte_test, 1, MPI_FLOAT, MPI_STATUS_IGNORE);
   if(abs(byte_test - 6.543210) > 1e-7) {
-    if(myid == 0) printf("ERROR byte_test failed! %f\n",byte_test);
+    if(myid == 0) printf("ERROR byte_test failed! %f\n", byte_test);
     return 1;
   }
 
-  int nelr = nelgt/np;
-  for(i = 0; i < nelgt%np; ++i) if(np-i == myid) nelr++;
+  int nelr = nelgt / np;
+  for(i = 0; i < nelgt % np; ++i) if(np - i == myid) nelr++;
 
   int nelr_;
   MPI_Scan(&nelr, &nelr_, 1, MPI_INT, MPI_SUM, comm);
   long long off = sizeof(hdr) + sizeof(int);
-  off += (long long)(nelr_ - nelr)*(nv+1)*sizeof(int);
-  MPI_File_set_view(fh, off, MPI_BYTE, MPI_BYTE, "native", MPI_INFO_NULL); 
+  off += (long long)(nelr_ - nelr) * (nv + 1) * sizeof(int);
+  MPI_File_set_view(fh, off, MPI_BYTE, MPI_BYTE, "native", MPI_INFO_NULL);
 
-  int *buf = (int*) malloc(nelr*(nv+1) * sizeof(int)); 
-  ierr = MPI_File_read_all(fh, buf, (nv+1)*nelr, MPI_INT, MPI_STATUS_IGNORE);
+  int *buf = (int*) malloc(nelr * (nv + 1) * sizeof(int));
+  ierr = MPI_File_read_all(fh, buf, (nv + 1) * nelr, MPI_INT,
+                           MPI_STATUS_IGNORE);
   if(ierr != 0) {
     if(myid == 0) printf("ERROR: failure while reading!\n");
     return 1;
@@ -53,10 +54,10 @@ int conRead(char *fname, struct con *c, MPI_Comm comm){
   MPI_File_close(&fh);
 
   c->el = (long long*) malloc(nelr * sizeof(long long));
-  c->vl = (long long*) malloc(nv*nelr * sizeof(long long));
+  c->vl = (long long*) malloc(nv * nelr * sizeof(long long));
   for(i = 0; i < nelr; ++i) {
-    c->el[i] = buf[i*(nv+1)];
-    for(j = 0; j < nv; ++j) c->vl[i*nv + j] = buf[i*(nv+1)+1+j]; 
+    c->el[i] = buf[i * (nv + 1)];
+    for(j = 0; j < nv; ++j) c->vl[i * nv + j] = buf[i * (nv + 1) + 1 + j];
   }
 
   c->nv = nv;

tests/gmsh/gmsh-test.c

@@ -6,8 +6,7 @@
 #include <genmap-impl.h>
 
 
-int main(int argc, char **argv)
-{
+int main(int argc, char **argv) {
 #if defined(GENMAP_MPI)
   MPI_Init(&argc, &argv);
 #else