LoopControl: support vectorized loops in Fortran interface

3 files changed, 62 insertions, 40 deletions

diff --git a/Carpet/LoopControl/src/loopcontrol.F90 b/Carpet/LoopControl/src/loopcontrol.F90
index 86d3abd30..522e75857 100644
--- a/Carpet/LoopControl/src/loopcontrol.F90
+++ b/Carpet/LoopControl/src/loopcontrol.F90
@@ -4,15 +4,16 @@ module loopcontrol
   
   interface
      
-     subroutine lc_statmap_init (lc_lm, name)
+     subroutine lc_statmap_init (initialised, lc_lm, name)
        use loopcontrol_types
        implicit none
+       integer, intent(out) :: initialised
        type (lc_statmap_t) :: lc_lm
        character(*)        :: name
      end subroutine lc_statmap_init
      
      subroutine lc_control_init (lc_lc, lc_lm, &
-          imin,jmin,kmin, imax,jmax,kmax, ilsh,jlsh,klsh)
+          imin,jmin,kmin, imax,jmax,kmax, ilsh,jlsh,klsh,di)
        use loopcontrol_types
        implicit none
        type (lc_control_t) :: lc_lc
@@ -20,6 +21,7 @@ module loopcontrol
        integer, intent(in) :: imin, jmin, kmin
        integer, intent(in) :: imax, jmax, kmax
        integer, intent(in) :: ilsh, jlsh, klsh
+       integer, intent(in) :: di
      end subroutine lc_control_init
      
      subroutine lc_control_finish (lc_lc)
diff --git a/Carpet/LoopControl/src/loopcontrol_fortran.h b/Carpet/LoopControl/src/loopcontrol_fortran.h
index 94020ac28..dcde8e149 100644
--- a/Carpet/LoopControl/src/loopcontrol_fortran.h
+++ b/Carpet/LoopControl/src/loopcontrol_fortran.h
@@ -7,11 +7,12 @@
 
 #define LC_DECLARE3(name, i,j,k)                        &&\
 type (lc_statmap_t), save :: name/**/_lm                &&\
-logical, save :: name/**/_initialised = .false.         &&\
+integer, save :: name/**/_initialised = 0               &&\
 type (lc_control_t) :: name/**/_lc                      &&\
 integer :: name/**/_ii, name/**/_jj, name/**/_kk        &&\
 integer :: name/**/_imin, name/**/_jmin, name/**/_kmin  &&\
 integer :: name/**/_imax, name/**/_jmax, name/**/_kmax  &&\
+integer :: name/**/_ipos, name/**/_ioffset, name/**/_di &&\
 integer :: i, j, k
 
 
@@ -19,51 +20,69 @@ integer :: i, j, k
 #define LC_PRIVATE3(name) \
 name/**/_lc, \
 name/**/_imin, name/**/_jmin, name/**/_kmin, \
-name/**/_imax, name/**/_jmax, name/**/_kmax
-
-
-
-#define LC_LOOP3(name, i,j,k, imin,jmin,kmin, imax,jmax,kmax, ilsh,jlsh,klsh)                   &&\
-if (.not. name/**/_initialised) then                                                            &&\
-!$omp single                                                                                    &&\
-   call lc_statmap_init (name/**/_lm, "name")                                                   &&\
-!$omp end single                                                                                &&\
-!$omp single                                                                                    &&\
-   /* Set this flag only after initialising */                                                  &&\
-   name/**/_initialised = .true.                                                                &&\
-!$omp end single                                                                                &&\
-end if                                                                                          &&\
-call lc_control_init (name/**/_lc, name/**/_lm, imin,jmin,kmin, imax,jmax,kmax, ilsh,jlsh,klsh) &&\
-                                                                                                &&\
-/* Coarse loop */                                                                               &&\
-do name/**/_kk = name/**/_lc%kkmin + 1, name/**/_lc%kkmax, name/**/_lc%kkstep                   &&\
-   name/**/_kmin = name/**/_kk + name/**/_lc%kkkkmax                                            &&\
-   name/**/_kmax = min (name/**/_kmin - 1 + name/**/_lc%kkkkstep, name/**/_lc%kkmax)            &&\
-   do name/**/_jj = name/**/_lc%jjmin + 1, name/**/_lc%jjmax, name/**/_lc%jjstep                &&\
-      name/**/_jmin = name/**/_jj + name/**/_lc%jjjjmax                                         &&\
-      name/**/_jmax = min (name/**/_jj - 1 + name/**/_lc%jjjjmax, name/**/_lc%jjmax)            &&\
-      do name/**/_ii = name/**/_lc%iimin + 1, name/**/_lc%iimax, name/**/_lc%iistep             &&\
-         name/**/_imin = name/**/_ii + name/**/_lc%iiiimax                                      &&\
-         name/**/_imax = min (name/**/_ii - 1 + name/**/_lc%iiiimax, name/**/_lc%iimax)         &&\
-                                                                                                &&\
-         /* Fine loop */                                                                        &&\
-         do k = name/**/_kmin, name/**/_kmax                                                    &&\
-            do j = name/**/_jmin, name/**/_jmax                                                 &&\
-               do i = name/**/_imin, name/**/_imax
-                  
-
-
-#define LC_ENDLOOP3(name)                       &&\
+name/**/_imax, name/**/_jmax, name/**/_kmax, \
+name/**/_ipos, name/**/_ioffset, name/**/_di 
+
+
+
+#define LC_LOOP3(name, i,j,k, imin,jmin,kmin, imax,jmax,kmax, ilsh,jlsh,klsh) \
+  LC_LOOP3VEC(name, i,j,k, imin,jmin,kmin, imax,jmax,kmax, ilsh,jlsh,klsh, 1)
+#define LC_ENDLOOP3(name)                       \
+  LC_ENDLOOP3VEC(name)
+
+#define LC_LOOP3VEC(name, i,j,k, imin_,jmin_,kmin_, imax_,jmax_,kmax_, ilsh_,jlsh_,klsh_, di_) &&\
+    if (name/**/_initialised .eq. 0) then                                                      &&\
+      call lc_statmap_init (name/**/_initialised, name/**/_lm, "name")                         &&\
+    end if                                                                                     &&\
+    name/**/_di = (di_)                                                                        &&\
+    call lc_control_init (name/**/_lc, name/**/_lm,                                              \
+                     (imin_), (jmin_), (kmin_),                                                  \
+                     (imax_), (jmax_), (kmax_),                                                  \
+                     (ilsh_), (jlsh_), (klsh_),                                                  \
+                     name/**/_di)                                                              &&\
+                                                                                               &&\
+    /* Coarse loop */                                                                          &&\
+name/**/_lc_loop3vec:                                                                            \
+     do name/**/_kk = name/**/_lc%kkmin + 1, name/**/_lc%kkmax, name/**/_lc%kkstep             &&\
+      name/**/_kmin = name/**/_kk - 1 + name/**/_lc%kkkkmin                                    &&\
+      name/**/_kmax = min (name/**/_kk - 1 + name/**/_lc%kkkkmax, name/**/_lc%kkmax)           &&\
+                                                                                               &&\
+      do name/**/_jj = name/**/_lc%jjmin + 1, name/**/_lc%jjmax, name/**/_lc%jjstep            &&\
+        name/**/_jmin = name/**/_jj - 1 + name/**/_lc%jjjjmin                                  &&\
+        name/**/_jmax = min (name/**/_jj - 1 + name/**/_lc%jjjjmax, name/**/_lc%jjmax)         &&\
+                                                                                               &&\
+        do name/**/_ii = name/**/_lc%iimin + 1, name/**/_lc%iimax, name/**/_lc%iistep          &&\
+          name/**/_imin = name/**/_ii - 1 + name/**/_lc%iiiimin                                &&\
+          name/**/_imax = min (name/**/_ii - 1 + name/**/_lc%iiiimax, name/**/_lc%iimax)       &&\
+                                                                                               &&\
+          /* Fine loop */                                                                      &&\
+          do k = name/**/_kmin + 1, name/**/_kmax                                              &&\
+            do j = name/**/_jmin + 1, name/**/_jmax                                            &&\
+              LC_PRELOOP_STATEMENTS                                                            &&\
+                name/**/_ipos =                                                                  \
+                  name/**/_imin + name/**/_lc%ilsh * (j - 1 + name/**/_lc%jlsh * (k-1))        &&\
+                /* round down to the next multiple of lc_di, assuming that lc_di is a power */ &&\
+                /* of 2, and integers are encoded as two's-complement */                       &&\
+                name/**/_ioffset = iand(name/**/_ipos, - name/**/_di) - name/**/_ipos          &&\
+                do i = name/**/_imin + name/**/_ioffset + 1, name/**/_imax, name/**/_di
+
+
+
+#define LC_ENDLOOP3VEC(name)                    &&\
                end do                           &&\
             end do                              &&\
+            LC_POSTLOOP_STATEMENTS              &&\
          end do                                 &&\
                                                 &&\
       end do                                    &&\
    end do                                       &&\
-end do                                          &&\
+end do name/**/_lc_loop3vec                     &&\
                                                 &&\
 call lc_control_finish (name/**/_lc)
 
-
+/* Pre- and post loop statements are inserted around the innermost
+   loop, which is executed serially. By default these are empty. */
+#define LC_PRELOOP_STATEMENTS
+#define LC_POSTLOOP_STATEMENTS
 
 #endif  /* #ifndef LOOPCONTROL_FORTRAN_H */
diff --git a/Carpet/LoopControl/src/loopcontrol_types.F90 b/Carpet/LoopControl/src/loopcontrol_types.F90
index 843d6bdd0..76fcaf910 100644
--- a/Carpet/LoopControl/src/loopcontrol_types.F90
+++ b/Carpet/LoopControl/src/loopcontrol_types.F90
@@ -28,6 +28,7 @@ module loopcontrol_types
      integer imin, jmin, kmin
      integer imax, jmax, kmax
      integer ilsh, jlsh, klsh
+     integer di
      
      ! Control settings for thread parallelism (useful for debugging)
      integer iiimin, jjjmin, kkkmin