Update to gslib.h and allow to set a communicator on input (Nek5000#603)

core/comm_mpi.f

@@ -1,11 +1,13 @@
 c-----------------------------------------------------------------------
-      subroutine setupcomm()
+      subroutine setupcomm(comm,newcomm,newcommg,path_in,session_in)
       include 'mpif.h'
       include 'SIZE'
       include 'PARALLEL' 
       include 'TSTEP' 
       include 'INPUT'
 
+      integer comm, newcomm, newcommg
+      character session_in*(*), path_in*(*)
       logical flag
 
       common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
@@ -18,48 +20,58 @@ subroutine setupcomm()
 
       integer*8 ntags
 
-      ! Init MPI
       call mpi_initialized(mpi_is_initialized, ierr)
       if (.not.mpi_is_initialized) call mpi_init(ierr)
-      call mpi_comm_size(MPI_COMM_WORLD,np_global,ierr)
-      call mpi_comm_rank(MPI_COMM_WORLD,nid_global,ierr)
+
+      call mpi_comm_dup(comm,newcommg,ierr)
+      newcomm = newcommg
+      nekcomm = newcommg 
+
+      call mpi_comm_size(nekcomm,np_global,ierr)
+      call mpi_comm_rank(nekcomm,nid_global,ierr)
 
       ! check upper tag size limit
-      call mpi_comm_get_attr(MPI_COMM_WORLD,MPI_TAG_UB,ntags,flag,ierr)
+      call mpi_comm_get_attr(nekcomm,MPI_TAG_UB,ntags,flag,ierr)
       if (ntags .lt. np_global) then
          if(nid_global.eq.0) write(6,*) 'ABORT: MPI_TAG_UB too small!'
          call exitt
       endif
 
       ! set defaults
       nid         = nid_global
-      nekcomm     = MPI_COMM_WORLD
-      iglobalcomm = MPI_COMM_WORLD 
       ifneknek    = .false.
       ifneknekc   = .false. ! session are uncoupled
-      ifneknekm   = .false. ! not moving
       nsessions   = 1
 
       ierr = 0
       nlin = 0
       if (nid .eq. 0) then
-         write(6,*) 'Reading session file ...'
-         open (unit=8,file='SESSION.NAME',status='old',err=24)
- 21      read (8,*,END=22)
-         nlin = nlin + 1 
-         goto 21
- 22      rewind(8)
-         if (nlin.gt.2) read(8,*,err=24) nsessions
-         if (nsessions.gt.1) read(8,*,err=24) ifneknekc
-         do n=0,nsessions-1
-            call blank(session_mult(n),132)
-            call blank(path_mult(n)   ,132)
-            read(8,11,err=24) session_mult(n)
-            read(8,11,err=24) path_mult(n)
-            if (nsessions.gt.1) read(8,*,err=24)  npsess(n)
-         enddo
- 11      format(a132)
-         close(8)
+         l = ltrunc(session_in,len(session_in))
+         if (l .gt. 0) then
+            call blank(session_mult(0),132)
+            call chcopy(session_mult(0), session_in, l)
+            l = ltrunc(path_in,len(path_in))
+            call blank(path_mult(0)   ,132)
+            call chcopy(path_mult(0), path_in, l)
+         else
+           write(6,*) 'Reading session file ...'
+           open (unit=8,file='SESSION.NAME',status='old',err=24)
+ 21        read (8,*,END=22)
+           nlin = nlin + 1 
+           goto 21
+ 22        rewind(8)
+           if (nlin.gt.2) read(8,*,err=24) nsessions
+           if (nsessions.gt.1) read(8,*,err=24) ifneknekc
+           do n=0,nsessions-1
+              call blank(session_mult(n),132)
+              call blank(path_mult(n)   ,132)
+              read(8,11,err=24) session_mult(n)
+              read(8,11,err=24) path_mult(n)
+              if (nsessions.gt.1) read(8,*,err=24)  npsess(n)
+           enddo
+ 11        format(a132)
+           close(8)
+         endif
          write(6,*) 'Number of sessions:',nsessions
          goto 23
  24      ierr = 1
@@ -68,25 +80,22 @@ subroutine setupcomm()
       call err_chk(ierr,' Error while reading SESSION.NAME!$')
 
       call bcast(nsessions,ISIZE)
+      if (nsessions .gt. nsessmax) 
+     &   call exitti('nsessmax in SIZE too low!$',nsessmax)
+      if (nsessions .gt. 1) ifneknek = .true.
+
       call bcast(ifneknekc,LSIZE)
       do n = 0,nsessions-1
          call bcast(npsess(n),ISIZE)
          call bcast(session_mult(n),132*CSIZE)
          call bcast(path_mult(n),132*CSIZE)
       enddo
 
-      if (nsessions .gt. 1) ifneknek = .true.
-
-      if (nsessions .gt. nsessmax) 
-     &   call exitti('nsessmax in SIZE too low!$',nsessmax)
-
       ! single session run
       if (.not.ifneknek) then
          ifneknekc = .false.
          session   = session_mult(0)
          path      = path_mult(0)
-         call mpi_comm_dup(mpi_comm_world,iglobalcomm,ierr)
-         intracomm = iglobalcomm
          return
       endif
 
@@ -107,31 +116,15 @@ subroutine setupcomm()
          if (nid_global.ge.nid_global_root(n).and.
      &       nid_global.lt.nid_global_root_next) idsess = n
       enddo
-
-      call mpi_comm_split(mpi_comm_world,idsess,nid,intracomm,ierr)
+      call mpi_comm_split(comm,idsess,nid,newcomm,ierr)
 
       session = session_mult(idsess)
       path    = path_mult   (idsess)
 
-      ! setup intercommunication 
       if (ifneknekc) then
          if (nsessions.gt.2) call exitti(
      &     'More than 2 coupled sessions are currently not supported!$',
      $     nsessions)
-
-         if (idsess.eq.0) idsess_neighbor=1
-         if (idsess.eq.1) idsess_neighbor=0
-
-         call mpi_intercomm_create(intracomm,0,mpi_comm_world, 
-     &     nid_global_root(idsess_neighbor), 10,intercomm,ierr)
-
-         np_neighbor=npsess(idsess_neighbor)
-
-         ifhigh=.true.
-         call mpi_intercomm_merge(intercomm, ifhigh, iglobalcomm, ierr)
-
-         ngeom = 2  ! Initialize NEKNEK interface subiterations to 2.
-         ninter = 1 ! Initialize NEKNEK interface extrapolation order to 1.
       endif 
 
       return

core/crs_amg.c

@@ -5,18 +5,7 @@
 #include <math.h>
 #include <sys/stat.h>
 #include <limits.h>
-#include "c99.h"
-#include "name.h"
-#include "types.h"
-#include "fail.h"
-#include "mem.h"
-#include "sort.h"
-#include "sarray_sort.h"
-#include "gs_defs.h"
-#include "comm.h"
-#include "crystal.h"
-#include "sarray_transfer.h"
-#include "gs.h"
+#include "gslib.h"
 
 #define crs_setup PREFIXED_NAME(crs_amg_setup)
 #define crs_solve PREFIXED_NAME(crs_amg_solve)

core/crs_xxt.c

@@ -5,18 +5,7 @@
 #include <float.h>
 #include <string.h>
 #include <math.h>
-#include "c99.h"
-#include "name.h"
-#include "fail.h"
-#include "types.h"
-#include "tensor.h"
-#include "gs_defs.h"
-#include "comm.h"
-#include "mem.h"
-#include "sort.h"
-#include "sarray_sort.h"
-#include "sparse_cholesky.h"
-#include "gs.h"
+#include "gslib.h"
 
 #define crs_setup PREFIXED_NAME(crs_xxt_setup)
 #define crs_solve PREFIXED_NAME(crs_xxt_solve)
@@ -51,6 +40,166 @@ static unsigned lg(uint v)
 #undef CHECK
 }
 
+/* factors: L is in CSR format
+            D is a diagonal matrix stored as a vector
+   actual factorization is:
+
+                  -1      T
+     A   = (I-L) D   (I-L)
+
+      -1        -T        -1
+     A   = (I-L)   D (I-L)
+
+   (triangular factor is unit diagonal; the diagonal is not stored)
+*/
+struct sparse_cholesky {
+  uint n, *Lrp, *Lj;
+  double *L, *D;
+};
+
+/*
+  symbolic factorization: finds the sparsity structure of L
+
+  uses the concept of elimination tree:
+    the parent of node j is node i when L(i,j) is the first
+      non-zero in column j below the diagonal (i>j)
+    L's structure is discovered row-by-row; the first time
+      an entry in column j is set, it must be the parent
+
+  the nonzeros in L are the nonzeros in A + paths up the elimination tree
+
+  linear in the number of nonzeros of L
+*/
+static void factor_symbolic(uint n, const uint *Arp, const uint *Aj,
+                            struct sparse_cholesky *out, buffer *buf)
+{
+  uint *visit = tmalloc(uint,2*n), *parent = visit+n;
+  uint *Lrp, *Lj;
+  uint i,nz=0;
+
+  out->n=n;
+
+  for(i=0;i<n;++i) {
+    uint p=Arp[i], pe=Arp[i+1];
+    visit[i]=i, parent[i]=n;
+    for(;p!=pe;++p) {
+      uint j=Aj[p]; if(j>=i) break;
+      for(;visit[j]!=i;j=parent[j]) {
+        ++nz, visit[j]=i;
+        if(parent[j]==n) { parent[j]=i; break; }
+      }
+    }
+  }
+
+  Lrp=out->Lrp=tmalloc(uint,n+1+nz);
+  Lj =out->Lj =Lrp+n+1;
+
+  Lrp[0]=0;
+  for(i=0;i<n;++i) {
+    uint p=Arp[i], pe=Arp[i+1], count=0, *Ljr=&Lj[Lrp[i]];
+    visit[i]=i;
+    for(;p!=pe;++p) {
+      uint j=Aj[p]; if(j>=i) break;
+      for(;visit[j]!=i;j=parent[j]) Ljr[count++]=j, visit[j]=i;
+    }
+    sortv(Ljr, Ljr,count,sizeof(uint), buf);
+    Lrp[i+1]=Lrp[i]+count;
+  }
+  free(visit);
+}
+
+/*
+  numeric factorization:
+  
+  L is built row-by-row, using:    ( ' indicates transpose )
+
+  
+  [ A  r ]  = [ (I-L)   ] [ D^(-1)  ] [ (I-L)' -s ]
+  [ r' a ]    [  -s'  1 ] [     1/d ] [         1 ]
+            
+            = [ A   (I-L) D^(-1) (-s)  ]
+              [ r'  s' D^(-1) s + 1/d  ]
+              
+  so, if r' is the next row of A, up to but excluding the diagonal,
+  then the next row of L, s', obeys
+  
+     r = - (I-L) D^(-1) s
+ 
+  let y = (I-L)^(-1) (-r)
+  then s = D y, and d = 1/(s' y)
+  
+*/
+static void factor_numeric(uint n, const uint *Arp, const uint *Aj,
+                           const double *A,
+                           struct sparse_cholesky *out,
+                           uint *visit, double *y)
+{
+  const uint *Lrp=out->Lrp, *Lj=out->Lj;
+  double *D, *L;
+  uint i;
+
+  D=out->D=tmalloc(double,n+Lrp[n]);
+  L=out->L=D+n;
+
+  for(i=0;i<n;++i) {
+    uint p,pe; double a=0;
+    visit[i]=n;
+    for(p=Lrp[i],pe=Lrp[i+1];p!=pe;++p) {
+      uint j=Lj[p]; y[j]=0, visit[j]=i;
+    }
+    for(p=Arp[i],pe=Arp[i+1];p!=pe;++p) {
+      uint j=Aj[p];
+      if(j>=i) { if(j==i) a=A[p]; break; }
+      y[j]=-A[p];
+    }
+    for(p=Lrp[i],pe=Lrp[i+1];p!=pe;++p) {
+      uint q,qe,j=Lj[p]; double lij,yj=y[j];
+      for(q=Lrp[j],qe=Lrp[j+1];q!=qe;++q) {
+        uint k=Lj[q]; if(visit[k]==i) yj+=L[q]*y[k];
+      }
+      y[j]=yj;
+      L[p]=lij=D[j]*yj;
+      a-=yj*lij;
+    }
+    D[i]=1/a;
+  }
+}
+
+/* x = A^(-1) b;  works when x and b alias */
+void sparse_cholesky_solve(
+  double *x, const struct sparse_cholesky *fac, double *b)
+{
+  const uint n=fac->n, *Lrp=fac->Lrp, *Lj=fac->Lj;
+  const double *L=fac->L, *D=fac->D;
+  uint i, p,pe;
+  for(i=0;i<n;++i) {
+    double xi = b[i];
+    for(p=Lrp[i],pe=Lrp[i+1];p!=pe;++p) xi+=L[p]*x[Lj[p]];
+    x[i]=xi;
+  }
+  for(i=0;i<n;++i) x[i]*=D[i];
+  for(i=n;i;) {
+    double xi = x[--i];
+    for(p=Lrp[i],pe=Lrp[i+1];p!=pe;++p) x[Lj[p]]+=L[p]*xi;
+  }
+}
+
+void sparse_cholesky_factor(uint n, const uint *Arp, const uint *Aj,
+                            const double *A,
+                            struct sparse_cholesky *out, buffer *buf)
+{
+  const uint n_uints_as_dbls = (n*sizeof(uint)+sizeof(double)-1)/sizeof(double);
+  buffer_reserve(buf,(n_uints_as_dbls+n)*sizeof(double));
+  factor_symbolic(n,Arp,Aj,out,buf);
+  factor_numeric(n,Arp,Aj,A,out,buf->ptr,n_uints_as_dbls+(double*)buf->ptr);
+}
+
+void sparse_cholesky_free(struct sparse_cholesky *fac)
+{
+  free(fac->Lrp); fac->Lj=fac->Lrp=0;
+  free(fac->D);   fac->L =fac->D  =0;
+}
+
 struct csr_mat {
   uint n, *Arp, *Aj; double *A;
 };
@@ -848,4 +997,3 @@ void crs_free(struct xxt *data)
   free(data->vl);
   free(data);
 }
-

core/drive.f

@@ -1,6 +1,10 @@
       program NEKTON
-c
-      call nek_init(intracomm)
+
+      include 'mpif.h'
+      integer comm
+      comm = MPI_COMM_WORLD
+
+      call nek_init(comm)
       call nek_solve()
       call nek_end()