[cig-commits] [commit] devel: Combine duplicates of sort_array_coordinates. (a918633)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu May 1 00:50:26 PDT 2014
Repository : ssh://geoshell/specfem3d
On branch : devel
Link : https://github.com/geodynamics/specfem3d/compare/cb32c88d6155d7974561a6f72fc17aea596e2c4d...50aa953c1db3f565d76415f5305410a529996b75
>---------------------------------------------------------------
commit a918633752b4a52942cbfee0aa423ebd6e002010
Author: Elliott Sales de Andrade <esalesde at physics.utoronto.ca>
Date: Thu Jan 23 16:42:16 2014 -0500
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.
>---------------------------------------------------------------
a918633752b4a52942cbfee0aa423ebd6e002010
src/auxiliaries/rules.mk | 1 +
src/decompose_mesh/fault_scotch.f90 | 168 ++--------------------------------
src/decompose_mesh/rules.mk | 1 +
src/generate_databases/get_MPI.f90 | 2 +-
src/meshfem3D/rules.mk | 1 +
src/shared/get_global.f90 | 168 ++--------------------------------
src/shared/sort_array_coordinates.f90 | 92 +++++++++----------
7 files changed, 57 insertions(+), 376 deletions(-)
diff --git a/src/auxiliaries/rules.mk b/src/auxiliaries/rules.mk
index dece7f3..647da76 100644
--- a/src/auxiliaries/rules.mk
+++ b/src/auxiliaries/rules.mk
@@ -204,6 +204,7 @@ create_movie_shakemap_AVS_DX_GMT_auxiliaries_OBJECTS = \
$O/param_reader.cc.o \
$O/read_parameter_file.shared.o \
$O/read_value_parameters.shared.o \
+ $O/sort_array_coordinates.shared.o \
$O/utm_geo.shared.o \
$(EMPTY_MACRO)
diff --git a/src/decompose_mesh/fault_scotch.f90 b/src/decompose_mesh/fault_scotch.f90
index 074dd01..fdeab7c 100644
--- a/src/decompose_mesh/fault_scotch.f90
+++ b/src/decompose_mesh/fault_scotch.f90
@@ -297,73 +297,20 @@ CONTAINS
double precision, intent(in) :: xyz_c(3,nspec)
double precision, dimension(nspec) :: work,xp,yp,zp
- integer, dimension(nspec) :: ind,ninseg,iwork
+ integer, dimension(nspec) :: ind,ninseg,iwork,ibool,iglob
+ integer :: nglob
logical :: ifseg(nspec)
- integer :: ispec,i,j
- integer :: nseg,ioff,iseg
- double precision :: SMALLVALTOL
+ double precision :: xtol
xp=xyz_c(1,:)
yp=xyz_c(2,:)
zp=xyz_c(3,:)
! define geometrical tolerance based upon typical size of the model
- SMALLVALTOL = 1.d-10 * maxval( maxval(xyz_c,2) - minval(xyz_c,2) )
+ xtol = 1.d-10 * maxval( maxval(xyz_c,2) - minval(xyz_c,2) )
- ! establish initial pointers
- do ispec=1,nspec
- locval(ispec)=ispec
- enddo
-
- ifseg(:)=.false.
-
- nseg=1
- ifseg(1)=.true.
- ninseg(1)=nspec
-
- 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,nspec
- if(dabs(xp(i)-xp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
- enddo
- else if(j == 2) then
- do i=2,nspec
- if(dabs(yp(i)-yp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
- enddo
- else
- do i=2,nspec
- if(dabs(zp(i)-zp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
- enddo
- endif
-
- ! count up number of different segments
- nseg=0
- do i=1,nspec
- if(ifseg(i)) then
- nseg=nseg+1
- ninseg(nseg)=1
- else
- ninseg(nseg)=ninseg(nseg)+1
- endif
- enddo
- enddo
+ call sort_array_coordinates(nspec,xp,yp,zp,ibool,iglob,locval,ifseg, &
+ nglob,ind,ninseg,iwork,work,xtol)
end subroutine lex_order
@@ -590,108 +537,5 @@ CONTAINS
end subroutine write_fault_database
-
-! -----------------------------------
-
-! sorting routines put in same file to allow for inlining
-
- subroutine rank(A,IND,N)
-!
-! Use Heap Sort (Numerical Recipes)
-!
- implicit none
-
- integer n
- double precision A(n)
- integer IND(n)
-
- 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
- 100 CONTINUE
- 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
- if (ir == 1) then
- ind(1)=indx
- return
- endif
- endif
- i=l
- j=l+l
- 200 CONTINUE
- IF (J <= IR) THEN
- 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
- goto 200
- endif
- IND(I)=INDX
- goto 100
-
- 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 IND(n)
- integer IA(n),IW(n)
- double precision A(n),B(n),C(n),W(n)
-
- 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
-
-! ------------------------------------------------------------------
-
-
end module fault_scotch
diff --git a/src/decompose_mesh/rules.mk b/src/decompose_mesh/rules.mk
index 9189c1a..344f7ed 100644
--- a/src/decompose_mesh/rules.mk
+++ b/src/decompose_mesh/rules.mk
@@ -53,6 +53,7 @@ decompose_mesh_SHARED_OBJECTS = \
$O/param_reader.cc.o \
$O/read_parameter_file.shared.o \
$O/read_value_parameters.shared.o \
+ $O/sort_array_coordinates.shared.o \
$(EMPTY_MACRO)
diff --git a/src/generate_databases/get_MPI.f90 b/src/generate_databases/get_MPI.f90
index f74827e..7bf1126 100644
--- a/src/generate_databases/get_MPI.f90
+++ b/src/generate_databases/get_MPI.f90
@@ -135,7 +135,7 @@
ibool_interfaces_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
reorder_interface_ext_mesh,locval,ifseg, &
nibool_interfaces_ext_mesh_true(iinterface), &
- ind_ext_mesh,ninseg_ext_mesh,iwork_ext_mesh,work_ext_mesh)
+ ind_ext_mesh,ninseg_ext_mesh,iwork_ext_mesh,work_ext_mesh,SMALLVAL_TOL)
! checks that number of MPI points are still the same
num_points1 = num_points1 + nibool_interfaces_ext_mesh(iinterface)
diff --git a/src/meshfem3D/rules.mk b/src/meshfem3D/rules.mk
index a9a8c99..c899b38 100644
--- a/src/meshfem3D/rules.mk
+++ b/src/meshfem3D/rules.mk
@@ -71,6 +71,7 @@ meshfem3D_SHARED_OBJECTS = \
$O/read_topo_bathy_file.shared.o \
$O/read_value_parameters.shared.o \
$O/safe_alloc_mod.shared.o \
+ $O/sort_array_coordinates.shared.o \
$O/utm_geo.shared.o \
$(EMPTY_MACRO)
diff --git a/src/shared/get_global.f90 b/src/shared/get_global.f90
index c89991f..7f6a4ab 100644
--- a/src/shared/get_global.f90
+++ b/src/shared/get_global.f90
@@ -48,7 +48,7 @@
integer ispec,i,j,ier
integer ieoff,ilocnum,nseg,ioff,iseg,ig
- integer, dimension(:), allocatable :: ind,ninseg,iwork
+ integer, dimension(:), allocatable :: ind,ninseg,iwork,idummy
double precision, dimension(:), allocatable :: work
! geometry tolerance parameter to calculate number of independent grid points
@@ -67,177 +67,21 @@
if( ier /= 0 ) stop 'error allocating array iwork'
allocate(work(npointot),stat=ier)
if( ier /= 0 ) stop 'error allocating array work'
+ allocate(idummy(npointot),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array idummy'
-! establish initial pointers
- do ilocnum=1,npointot
- locval(ilocnum) = ilocnum
- 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,SMALLVALTOL)
! deallocate arrays
deallocate(ind)
deallocate(ninseg)
deallocate(iwork)
deallocate(work)
+ deallocate(idummy)
end subroutine get_global
-! -----------------------------------
-
-! sorting routines put in same file to allow for inlining
-
- subroutine rank(A,IND,N)
-!
-! Use Heap Sort (Numerical Recipes)
-!
- implicit none
-
- integer n
- double precision A(n)
- integer IND(n)
-
- 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
- 100 CONTINUE
- 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
- if (ir == 1) then
- ind(1)=indx
- return
- endif
- endif
- i=l
- j=l+l
- 200 CONTINUE
- IF (J <= IR) THEN
- 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
- goto 200
- endif
- IND(I)=INDX
- goto 100
- 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 IND(n)
- integer IA(n),IW(n)
- double precision A(n),B(n),C(n),W(n)
-
- 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
-
! ------------------------------------------------------------------
diff --git a/src/shared/sort_array_coordinates.f90 b/src/shared/sort_array_coordinates.f90
index 25c5e75..a21e7f6 100644
--- a/src/shared/sort_array_coordinates.f90
+++ b/src/shared/sort_array_coordinates.f90
@@ -28,7 +28,7 @@
! subroutines to sort MPI buffers to assemble between chunks
subroutine sort_array_coordinates(npointot,x,y,z,ibool,iglob,locval,ifseg, &
- nglob,ind,ninseg,iwork,work)
+ nglob,ind,ninseg,iwork,work,xtol)
! this routine MUST be in double precision to avoid sensitivity
! to roundoff errors in the coordinates of the points
@@ -57,9 +57,6 @@
locval(ipoin)=ipoin
enddo
-! define a tolerance, normalized radius is 1., so let's use a small value
- xtol = SMALLVAL_TOL
-
ifseg(:)=.false.
nseg=1
@@ -67,55 +64,54 @@
ninseg(1)=npointot
do j=1,NDIM
-
! sort within each segment
- ioff=1
- do iseg=1,nseg
- if(j == 1) then
+ ioff=1
+ do iseg=1,nseg
+ if(j == 1) then
- call rank_buffers(x(ioff),ind,ninseg(iseg))
+ 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))
+ call rank_buffers(y(ioff),ind,ninseg(iseg))
- else
+ else
- call rank_buffers(z(ioff),ind,ninseg(iseg))
-
- endif
+ call rank_buffers(z(ioff),ind,ninseg(iseg))
- call swap_all_buffers(ibool(ioff),locval(ioff), &
- x(ioff),y(ioff),z(ioff),iwork,work,ind,ninseg(iseg))
+ endif
- ioff=ioff+ninseg(iseg)
- enddo
+ call swap_all_buffers(ibool(ioff),locval(ioff), &
+ x(ioff),y(ioff),z(ioff),iwork,work,ind,ninseg(iseg))
-! check for jumps in current coordinate
- if(j == 1) then
- do i=2,npointot
- if(dabs(x(i)-x(i-1)) > xtol) ifseg(i)=.true.
- enddo
- else if(j == 2) then
- do i=2,npointot
- if(dabs(y(i)-y(i-1)) > xtol) ifseg(i)=.true.
+ ioff=ioff+ninseg(iseg)
enddo
- else
- do i=2,npointot
- if(dabs(z(i)-z(i-1)) > xtol) 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
+ ! check for jumps in current coordinate
+ if (j == 1) then
+ do i=2,npointot
+ if(dabs(x(i)-x(i-1)) > xtol) ifseg(i)=.true.
+ enddo
+ else if (j == 2) then
+ do i=2,npointot
+ if(dabs(y(i)-y(i-1)) > xtol) ifseg(i)=.true.
+ enddo
else
- ninseg(nseg)=ninseg(nseg)+1
+ do i=2,npointot
+ if(dabs(z(i)-z(i-1)) > xtol) ifseg(i)=.true.
+ enddo
endif
- enddo
+
+ ! 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)
@@ -205,29 +201,23 @@
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