[cig-commits] [commit] devel, master: Combine duplicates of sort_array_coordinates. (2419de2)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:20:13 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
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