Optimize i7trivox.F: replace loops with intrinsics and fix boundary detection bug

engine/source/interfaces/intsort/i7trivox.F

@@ -214,10 +214,10 @@ SUBROUTINE I7TRIVOX(
       IF(ITASK==0.AND.TOTAL_NB_NRTM>0)THEN
         !==========LAGRANGIAN and STAGGERED ALE==========!
         TEST_V0 = FIRST_TEST(NIN)
+        IIX(1:NSN) = 0
+        IIY(1:NSN) = 0
+        IIZ(1:NSN) = 0
         DO I=1,NSN
-          IIX(I)=0
-          IIY(I)=0
-          IIZ(I)=0
           IF(STFN(I) == ZERO)CYCLE
           J=NSV(I)
           !C Optimization // search for nodes within xmin xmax of
@@ -228,17 +228,14 @@ SUBROUTINE I7TRIVOX(
           IF(X(2,J) > YMAX)  CYCLE
           IF(X(3,J) < ZMIN)  CYCLE
           IF(X(3,J) > ZMAX)  CYCLE
-          IIX(I)=INT(NBX*(X(1,J)-XMINB)/(XMAXB-XMINB))
-          IIY(I)=INT(NBY*(X(2,J)-YMINB)/(YMAXB-YMINB))
-          IIZ(I)=INT(NBZ*(X(3,J)-ZMINB)/(ZMAXB-ZMINB))
-          IIX(I)=MAX(1,2+MIN(NBX,IIX(I)))
-          IIY(I)=MAX(1,2+MIN(NBY,IIY(I)))
-          IIZ(I)=MAX(1,2+MIN(NBZ,IIZ(I)))
+          IIX(I)=MAX(1,MIN(NBX+2,2+INT(NBX*(X(1,J)-XMINB)/(XMAXB-XMINB))))
+          IIY(I)=MAX(1,MIN(NBY+2,2+INT(NBY*(X(2,J)-YMINB)/(YMAXB-YMINB))))
+          IIZ(I)=MAX(1,MIN(NBZ+2,2+INT(NBZ*(X(3,J)-ZMINB)/(ZMAXB-ZMINB))))
           FIRST = VOXEL(IIX(I),IIY(I),IIZ(I))
           IF(TEST_V0) THEN
 ! count the number of SECONDARY nodes outiside the reduced box
-            IF(IIX(I) ==  1    .OR.  IIY(I) ==  1    .OR. IIZ(I) == 1 .AND.
-     .         IIX(I) ==  NBX+2 .OR. IIY(I) == NBY+2 .OR. IIZ(I) == NBZ+2) THEN
+            IF(IIX(I)==1 .OR. IIY(I)==1 .OR. IIZ(I)==1 .OR.
+     .         IIX(I)==NBX+2 .OR. IIY(I)==NBY+2 .OR. IIZ(I)==NBZ+2) THEN
               CPT_VOX0 = CPT_VOX0 +1
             ENDIF
           ENDIF
@@ -277,18 +274,15 @@ SUBROUTINE I7TRIVOX(
 
           REMOTE_S_NODE = REMOTE_S_NODE + 1
           LIST_REMOTE_S_NODE( REMOTE_S_NODE ) = J
-          IIX(NSN+J)=INT(NBX*(XREM(1,J)-XMINB)/(XMAXB-XMINB))
-          IIY(NSN+J)=INT(NBY*(XREM(2,J)-YMINB)/(YMAXB-YMINB))
-          IIZ(NSN+J)=INT(NBZ*(XREM(3,J)-ZMINB)/(ZMAXB-ZMINB))
-          IIX(NSN+J)=MAX(1,2+MIN(NBX,IIX(NSN+J)))
-          IIY(NSN+J)=MAX(1,2+MIN(NBY,IIY(NSN+J)))
-          IIZ(NSN+J)=MAX(1,2+MIN(NBZ,IIZ(NSN+J)))
+          IIX(NSN+J)=MAX(1,MIN(NBX+2,2+INT(NBX*(XREM(1,J)-XMINB)/(XMAXB-XMINB))))
+          IIY(NSN+J)=MAX(1,MIN(NBY+2,2+INT(NBY*(XREM(2,J)-YMINB)/(YMAXB-YMINB))))
+          IIZ(NSN+J)=MAX(1,MIN(NBZ+2,2+INT(NBZ*(XREM(3,J)-ZMINB)/(ZMAXB-ZMINB))))
 
           FIRST = VOXEL(IIX(NSN+J),IIY(NSN+J),IIZ(NSN+J))
           IF(TEST_V0) THEN
 ! count the number of SECONDARY nodes outiside the reduced box
-            IF(IIX(J+NSN) ==  1     .OR. IIY(J+NSN) ==  1     .OR. IIZ(J+NSN) == 1 .AND.
-     .         IIX(J+NSN) ==  NBX+2 .OR. IIY(J+NSN) ==  NBY+2 .OR. IIZ(J+NSN) == NBZ+2) THEN
+            IF(IIX(J+NSN)==1 .OR. IIY(J+NSN)==1 .OR. IIZ(J+NSN)==1 .OR.
+     .         IIX(J+NSN)==NBX+2 .OR. IIY(J+NSN)==NBY+2 .OR. IIZ(J+NSN)==NBZ+2) THEN
               CPT_VOX0 = CPT_VOX0 +1
             ENDIF
           ENDIF
@@ -359,24 +353,24 @@ SUBROUTINE I7TRIVOX(
         M3 = IRECT(3,NE)
         M4 = IRECT(4,NE)
 
+        ! Extract coordinates and compute bounds using intrinsic functions
         XX1=X(1,M1)
         XX2=X(1,M2)
         XX3=X(1,M3)
         XX4=X(1,M4)
-        XMAXE=MAX(XX1,XX2,XX3,XX4)
-        XMINE=MIN(XX1,XX2,XX3,XX4)
-
         YY1=X(2,M1)
         YY2=X(2,M2)
         YY3=X(2,M3)
         YY4=X(2,M4)
-        YMAXE=MAX(YY1,YY2,YY3,YY4)
-        YMINE=MIN(YY1,YY2,YY3,YY4)
-
         ZZ1=X(3,M1)
         ZZ2=X(3,M2)
         ZZ3=X(3,M3)
         ZZ4=X(3,M4)
+
+        XMAXE=MAX(XX1,XX2,XX3,XX4)
+        XMINE=MIN(XX1,XX2,XX3,XX4)
+        YMAXE=MAX(YY1,YY2,YY3,YY4)
+        YMINE=MIN(YY1,YY2,YY3,YY4)
         ZMAXE=MAX(ZZ1,ZZ2,ZZ3,ZZ4)
         ZMINE=MIN(ZZ1,ZZ2,ZZ3,ZZ4)
 
@@ -388,37 +382,29 @@ SUBROUTINE I7TRIVOX(
 
         !index of voxels occupied by the facet
         IF(NBX>1) THEN
-          IX1=INT(NBX*(XMINE-AAA-XMINB)/(XMAXB-XMINB))
-          IX2=INT(NBX*(XMAXE+AAA-XMINB)/(XMAXB-XMINB))
+          IX1=MAX(1,MIN(NBX+2,2+INT(NBX*(XMINE-AAA-XMINB)/(XMAXB-XMINB))))
+          IX2=MAX(1,MIN(NBX+2,2+INT(NBX*(XMAXE+AAA-XMINB)/(XMAXB-XMINB))))
         ELSE
-          IX1=-2
-          IX2=1
+          IX1=1
+          IX2=NBX+2
         ENDIF
 
         IF(NBY>1) THEN
-          IY1=INT(NBY*(YMINE-AAA-YMINB)/(YMAXB-YMINB))
-          IY2=INT(NBY*(YMAXE+AAA-YMINB)/(YMAXB-YMINB))
+          IY1=MAX(1,MIN(NBY+2,2+INT(NBY*(YMINE-AAA-YMINB)/(YMAXB-YMINB))))
+          IY2=MAX(1,MIN(NBY+2,2+INT(NBY*(YMAXE+AAA-YMINB)/(YMAXB-YMINB))))
         ELSE
-          IY1=-2
-          IY2=1
+          IY1=1
+          IY2=NBY+2
         ENDIF
 
         IF(NBZ>1) THEN
-          IZ1=INT(NBZ*(ZMINE-AAA-ZMINB)/(ZMAXB-ZMINB))
-          IZ2=INT(NBZ*(ZMAXE+AAA-ZMINB)/(ZMAXB-ZMINB))
+          IZ1=MAX(1,MIN(NBZ+2,2+INT(NBZ*(ZMINE-AAA-ZMINB)/(ZMAXB-ZMINB))))
+          IZ2=MAX(1,MIN(NBZ+2,2+INT(NBZ*(ZMAXE+AAA-ZMINB)/(ZMAXB-ZMINB))))
         ELSE
-          IZ1=-2
-          IZ2=1
+          IZ1=1
+          IZ2=NBZ+2
         ENDIF
 
-        IX1=MAX(1,2+MIN(NBX,IX1))
-        IY1=MAX(1,2+MIN(NBY,IY1))
-        IZ1=MAX(1,2+MIN(NBZ,IZ1))
-
-        IX2=MAX(1,2+MIN(NBX,IX2))
-        IY2=MAX(1,2+MIN(NBY,IY2))
-        IZ2=MAX(1,2+MIN(NBZ,IZ2))
-
         !nbpelem = 0
         !nnpelem = 0
         !nnr0pelem = 0

engine/source/mpi/generic/spmd_alltoall.F90

@@ -96,9 +96,9 @@ subroutine spmd_alltoall_ints(sendbuf, sendcount, recvbuf, recvcount, comm)
           tag = 0
           call spmd_in(tag)
           if (present(comm)) then
-            call MPI_Alltoall(sendbuf, sendcount, MPI_INT, recvbuf, recvcount, MPI_INT, comm, ierr)
+            call MPI_Alltoall(sendbuf, sendcount, MPI_INTEGER, recvbuf, recvcount, MPI_INTEGER, comm, ierr)
           else
-            call MPI_Alltoall(sendbuf, sendcount, MPI_INT, recvbuf, recvcount, MPI_INT, SPMD_COMM_WORLD, ierr)
+            call MPI_Alltoall(sendbuf, sendcount, MPI_INTEGER, recvbuf, recvcount, MPI_INTEGER, SPMD_COMM_WORLD, ierr)
           end if
           call spmd_out(tag,ierr)
 #else
@@ -195,10 +195,10 @@ subroutine spmd_alltoall_int(sendbuf, sendcount, recvbuf, recvcount, comm)
           tag = 0
           call spmd_in(tag)
           if (present(comm)) then
-            call MPI_Alltoall(sendbuf, sendcount, MPI_INT, recvbuf, recvcount, MPI_INT, comm, ierr)
+            call MPI_Alltoall(sendbuf, sendcount, MPI_INTEGER, recvbuf, recvcount, MPI_INTEGER, comm, ierr)
           else
-            call MPI_Alltoall(sendbuf, sendcount, MPI_INT, recvbuf, recvcount,&
-              MPI_INT, SPMD_COMM_WORLD, ierr)
+            call MPI_Alltoall(sendbuf, sendcount, MPI_INTEGER, recvbuf, recvcount,&
+              MPI_INTEGER, SPMD_COMM_WORLD, ierr)
           end if
           call spmd_out(tag,ierr)
 #endif

engine/source/mpi/spmd_allgatherv.F90

@@ -98,9 +98,9 @@ subroutine spmd_allgatherv_ints(sendbuf, sendcount, recvbuf, recvcounts, displs,
           tag = 0
           call spmd_in(tag)
           if (present(comm)) then
-            call MPI_Allgatherv(sendbuf, sendcount, MPI_INT, recvbuf, recvcounts, displs, MPI_INT, comm, ierr)
+            call MPI_Allgatherv(sendbuf, sendcount, MPI_INTEGER, recvbuf, recvcounts, displs, MPI_INTEGER, comm, ierr)
           else
-            call MPI_Allgatherv(sendbuf, sendcount, MPI_INT, recvbuf, recvcounts, displs, MPI_INT, SPMD_COMM_WORLD, ierr)
+            call MPI_Allgatherv(sendbuf, sendcount, MPI_INTEGER, recvbuf, recvcounts, displs, MPI_INTEGER, SPMD_COMM_WORLD, ierr)
           end if
           call spmd_out(tag,ierr)
 #else
@@ -176,6 +176,8 @@ subroutine spmd_allgatherv_double(sendbuf, sendcount, recvbuf, recvcounts, displ
               recvcounts, displs, MPI_DOUBLE_PRECISION, SPMD_COMM_WORLD, ierr)
           end if
           call spmd_out(tag,ierr)
+#else
+          recvbuf(1) = sendbuf
 #endif
         end subroutine spmd_allgatherv_double
 !||====================================================================
@@ -202,12 +204,14 @@ subroutine spmd_allgatherv_int(sendbuf, sendcount, recvbuf, recvcounts, displs,
           tag = 0
           call spmd_in(tag)
           if (present(comm)) then
-            call MPI_Allgatherv(sendbuf, sendcount, MPI_INT, recvbuf, recvcounts, displs, MPI_INT, comm, ierr)
+            call MPI_Allgatherv(sendbuf, sendcount, MPI_INTEGER, recvbuf, recvcounts, displs, MPI_INTEGER, comm, ierr)
           else
-            call MPI_Allgatherv(sendbuf, sendcount, MPI_INT, recvbuf, recvcounts, &
-              displs, MPI_INT, SPMD_COMM_WORLD, ierr)
+            call MPI_Allgatherv(sendbuf, sendcount, MPI_INTEGER, recvbuf, recvcounts, &
+              displs, MPI_INTEGER, SPMD_COMM_WORLD, ierr)
           end if
           call spmd_out(tag,ierr)
+#else
+          recvbuf(1) = sendbuf
 #endif
         end subroutine spmd_allgatherv_int
 !||====================================================================
@@ -239,6 +243,8 @@ subroutine spmd_allgatherv_real(sendbuf, sendcount, recvbuf, recvcounts, displs,
             call MPI_Allgatherv(sendbuf, sendcount, MPI_REAL, recvbuf, recvcounts, displs, MPI_REAL, SPMD_COMM_WORLD, ierr)
           end if
           call spmd_out(tag,ierr)
+#else
+          recvbuf(1) = sendbuf
 #endif
         end subroutine spmd_allgatherv_real
 

engine/source/mpi/spmd_mod.F90

@@ -503,7 +503,7 @@ subroutine spmd_reduce_ints(sendbuf, recvbuf, buf_count, operation, root, comm)
             used_comm = SPMD_COMM_WORLD
           end if
 
-          call MPI_Reduce(sendbuf, recvbuf, buf_count, MPI_REAL, mpi_op, root, used_comm, ierr)
+          call MPI_Reduce(sendbuf, recvbuf, buf_count, MPI_INTEGER, mpi_op, root, used_comm, ierr)
           call spmd_out(TAG_REDUCE,ierr)
 #else
           recvbuf(1:buf_count) = sendbuf(1:buf_count)
@@ -539,7 +539,7 @@ subroutine spmd_reduce_doubles(sendbuf, recvbuf, buf_count, operation, root, com
             used_comm = SPMD_COMM_WORLD
           end if
 
-          call MPI_Reduce(sendbuf, recvbuf, buf_count, MPI_REAL, mpi_op, root, used_comm, ierr)
+          call MPI_Reduce(sendbuf, recvbuf, buf_count, MPI_DOUBLE_PRECISION, mpi_op, root, used_comm, ierr)
           call spmd_out(TAG_REDUCE,ierr)
 #else
           recvbuf(1:buf_count) = sendbuf(1:buf_count)
@@ -718,7 +718,7 @@ subroutine spmd_reduce_int(sendbuf, recvbuf, buf_count, operation, root, comm)
             used_comm = SPMD_COMM_WORLD
           end if
 
-          call MPI_Reduce(sendbuf, recvbuf, buf_count, MPI_REAL, mpi_op, root, used_comm, ierr)
+          call MPI_Reduce(sendbuf, recvbuf, buf_count, MPI_INTEGER, mpi_op, root, used_comm, ierr)
           call spmd_out(TAG_REDUCE,ierr)
 #else
           recvbuf = sendbuf  ! In case MPI is not defined, just copy the value
@@ -754,7 +754,7 @@ subroutine spmd_reduce_double(sendbuf, recvbuf, buf_count, operation, root, comm
             used_comm = SPMD_COMM_WORLD
           end if
 
-          call MPI_Reduce(sendbuf, recvbuf, buf_count, MPI_REAL, mpi_op, root, used_comm, ierr)
+          call MPI_Reduce(sendbuf, recvbuf, buf_count, MPI_DOUBLE_PRECISION, mpi_op, root, used_comm, ierr)
           call spmd_out(TAG_REDUCE,ierr)
 #else
           recvbuf = sendbuf  ! In case MPI is not defined, just copy the value
@@ -862,7 +862,7 @@ subroutine spmd_allreduce_real(sendbuf, recvbuf, buf_count, operation, comm)
             used_comm = SPMD_COMM_WORLD
           end if
 
-          call MPI_Allreduce(sendbuf, recvbuf, buf_count, MPI_DOUBLE_PRECISION, mpi_op, used_comm, ierr)
+          call MPI_Allreduce(sendbuf, recvbuf, buf_count, MPI_REAL, mpi_op, used_comm, ierr)
           call spmd_out(TAG_ALLREDUCE,ierr)
 #else
           recvbuf = sendbuf  ! In case MPI is not defined, just copy the value

engine/source/mpi/spmd_pack.F90

@@ -119,8 +119,8 @@ subroutine spmd_pack_doubles(inbuf, incount, outbuf, outsize, position, comm)
           integer, intent(in), optional :: comm
           double precision, dimension(incount), intent(in) :: inbuf
           integer, dimension(outsize), intent(inout) :: outbuf
-#ifdef MPI
           integer :: ierr, tag
+#ifdef MPI
           tag = 0  ! Default tag for pack operations
           call spmd_in(tag)
           if( present(comm) ) then

engine/source/mpi/spmd_recv.F90

@@ -87,7 +87,7 @@ subroutine spmd_recv_reals2D(buf, buf_count, source, tag, comm)
           implicit none
 #include "spmd.inc"
           integer, intent(in) :: buf_count, source, tag
-          real, dimension(buf_count,1), intent(inout) :: buf
+          real, dimension(1,buf_count), intent(inout) :: buf
           integer, intent(in), optional :: comm
 #ifdef MPI
           integer :: ierr
@@ -175,7 +175,7 @@ subroutine spmd_recv_doubles2D(buf, buf_count, source, tag,  comm)
           implicit none
 #include "spmd.inc"
           integer, intent(in) :: buf_count, source, tag
-          double precision, dimension(buf_count,1), intent(inout) :: buf
+          double precision, dimension(1,buf_count), intent(inout) :: buf
           integer, intent(in), optional :: comm
           integer :: ierr
 #ifdef MPI

engine/source/mpi/spmd_send.F90

@@ -88,7 +88,7 @@ subroutine spmd_send_reals2D(buf, buf_count, dest, tag,  comm)
 #include "spmd.inc"
           integer, intent(in) :: buf_count, dest, tag
           integer, intent(in), optional :: comm
-          real, dimension(1,buf_count), intent(in) :: buf
+          real, dimension(buf_count,1), intent(in) :: buf
           integer :: ierr
 #ifdef MPI
           call spmd_in(tag)
@@ -259,7 +259,7 @@ subroutine spmd_send_doubles2D(buf, buf_count, dest, tag, comm)
 #include "spmd.inc"
           integer, intent(in) :: buf_count, dest, tag
           integer, intent(in), optional :: comm
-          double precision, intent(in) :: buf(1,buf_count)
+          double precision, intent(in) :: buf(buf_count,1)
 #ifdef MPI
           integer :: ierr
           ! the MPI datatype for double precision is MPI_DOUBLE_PRECISION

engine/source/mpi/spmd_unpack.F90

@@ -116,8 +116,8 @@ subroutine spmd_unpack_doubles(inbuf, insize, position, outbuf, outcount, comm)
           integer, intent(in), optional :: comm
           integer, dimension(insize), intent(in) :: inbuf
           double precision, dimension(outcount), intent(out) :: outbuf
-#ifdef MPI
           integer :: ierr, tag
+#ifdef MPI
           tag = 0  ! Default tag for unpack operations
           call spmd_in(tag)
           if( present(comm) ) then

engine/source/mpi/spmd_wait.F90

@@ -99,7 +99,7 @@ subroutine spmd_wait(request, status)
           use spmd_error_mod, only: spmd_in, spmd_out
           implicit none
 #include "spmd.inc"
-          integer, intent(in) :: request
+          integer, intent(inout) :: request
           integer, dimension(MPI_STATUS_SIZE), optional, intent(inout) :: status
 #ifdef MPI
           integer :: ierr
@@ -113,9 +113,6 @@ subroutine spmd_wait(request, status)
 #endif
         end subroutine spmd_wait
 ! ======================================================================================================================
-
-! ======================================================================================================================
-! ======================================================================================================================
 !||====================================================================
 !||    spmd_waitany                  ../engine/source/mpi/spmd_wait.F90
 !||--- called by ------------------------------------------------------