From 4eae2b3e7e25a211904722cae3bb52ae5b22456e Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Mon, 19 Feb 2018 08:43:52 +0000 Subject: [PATCH] Initial steps for long indices. --- test/idx/Makefile | 36 ++++++++ test/idx/tryidxijk.f90 | 19 +++++ util/psb_partidx_mod.F90 | 176 ++++++++++++++++++++++++++++++++++++++- 3 files changed, 229 insertions(+), 2 deletions(-) create mode 100644 test/idx/Makefile create mode 100644 test/idx/tryidxijk.f90 diff --git a/test/idx/Makefile b/test/idx/Makefile new file mode 100644 index 00000000..0ecc3634 --- /dev/null +++ b/test/idx/Makefile @@ -0,0 +1,36 @@ +INSTALLDIR=../.. +INCDIR=$(INSTALLDIR)/include +MODDIR=$(INSTALLDIR)/modules/ +include $(INCDIR)/Make.inc.psblas +# +# Libraries used +LIBDIR=$(INSTALLDIR)/lib +PSBLAS_LIB= -L$(LIBDIR) -lpsb_util -lpsb_krylov -lpsb_prec -lpsb_base +LDLIBS=$(PSBLDLIBS) +# +# Compilers and such +# +CCOPT= -g +FINCLUDES=$(FMFLAG)$(MODDIR) $(FMFLAG). + + +EXEDIR=./runs + +all: tryidxijk + +tryidxijk: tryidxijk.o + $(FLINK) tryidxijk.o -o tryidxijk $(PSBLAS_LIB) $(LDLIBS) + /bin/mv tryidxijk $(EXEDIR) + + + +clean: + /bin/rm -f tryidxijk.o *$(.mod) \ + $(EXEDIR)/tryidxijk +verycleanlib: + (cd ../..; make veryclean) +lib: + (cd ../../; make library) + + + diff --git a/test/idx/tryidxijk.f90 b/test/idx/tryidxijk.f90 new file mode 100644 index 00000000..27e1c7c4 --- /dev/null +++ b/test/idx/tryidxijk.f90 @@ -0,0 +1,19 @@ +program tryidxijk + use psb_base_mod + use psb_util_mod + + integer(psb_lpk_) :: idx,idxm + integer(psb_ipk_) :: nx,ny,nz + integer(psb_ipk_) :: i,j,k, sidx + + idxm = 1000 + idxm = idxm*2000*1000 + nx = 2000 + ny = 2000 + nz = 2000 + do idx = idxm+300*1000*1000, idxm+300*1000*1000+50000 + call idx2ijk(i,j,k,idx,nx,ny,nz) + sidx = idx + write(*,*) 'idx2ijk: ',idx,i,j,k, sidx + end do +end program tryidxijk diff --git a/util/psb_partidx_mod.F90 b/util/psb_partidx_mod.F90 index 9b294e8f..433f87fb 100644 --- a/util/psb_partidx_mod.F90 +++ b/util/psb_partidx_mod.F90 @@ -41,11 +41,15 @@ module psb_partidx_mod use psb_base_mod, only : psb_ipk_ interface idx2ijk - module procedure idx2ijk3d, idx2ijkv, idx2ijk2d + module procedure idx2ijk3d, idx2ijkv, idx2ijk2d,& + & lidx2ijk3d, lidx2ijkv, lidx2ijk2d,& + & lidx2lijk3d, lidx2lijkv, lidx2lijk2d end interface idx2ijk interface ijk2idx - module procedure ijk2idx3d, ijk2idxv, ijk2idx2d + module procedure ijk2idx3d, ijk2idxv, ijk2idx2d!,& +!!$ & ijk2idx3d, ijk2idxv, ijk2lidx2d,& +!!$ & lijk2lidx3d, lijk2lidxv, lijk2lidx2d end interface ijk2idx @@ -141,6 +145,174 @@ contains end do end subroutine idx2ijkv + + subroutine lidx2ijk3d(i,j,k,idx,nx,ny,nz,base) + use psb_base_mod, only : psb_ipk_, psb_lpk_ + implicit none + integer(psb_ipk_), intent(out) :: i,j,k + integer(psb_lpk_), intent(in) :: idx + integer(psb_ipk_), intent(in) :: nx,ny,nz + integer(psb_ipk_), intent(in), optional :: base + + integer(psb_ipk_) :: coords(3) + + call idx2ijk(coords,idx,[nx,ny,nz],base) + + k = coords(3) + j = coords(2) + i = coords(1) + + end subroutine lidx2ijk3d + + subroutine lidx2ijk2d(i,j,idx,nx,ny,base) + use psb_base_mod, only : psb_ipk_, psb_lpk_ + implicit none + integer(psb_ipk_), intent(out) :: i,j + integer(psb_lpk_), intent(in) :: idx + integer(psb_ipk_), intent(in) :: nx,ny + integer(psb_ipk_), intent(in), optional :: base + + integer(psb_ipk_) :: coords(2) + + call idx2ijk(coords,idx,[nx,ny],base) + + j = coords(2) + i = coords(1) + + end subroutine lidx2ijk2d + + ! + ! Given a global index IDX and the domain size (NX,NY,NZ) + ! compute the point coordinates (I,J,K) + ! Optional argument: base 0 or 1, default 1 + ! + ! This mapping is equivalent to a loop nesting: + ! idx = base + ! do i=1,nx + ! do j=1,ny + ! do k=1,nz + ! ijk2idx(i,j,k) = idx + ! idx = idx + 1 + subroutine lidx2ijkv(coords,idx,dims,base) + use psb_base_mod, only : psb_ipk_, psb_lpk_ + implicit none + integer(psb_ipk_), intent(out) :: coords(:) + integer(psb_lpk_), intent(in) :: idx + integer(psb_ipk_), intent(in) :: dims(:) + integer(psb_ipk_), intent(in), optional :: base + + integer(psb_lpk_) :: base_, idx_ + integer(psb_ipk_) :: i, sz + if (present(base)) then + base_ = base + else + base_ = 1 + end if + + idx_ = idx - base_ + + if (size(coords) < size(dims)) then + write(0,*) 'Error: size mismatch ',size(coords),size(dims) + coords = 0 + return + end if + + ! + ! This code is equivalent to (3D case) + ! k = mod(idx_,nz) + base_ + ! j = mod(idx_/nz,ny) + base_ + ! i = mod(idx_/(nx*ny),nx) + base_ + ! + do i=size(dims),1,-1 + coords(i) = mod(idx_,dims(i)) + base_ + idx_ = idx_ / dims(i) + end do + + end subroutine lidx2ijkv + + subroutine lidx2lijk3d(i,j,k,idx,nx,ny,nz,base) + use psb_base_mod, only : psb_ipk_, psb_lpk_ + implicit none + integer(psb_lpk_), intent(out) :: i,j,k + integer(psb_lpk_), intent(in) :: idx + integer(psb_lpk_), intent(in) :: nx,ny,nz + integer(psb_ipk_), intent(in), optional :: base + + integer(psb_lpk_) :: coords(3) + + call idx2ijk(coords,idx,[nx,ny,nz],base) + + k = coords(3) + j = coords(2) + i = coords(1) + + end subroutine lidx2lijk3d + + subroutine lidx2lijk2d(i,j,idx,nx,ny,base) + use psb_base_mod, only : psb_ipk_, psb_lpk_ + implicit none + integer(psb_lpk_), intent(out) :: i,j + integer(psb_lpk_), intent(in) :: idx + integer(psb_lpk_), intent(in) :: nx,ny + integer(psb_ipk_), intent(in), optional :: base + + integer(psb_lpk_) :: coords(2) + + call idx2ijk(coords,idx,[nx,ny],base) + + j = coords(2) + i = coords(1) + + end subroutine lidx2lijk2d + + ! + ! Given a global index IDX and the domain size (NX,NY,NZ) + ! compute the point coordinates (I,J,K) + ! Optional argument: base 0 or 1, default 1 + ! + ! This mapping is equivalent to a loop nesting: + ! idx = base + ! do i=1,nx + ! do j=1,ny + ! do k=1,nz + ! ijk2idx(i,j,k) = idx + ! idx = idx + 1 + subroutine lidx2lijkv(coords,idx,dims,base) + use psb_base_mod, only : psb_ipk_, psb_lpk_ + implicit none + integer(psb_lpk_), intent(out) :: coords(:) + integer(psb_lpk_), intent(in) :: idx + integer(psb_lpk_), intent(in) :: dims(:) + integer(psb_ipk_), intent(in), optional :: base + + integer(psb_lpk_) :: base_, idx_ + integer(psb_lpk_) :: i, sz + if (present(base)) then + base_ = base + else + base_ = 1 + end if + + idx_ = idx - base_ + + if (size(coords) < size(dims)) then + write(0,*) 'Error: size mismatch ',size(coords),size(dims) + coords = 0 + return + end if + + ! + ! This code is equivalent to (3D case) + ! k = mod(idx_,nz) + base_ + ! j = mod(idx_/nz,ny) + base_ + ! i = mod(idx_/(nx*ny),nx) + base_ + ! + do i=size(dims),1,-1 + coords(i) = mod(idx_,dims(i)) + base_ + idx_ = idx_ / dims(i) + end do + + end subroutine lidx2lijkv ! ! Given a triple (I,J,K) and the domain size (NX,NY,NZ)