[cig-commits] [commit] devel: Combine duplicates of sort_array_coordinates. (2419de2)

cig_noreply at geodynamics.org cig_noreply at geodynamics.org
Wed Jun 4 13:02:28 PDT 2014


Repository : https://github.com/geodynamics/specfem3d_globe

On branch  : devel
Link       : https://github.com/geodynamics/specfem3d_globe/compare/754300e5a26e2a5d17069e661645fe819a7857a1...cd83b686a7c6014e132933f685157d43151a665e

>---------------------------------------------------------------

commit 2419de2008510cdbc6c400cf6039ed80ea5a43e0
Author: Elliott Sales de Andrade <esalesde at physics.utoronto.ca>
Date:   Sun Jun 1 01:35:03 2014 -0400

    Combine duplicates of sort_array_coordinates.
    
    This unfortunately causes a bit of extra allocation for unused
    parameters. Hopefully I can find a workaraound for that.


>---------------------------------------------------------------

2419de2008510cdbc6c400cf6039ed80ea5a43e0
 src/auxiliaries/rules.mk                           |   1 +
 src/meshfem3D/rules.mk                             |   2 +-
 src/shared/get_global.f90                          | 193 +--------------------
 src/shared/rules.mk                                |   1 +
 .../sort_array_coordinates.f90                     |  34 ++--
 5 files changed, 25 insertions(+), 206 deletions(-)

diff --git a/src/auxiliaries/rules.mk b/src/auxiliaries/rules.mk
index 583b616..e6ea1e6 100644
--- a/src/auxiliaries/rules.mk
+++ b/src/auxiliaries/rules.mk
@@ -77,6 +77,7 @@ auxiliaries_SHARED_OBJECTS = \
 	$O/read_value_parameters.shared.o \
 	$O/reduce.shared.o \
 	$O/rthetaphi_xyz.shared.o \
+	$O/sort_array_coordinates.shared.o \
 	$(EMPTY_MACRO)
 
 
diff --git a/src/meshfem3D/rules.mk b/src/meshfem3D/rules.mk
index 500fc69..a6d8f4c 100644
--- a/src/meshfem3D/rules.mk
+++ b/src/meshfem3D/rules.mk
@@ -107,7 +107,6 @@ meshfem3D_OBJECTS = \
 	$O/setup_inner_outer.check.o \
 	$O/setup_model.check.o \
 	$O/setup_MPI_interfaces.check.o \
-	$O/sort_array_coordinates.check.o \
 	$O/stretching_function.check.o \
 	$O/test_MPI_interfaces.check.o \
 	$O/write_AVS_DX_global_chunks_data.check.o \
@@ -185,6 +184,7 @@ meshfem3D_SHARED_OBJECTS = \
 	$O/reduce.shared.o \
 	$O/rthetaphi_xyz.shared.o \
 	$O/save_header_file.shared.o \
+	$O/sort_array_coordinates.shared.o \
 	$O/spline_routines.shared.o \
 	$O/write_VTK_file.shared.o \
 	$(EMPTY_MACRO)
diff --git a/src/shared/get_global.f90 b/src/shared/get_global.f90
index 1d3e294..dc4b018 100644
--- a/src/shared/get_global.f90
+++ b/src/shared/get_global.f90
@@ -49,9 +49,8 @@
 
   ! local variables
   double precision, dimension(:), allocatable :: work
-  integer, dimension(:), allocatable :: ind,ninseg,iwork
-  integer :: i,j,ier
-  integer :: ieoff,ilocnum,nseg,ioff,iseg,ig
+  integer, dimension(:), allocatable :: ind,ninseg,iwork,idummy
+  integer :: ier
 
   ! dynamically allocate arrays
   allocate(ind(npointot),stat=ier)
@@ -60,196 +59,20 @@
   allocate(ninseg(npointot),stat=ier)
   if (ier /= 0) stop 'error while allocating ninseg'
 
+  allocate(idummy(npointot),stat=ier)
+  if (ier /= 0) stop 'error while allocating idummy'
+
   allocate(iwork(npointot),stat=ier)
   if (ier /= 0) stop 'error while allocating iwork'
 
   allocate(work(npointot),stat=ier)
   if (ier /= 0) stop 'error while allocating work'
 
-  ! establish initial pointers
-  do i=1,npointot
-    locval(i) = i
-  enddo
-
-  ifseg(:) = .false.
-
-  nseg = 1
-  ifseg(1) = .true.
-  ninseg(1) = npointot
-
-  do j=1,NDIM
-    ! sort within each segment
-    ioff=1
-    do iseg=1,nseg
-      if (j == 1) then
-        call rank(xp(ioff),ind,ninseg(iseg))
-      else if (j == 2) then
-        call rank(yp(ioff),ind,ninseg(iseg))
-      else
-        call rank(zp(ioff),ind,ninseg(iseg))
-      endif
-
-      call swap_all(locval(ioff),xp(ioff),yp(ioff),zp(ioff),iwork,work,ind,ninseg(iseg))
-
-      ioff=ioff+ninseg(iseg)
-    enddo
-
-    ! check for jumps in current coordinate
-    ! compare the coordinates of the points within a small tolerance
-    if (j == 1) then
-      do i=2,npointot
-        if (dabs(xp(i)-xp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
-      enddo
-    else if (j == 2) then
-      do i=2,npointot
-        if (dabs(yp(i)-yp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
-      enddo
-    else
-      do i=2,npointot
-        if (dabs(zp(i)-zp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
-      enddo
-    endif
-
-    ! count up number of different segments
-    nseg=0
-    do i=1,npointot
-      if (ifseg(i)) then
-        nseg=nseg+1
-        ninseg(nseg)=1
-      else
-        ninseg(nseg)=ninseg(nseg)+1
-      endif
-    enddo
-  enddo
-
-  ! assign global node numbers (now sorted lexicographically)
-  ig=0
-  do i=1,npointot
-    if (ifseg(i)) ig=ig+1
-    iglob(locval(i))=ig
-  enddo
-
-  nglob=ig
+  call sort_array_coordinates(npointot,xp,yp,zp,idummy,iglob,locval,ifseg, &
+                              nglob,ind,ninseg,iwork,work)
 
   ! deallocate arrays
-  deallocate(ind,ninseg,iwork,work)
-
-! ------------------------------------------------------------------
-
-! get_global internal procedures follow
-
-! sorting routines put in same file to allow for inlining
-
-  contains
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
-  subroutine rank(A,IND,N)
-!
-! Use Heap Sort (Numerical Recipes)
-!
-  implicit none
-
-  integer :: n
-  double precision,dimension(n) :: A
-  integer,dimension(n) :: IND
-
-  ! local parameters
-  integer :: i,j,l,ir,indx
-  double precision :: q
-
-  do j=1,n
-   IND(j)=j
-  enddo
-
-  if (n == 1) return
-
-  L = n/2 + 1
-  ir = n
-
-  do while (.true.)
-
-    IF (l > 1) THEN
-      l = l-1
-      indx = ind(l)
-      q = a(indx)
-    ELSE
-      indx = ind(ir)
-      q = a(indx)
-      ind(ir) = ind(1)
-      ir = ir-1
-
-      ! checks exit criterion
-      if (ir == 1) then
-         ind(1) = indx
-         return
-      endif
-    endif
-
-    i = l
-    j = l+l
-
-    do while (J <= IR)
-      IF (J < IR) THEN
-        IF (A(IND(j)) < A(IND(j+1))) j=j+1
-      endif
-      IF (q < A(IND(j))) THEN
-        IND(I) = IND(J)
-        I = J
-        J = J+J
-      ELSE
-        J = IR+1
-      endif
-    enddo
-
-    IND(I)=INDX
-  enddo
-
-  end subroutine rank
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
-  subroutine swap_all(IA,A,B,C,IW,W,ind,n)
-!
-! swap arrays IA, A, B and C according to addressing in array IND
-!
-  implicit none
-
-  integer :: n
-  integer,dimension(n) :: IND
-  integer,dimension(n) :: IA,IW
-  double precision,dimension(n) :: A,B,C,W
-
-  ! local parameter
-  integer :: i
-
-  IW(:) = IA(:)
-  W(:) = A(:)
-
-  do i=1,n
-    IA(i) = IW(ind(i))
-    A(i) = W(ind(i))
-  enddo
-
-  W(:) = B(:)
-
-  do i=1,n
-    B(i) = W(ind(i))
-  enddo
-
-  W(:) = C(:)
-
-  do i=1,n
-    C(i) = W(ind(i))
-  enddo
-
-  end subroutine swap_all
-
-! ------------------------------------------------------------------
+  deallocate(ind,ninseg,idummy,iwork,work)
 
   end subroutine get_global
 
diff --git a/src/shared/rules.mk b/src/shared/rules.mk
index 2f89a97..d441108 100644
--- a/src/shared/rules.mk
+++ b/src/shared/rules.mk
@@ -66,6 +66,7 @@ shared_OBJECTS = \
 	$O/reduce.shared.o \
 	$O/rthetaphi_xyz.shared.o \
 	$O/save_header_file.shared.o \
+	$O/sort_array_coordinates.shared.o \
 	$O/spline_routines.shared.o \
 	$O/write_c_binary.cc.o \
 	$O/write_VTK_file.shared.o \
diff --git a/src/meshfem3D/sort_array_coordinates.f90 b/src/shared/sort_array_coordinates.f90
similarity index 91%
rename from src/meshfem3D/sort_array_coordinates.f90
rename to src/shared/sort_array_coordinates.f90
index 864c003..fb6521b 100644
--- a/src/meshfem3D/sort_array_coordinates.f90
+++ b/src/shared/sort_array_coordinates.f90
@@ -69,9 +69,9 @@
     ! sort within each segment
     ioff=1
     do iseg=1,nseg
-      if(j == 1) then
+      if (j == 1) then
         call rank_buffers(x(ioff),ind,ninseg(iseg))
-      else if(j == 2) then
+      else if (j == 2) then
         call rank_buffers(y(ioff),ind,ninseg(iseg))
       else
         call rank_buffers(z(ioff),ind,ninseg(iseg))
@@ -85,24 +85,24 @@
 
     ! check for jumps in current coordinate
     ! define a tolerance, normalized radius is 1., so let's use a small value
-    if(j == 1) then
+    if (j == 1) then
       do i=2,npointot
-        if(dabs(x(i)-x(i-1)) > SMALLVALTOL ) ifseg(i)=.true.
+        if (dabs(x(i)-x(i-1)) > SMALLVALTOL) ifseg(i)=.true.
       enddo
-    else if(j == 2) then
+    else if (j == 2) then
       do i=2,npointot
-        if(dabs(y(i)-y(i-1)) > SMALLVALTOL ) ifseg(i)=.true.
+        if (dabs(y(i)-y(i-1)) > SMALLVALTOL) ifseg(i)=.true.
       enddo
     else
       do i=2,npointot
-        if(dabs(z(i)-z(i-1)) > SMALLVALTOL ) ifseg(i)=.true.
+        if (dabs(z(i)-z(i-1)) > SMALLVALTOL) ifseg(i)=.true.
       enddo
     endif
 
     ! count up number of different segments
     nseg=0
     do i=1,npointot
-      if(ifseg(i)) then
+      if (ifseg(i)) then
         nseg=nseg+1
         ninseg(nseg)=1
       else
@@ -115,7 +115,7 @@
   ! assign global node numbers (now sorted lexicographically)
   ig=0
   do i=1,npointot
-    if(ifseg(i)) ig=ig+1
+    if (ifseg(i)) ig=ig+1
     iglob(locval(i))=ig
   enddo
 
@@ -208,29 +208,23 @@
   ! local parameter
   integer :: i
 
-  do i=1,n
-    W(i)=A(i)
-    IW(i)=IA(i)
-  enddo
+  W(:) = A(:)
+  IW(:) = IA(:)
 
   do i=1,n
     A(i)=W(ind(i))
     IA(i)=IW(ind(i))
   enddo
 
-  do i=1,n
-    W(i)=B(i)
-    IW(i)=IB(i)
-  enddo
+  W(:) = B(:)
+  IW(:) = IB(:)
 
   do i=1,n
     B(i)=W(ind(i))
     IB(i)=IW(ind(i))
   enddo
 
-  do i=1,n
-    W(i)=C(i)
-  enddo
+  W(:) = C(:)
 
   do i=1,n
     C(i)=W(ind(i))



More information about the CIG-COMMITS mailing list