fdurastante/psblas3
1
0
Fork 0
mirror of https://github.com/sfilippone/psblas3.git synced 2026-10-16 10:52:08 +00:00
Code Issues Packages Projects Releases Wiki Activity

New psb_partidx_mod.F90.

Browse Source 
Merged into psb_util.
Updated all test/pargen progs.
This commit is contained in:
Salvatore Filippone
2018-02-18 09:42:35 +00:00
parent 470c6658f9
commit aaaf4c9f09
10 changed files with 773 additions and 322 deletions
Download Patch File Download Diff File Expand all files Collapse all files
base/modules/psb_const_mod.F90
+2
View File
@@ -41,6 +41,7 @@ module psb_const_mod
#endif
! This is always an 8-byte integer.
integer, parameter :: psb_long_int_k_ = int64
integer, parameter :: psb_lpk_ = psb_long_int_k_
! This is always a 4-byte integer, for MPI-related stuff
integer, parameter :: psb_mpik_ = int32
!
@@ -60,6 +61,7 @@ module psb_const_mod
! This is always an 8-byte integer.
integer, parameter :: longndig=12
integer, parameter :: psb_long_int_k_ = selected_int_kind(longndig)
integer, parameter :: psb_lpk_ = psb_long_int_k_
! This is always a 4-byte integer, for MPI-related stuff
integer, parameter :: psb_mpik_ = kind(1)
!
test/pargen/Makefile
+3 -3
View File
@@ -18,8 +18,8 @@ EXEDIR=./runs
all: psb_d_pde3d psb_s_pde3d psb_d_pde2d psb_s_pde2d
tryidxijk: tryidxijk.o
$(FLINK) tryidxijk.o -o tryidxijk $(PSBLAS_LIB) $(LDLIBS)
tryidxijk: tryidxijk.o
$(FLINK) tryidxijk.o -o tryidxijk $(PSBLAS_LIB) $(LDLIBS)
/bin/mv tryidxijk $(EXEDIR)
psb_d_pde3d: psb_d_pde3d.o
@@ -42,7 +42,7 @@ psb_s_pde2d: psb_s_pde2d.o
clean:
/bin/rm -f psb_d_pde3d.o psb_s_pde3d.o psb_d_pde2d.o psb_s_pde2d.o *$(.mod) \
/bin/rm -f psb_d_pde3d.o psb_s_pde3d.o psb_d_pde2d.o psb_s_pde2d.o partidx.o *$(.mod) \
$(EXEDIR)/psb_d_pde3d $(EXEDIR)/psb_s_pde3d $(EXEDIR)/psb_d_pde2d $(EXEDIR)/psb_s_pde2d
verycleanlib:
(cd ../..; make veryclean)
test/pargen/psb_d_pde2d.f90
+113 -29
View File
@@ -50,10 +50,12 @@
!
! Note that if b1=b2=c=0., the PDE is the Laplace equation.
!
! In this sample program the index space of the discretized
! computational domain is first numbered sequentially in a standard way,
! then the corresponding vector is distributed according to a BLOCK
! data distribution.
! There are three choices available for data distribution:
! 1. A simple BLOCK distribution
! 2. A ditribution based on arbitrary assignment of indices to processes,
! typically from a graph partitioner
! 3. A 2D distribution in which the unit square is partitioned
! into rectangles, each one assigned to a process.
!
module psb_d_pde2d_mod
@@ -91,8 +93,9 @@ contains
! the rhs.
!
subroutine psb_d_gen_pde2d(ictxt,idim,a,bv,xv,desc_a,afmt,&
& a1,a2,b1,b2,c,g,info,f,amold,vmold,imold,nrl,iv)
& a1,a2,b1,b2,c,g,info,f,amold,vmold,imold,partition,nrl,iv)
use psb_base_mod
use psb_util_mod
!
! Discretizes the partial differential equation
!
@@ -120,7 +123,7 @@ contains
class(psb_d_base_sparse_mat), optional :: amold
class(psb_d_base_vect_type), optional :: vmold
class(psb_i_base_vect_type), optional :: imold
integer(psb_ipk_), optional :: nrl, iv(:)
integer(psb_ipk_), optional :: partition, nrl,iv(:)
! Local variables.
@@ -129,9 +132,13 @@ contains
type(psb_d_coo_sparse_mat) :: acoo
type(psb_d_csr_sparse_mat) :: acsr
real(psb_dpk_) :: zt(nb),x,y,z
integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k
integer(psb_ipk_) :: m,n,nnz,nr,nt,glob_row,nlr,i,j,ii,ib,k, partition_
integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
integer(psb_ipk_) :: np, iam, nr, nt
! For 2D partition
integer(psb_ipk_) :: npx,npy,npdims(2),iamx,iamy,mynx,myny
integer(psb_ipk_), allocatable :: bndx(:),bndy(:)
! Process grid
integer(psb_ipk_) :: np, iam
integer(psb_ipk_) :: icoeff
integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:)
real(psb_dpk_), allocatable :: val(:)
@@ -161,6 +168,17 @@ contains
sqdeltah = deltah*deltah
deltah2 = 2.d0* deltah
if (present(partition)) then
if ((1<= partition).and.(partition <= 3)) then
partition_ = partition
else
write(*,*) 'Invalid partition choice ',partition,' defaulting to 3'
partition_ = 3
end if
else
partition_ = 3
end if
! initialize array descriptor and sparse matrix storage. provide an
! estimate of the number of non zeroes
@@ -169,7 +187,9 @@ contains
nnz = ((n*7)/(np))
if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n
if (.not.present(iv)) then
select case(partition_)
case(1)
! A BLOCK partition
if (present(nrl)) then
nr = nrl
else
@@ -189,27 +209,99 @@ contains
call psb_abort(ictxt)
return
end if
else
if (size(iv) /= m) then
write(psb_err_unit,*) iam, 'Initialization error IV',size(iv),m
!
! First example of use of CDALL: specify for each process a number of
! contiguous rows
!
call psb_cdall(ictxt,desc_a,info,nl=nr)
myidx = desc_a%get_global_indices()
nlr = size(myidx)
case(2)
! A partition defined by the user through IV
if (present(iv)) then
if (size(iv) /= m) then
write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
return
end if
else
write(psb_err_unit,*) iam, 'Initialization error: IV not present'
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
return
end if
end if
call psb_barrier(ictxt)
t0 = psb_wtime()
if (present(iv)) then
!
! Second example of use of CDALL: specify for each row the
! process that owns it
!
call psb_cdall(ictxt,desc_a,info,vg=iv)
else
call psb_cdall(ictxt,desc_a,info,nl=nr)
end if
myidx = desc_a%get_global_indices()
nlr = size(myidx)
case(3)
! A 2-dimensional partition
! A nifty MPI function will split the process list
npdims = 0
call mpi_dims_create(np,2,npdims,info)
npx = npdims(1)
npy = npdims(2)
allocate(bndx(0:npx),bndy(0:npy))
! We can reuse idx2ijk for process indices as well.
call idx2ijk(iamx,iamy,iam,npx,npy,base=0)
! Now let's split the 2D square in rectangles
call dist1Didx(bndx,idim,npx)
mynx = bndx(iamx+1)-bndx(iamx)
call dist1Didx(bndy,idim,npy)
myny = bndy(iamy+1)-bndy(iamy)
! How many indices do I own?
nlr = mynx*myny
allocate(myidx(nlr))
! Now, let's generate the list of indices I own
nr = 0
do i=bndx(iamx),bndx(iamx+1)-1
do j=bndy(iamy),bndx(iamy+1)-1
nr = nr + 1
call ijk2idx(myidx(nr),i,j,idim,idim)
end do
end do
if (nr /= nlr) then
write(psb_err_unit,*) iam,iamx,iamy, 'Initialization error: NR vs NLR ',&
& nr,nlr,mynx,myny
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
end if
!
! Third example of use of CDALL: specify for each process
! the set of global indices it owns.
!
call psb_cdall(ictxt,desc_a,info,vl=myidx)
case default
write(psb_err_unit,*) iam, 'Initialization error: should not get here'
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
return
end select
if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz)
! define rhs from boundary conditions; also build initial guess
if (info == psb_success_) call psb_geall(xv,desc_a,info)
if (info == psb_success_) call psb_geall(bv,desc_a,info)
call psb_barrier(ictxt)
talc = psb_wtime()-t0
@@ -233,9 +325,6 @@ contains
endif
myidx = desc_a%get_global_indices()
nlr = size(myidx)
! loop over rows belonging to current process in a block
! distribution.
@@ -249,13 +338,8 @@ contains
! local matrix pointer
glob_row=myidx(i)
! compute gridpoint coordinates
if (mod(glob_row,(idim)) == 0) then
ix = glob_row/(idim)
else
ix = glob_row/(idim)+1
endif
iy = (glob_row-(ix-1)*idim)
! x, y
call idx2ijk(ix,iy,glob_row,idim,idim)
! x, y coordinates
x = (ix-1)*deltah
y = (iy-1)*deltah
test/pargen/psb_d_pde3d.f90
+115 -100
View File
@@ -50,13 +50,14 @@
!
! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation.
!
! There are two choices available for data distribution:
! There are three choices available for data distribution:
! 1. A simple BLOCK distribution
! 2. A 3D distribution in which the unit cube is partitioned
! 2. A ditribution based on arbitrary assignment of indices to processes,
! typically from a graph partitioner
! 3. A 3D distribution in which the unit cube is partitioned
! into subcubes, each one assigned to a process.
!
!
!
module psb_d_pde3d_mod
@@ -76,9 +77,6 @@ module psb_d_pde3d_mod
module procedure psb_d_gen_pde3d
end interface psb_gen_pde3d
integer, private :: dims(3),coords(3)
integer, private, allocatable :: rk2coo(:,:), coo2rk(:,:,:)
contains
function d_null_func_3d(x,y,z) result(val)
@@ -90,78 +88,15 @@ contains
end function d_null_func_3d
!
! 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 idx2ijk(i,j,k,idx,nx,ny,nz,base)
integer(psb_ipk_), intent(out) :: i,j,k
integer(psb_ipk_), intent(in) :: idx,nx,ny,nz
integer(psb_ipk_), intent(in), optional :: base
integer(psb_ipk) :: base_, idx_
if (present(base)) then
base_ = base
else
base_ = 1
end if
idx_ = idx - base_
i = idx/(nx*ny)
j = (idx - i*nx*ny)/ny
k = (idx - i*nx*ny -j*ny)/nz
i = i + base_
j = j + base_
k = k + base_
end subroutine idx2ijk
!
! Given a triple (I,J,K) and the domain size (NX,NY,NZ)
! compute the global index IDX
! 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 ijk2idx(idx,i,j,k,nx,ny,nz,base)
integer(psb_ipk_), intent(out) :: idx,
integer(psb_ipk_), intent(in) :: i,j,k,nx,ny,nz
integer(psb_ipk_), intent(in), optional :: base
integer(psb_ipk) :: base_
if (present(base)) then
base_ = base
else
base_ = 1
end if
idx = ((i-base_)*nx*ny + (j-base_)*nx + k -base_) + base_
end subroutine ijk2idx
!
! subroutine to allocate and fill in the coefficient matrix and
! the rhs.
!
subroutine psb_d_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,&
& a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,imold,nrl,iv)
& a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,imold,partition,nrl,iv)
use psb_base_mod
use psb_util_mod
!
! Discretizes the partial differential equation
!
@@ -189,7 +124,7 @@ contains
class(psb_d_base_sparse_mat), optional :: amold
class(psb_d_base_vect_type), optional :: vmold
class(psb_i_base_vect_type), optional :: imold
integer(psb_ipk_), optional :: nrl,iv(:)
integer(psb_ipk_), optional :: partition, nrl,iv(:)
! Local variables.
@@ -198,9 +133,13 @@ contains
type(psb_d_coo_sparse_mat) :: acoo
type(psb_d_csr_sparse_mat) :: acsr
real(psb_dpk_) :: zt(nb),x,y,z
integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k
integer(psb_ipk_) :: m,n,nnz,nr,nt,glob_row,nlr,i,j,ii,ib,k, partition_
integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
integer(psb_ipk_) :: np, iam, nr, nt
! For 3D partition
integer(psb_ipk_) :: npx,npy,npz, npdims(3),iamx,iamy,iamz,mynx,myny,mynz
integer(psb_ipk_), allocatable :: bndx(:),bndy(:),bndz(:)
! Process grid
integer(psb_ipk_) :: np, iam
integer(psb_ipk_) :: icoeff
integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:)
real(psb_dpk_), allocatable :: val(:)
@@ -230,15 +169,28 @@ contains
sqdeltah = deltah*deltah
deltah2 = 2.d0* deltah
if (present(partition)) then
if ((1<= partition).and.(partition <= 3)) then
partition_ = partition
else
write(*,*) 'Invalid partition choice ',partition,' defaulting to 3'
partition_ = 3
end if
else
partition_ = 3
end if
! initialize array descriptor and sparse matrix storage. provide an
! estimate of the number of non zeroes
m = idim*idim*idim
n = m
nnz = ((n*9)/(np))
if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n
if (.not.present(iv)) then
select case(partition_)
case(1)
! A BLOCK partition
if (present(nrl)) then
nr = nrl
else
@@ -258,24 +210,99 @@ contains
call psb_abort(ictxt)
return
end if
else
if (size(iv) /= m) then
write(psb_err_unit,*) iam, 'Initialization error IV',size(iv),m
!
! First example of use of CDALL: specify for each process a number of
! contiguous rows
!
call psb_cdall(ictxt,desc_a,info,nl=nr)
myidx = desc_a%get_global_indices()
nlr = size(myidx)
case(2)
! A partition defined by the user through IV
if (present(iv)) then
if (size(iv) /= m) then
write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
return
end if
else
write(psb_err_unit,*) iam, 'Initialization error: IV not present'
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
return
end if
end if
call psb_barrier(ictxt)
t0 = psb_wtime()
if (present(iv)) then
!
! Second example of use of CDALL: specify for each row the
! process that owns it
!
call psb_cdall(ictxt,desc_a,info,vg=iv)
else
call psb_cdall(ictxt,desc_a,info,nl=nr)
end if
myidx = desc_a%get_global_indices()
nlr = size(myidx)
case(3)
! A 3-dimensional partition
! A nifty MPI function will split the process list
npdims = 0
call mpi_dims_create(np,3,npdims,info)
npx = npdims(1)
npy = npdims(2)
npz = npdims(3)
allocate(bndx(0:npx),bndy(0:npy),bndz(0:npz))
! We can reuse idx2ijk for process indices as well.
call idx2ijk(iamx,iamy,iamz,iam,npx,npy,npz,base=0)
! Now let's split the 3D cube in hexahedra
call dist1Didx(bndx,idim,npx)
mynx = bndx(iamx+1)-bndx(iamx)
call dist1Didx(bndy,idim,npy)
myny = bndy(iamy+1)-bndy(iamy)
call dist1Didx(bndz,idim,npz)
mynz = bndz(iamz+1)-bndz(iamz)
! How many indices do I own?
nlr = mynx*myny*mynz
allocate(myidx(nlr))
! Now, let's generate the list of indices I own
nr = 0
do i=bndx(iamx),bndx(iamx+1)-1
do j=bndy(iamy),bndx(iamy+1)-1
do k=bndz(iamz),bndz(iamz+1)-1
nr = nr + 1
call ijk2idx(myidx(nr),i,j,k,idim,idim,idim)
end do
end do
end do
if (nr /= nlr) then
write(psb_err_unit,*) iam,iamx,iamy,iamz, 'Initialization error: NR vs NLR ',&
& nr,nlr,mynx,myny,mynz
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
end if
!
! Third example of use of CDALL: specify for each process
! the set of global indices it owns.
!
call psb_cdall(ictxt,desc_a,info,vl=myidx)
case default
write(psb_err_unit,*) iam, 'Initialization error: should not get here'
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
return
end select
if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz)
! define rhs from boundary conditions; also build initial guess
if (info == psb_success_) call psb_geall(xv,desc_a,info)
@@ -303,8 +330,6 @@ contains
goto 9999
endif
myidx = desc_a%get_global_indices()
nlr = size(myidx)
! loop over rows belonging to current process in a block
! distribution.
@@ -319,18 +344,8 @@ contains
! local matrix pointer
glob_row=myidx(i)
! compute gridpoint coordinates
if (mod(glob_row,(idim*idim)) == 0) then
ix = glob_row/(idim*idim)
else
ix = glob_row/(idim*idim)+1
endif
if (mod((glob_row-(ix-1)*idim*idim),idim) == 0) then
iy = (glob_row-(ix-1)*idim*idim)/idim
else
iy = (glob_row-(ix-1)*idim*idim)/idim+1
endif
iz = glob_row-(ix-1)*idim*idim-(iy-1)*idim
! x, y, x coordinates
call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
! x, y, z coordinates
x = (ix-1)*deltah
y = (iy-1)*deltah
z = (iz-1)*deltah
test/pargen/psb_s_pde2d.f90
+135 -36
View File
@@ -50,10 +50,12 @@
!
! Note that if b1=b2=c=0., the PDE is the Laplace equation.
!
! In this sample program the index space of the discretized
! computational domain is first numbered sequentially in a standard way,
! then the corresponding vector is distributed according to a BLOCK
! data distribution.
! There are three choices available for data distribution:
! 1. A simple BLOCK distribution
! 2. A ditribution based on arbitrary assignment of indices to processes,
! typically from a graph partitioner
! 3. A 2D distribution in which the unit square is partitioned
! into rectangles, each one assigned to a process.
!
module psb_s_pde2d_mod
@@ -91,8 +93,9 @@ contains
! the rhs.
!
subroutine psb_s_gen_pde2d(ictxt,idim,a,bv,xv,desc_a,afmt,&
& a1,a2,b1,b2,c,g,info,f,amold,vmold,imold,nrl)
& a1,a2,b1,b2,c,g,info,f,amold,vmold,imold,partition,nrl,iv)
use psb_base_mod
use psb_util_mod
!
! Discretizes the partial differential equation
!
@@ -120,7 +123,7 @@ contains
class(psb_s_base_sparse_mat), optional :: amold
class(psb_s_base_vect_type), optional :: vmold
class(psb_i_base_vect_type), optional :: imold
integer(psb_ipk_), optional :: nrl
integer(psb_ipk_), optional :: partition, nrl,iv(:)
! Local variables.
@@ -129,9 +132,13 @@ contains
type(psb_s_coo_sparse_mat) :: acoo
type(psb_s_csr_sparse_mat) :: acsr
real(psb_spk_) :: zt(nb),x,y,z
integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k
integer(psb_ipk_) :: m,n,nnz,nr,nt,glob_row,nlr,i,j,ii,ib,k, partition_
integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
integer(psb_ipk_) :: np, iam, nr, nt
! For 2D partition
integer(psb_ipk_) :: npx,npy,npdims(2),iamx,iamy,mynx,myny
integer(psb_ipk_), allocatable :: bndx(:),bndy(:)
! Process grid
integer(psb_ipk_) :: np, iam
integer(psb_ipk_) :: icoeff
integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:)
real(psb_spk_), allocatable :: val(:)
@@ -161,6 +168,17 @@ contains
sqdeltah = deltah*deltah
deltah2 = 2.e0* deltah
if (present(partition)) then
if ((1<= partition).and.(partition <= 3)) then
partition_ = partition
else
write(*,*) 'Invalid partition choice ',partition,' defaulting to 3'
partition_ = 3
end if
else
partition_ = 3
end if
! initialize array descriptor and sparse matrix storage. provide an
! estimate of the number of non zeroes
@@ -169,32 +187,121 @@ contains
nnz = ((n*7)/(np))
if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n
if (present(nrl)) then
nr = nrl
else
!
! Using a simple BLOCK distribution.
!
nt = (m+np-1)/np
nr = max(0,min(nt,m-(iam*nt)))
end if
select case(partition_)
case(1)
! A BLOCK partition
if (present(nrl)) then
nr = nrl
else
!
! Using a simple BLOCK distribution.
!
nt = (m+np-1)/np
nr = max(0,min(nt,m-(iam*nt)))
end if
nt = nr
call psb_sum(ictxt,nt)
if (nt /= m) then
write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m
nt = nr
call psb_sum(ictxt,nt)
if (nt /= m) then
write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
return
end if
!
! First example of use of CDALL: specify for each process a number of
! contiguous rows
!
call psb_cdall(ictxt,desc_a,info,nl=nr)
myidx = desc_a%get_global_indices()
nlr = size(myidx)
case(2)
! A partition defined by the user through IV
if (present(iv)) then
if (size(iv) /= m) then
write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
return
end if
else
write(psb_err_unit,*) iam, 'Initialization error: IV not present'
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
return
end if
!
! Second example of use of CDALL: specify for each row the
! process that owns it
!
call psb_cdall(ictxt,desc_a,info,vg=iv)
myidx = desc_a%get_global_indices()
nlr = size(myidx)
case(3)
! A 2-dimensional partition
! A nifty MPI function will split the process list
npdims = 0
call mpi_dims_create(np,2,npdims,info)
npx = npdims(1)
npy = npdims(2)
allocate(bndx(0:npx),bndy(0:npy))
! We can reuse idx2ijk for process indices as well.
call idx2ijk(iamx,iamy,iam,npx,npy,base=0)
! Now let's split the 2D square in rectangles
call dist1Didx(bndx,idim,npx)
mynx = bndx(iamx+1)-bndx(iamx)
call dist1Didx(bndy,idim,npy)
myny = bndy(iamy+1)-bndy(iamy)
! How many indices do I own?
nlr = mynx*myny
allocate(myidx(nlr))
! Now, let's generate the list of indices I own
nr = 0
do i=bndx(iamx),bndx(iamx+1)-1
do j=bndy(iamy),bndx(iamy+1)-1
nr = nr + 1
call ijk2idx(myidx(nr),i,j,idim,idim)
end do
end do
if (nr /= nlr) then
write(psb_err_unit,*) iam,iamx,iamy, 'Initialization error: NR vs NLR ',&
& nr,nlr,mynx,myny
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
end if
!
! Third example of use of CDALL: specify for each process
! the set of global indices it owns.
!
call psb_cdall(ictxt,desc_a,info,vl=myidx)
case default
write(psb_err_unit,*) iam, 'Initialization error: should not get here'
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
return
end if
call psb_barrier(ictxt)
t0 = psb_wtime()
call psb_cdall(ictxt,desc_a,info,nl=nr)
return
end select
if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz)
! define rhs from boundary conditions; also build initial guess
if (info == psb_success_) call psb_geall(xv,desc_a,info)
if (info == psb_success_) call psb_geall(bv,desc_a,info)
call psb_barrier(ictxt)
talc = psb_wtime()-t0
@@ -218,9 +325,6 @@ contains
endif
myidx = desc_a%get_global_indices()
nlr = size(myidx)
! loop over rows belonging to current process in a block
! distribution.
@@ -234,13 +338,8 @@ contains
! local matrix pointer
glob_row=myidx(i)
! compute gridpoint coordinates
if (mod(glob_row,(idim)) == 0) then
ix = glob_row/(idim)
else
ix = glob_row/(idim)+1
endif
iy = (glob_row-(ix-1)*idim)
! x, y
call idx2ijk(ix,iy,glob_row,idim,idim)
! x, y coordinates
x = (ix-1)*deltah
y = (iy-1)*deltah
test/pargen/psb_s_pde3d.f90
+143 -45
View File
@@ -50,10 +50,12 @@
!
! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation.
!
! In this sample program the index space of the discretized
! computational domain is first numbered sequentially in a standard way,
! then the corresponding vector is distributed according to a BLOCK
! data distribution.
! There are three choices available for data distribution:
! 1. A simple BLOCK distribution
! 2. A ditribution based on arbitrary assignment of indices to processes,
! typically from a graph partitioner
! 3. A 3D distribution in which the unit cube is partitioned
! into subcubes, each one assigned to a process.
!
!
module psb_s_pde3d_mod
@@ -85,14 +87,15 @@ contains
end function s_null_func_3d
!
!
! subroutine to allocate and fill in the coefficient matrix and
! the rhs.
!
!
! subroutine to allocate and fill in the coefficient matrix and
! the rhs.
!
subroutine psb_s_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,&
& a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,imold,nrl)
& a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,imold,partition,nrl,iv)
use psb_base_mod
use psb_util_mod
!
! Discretizes the partial differential equation
!
@@ -120,7 +123,7 @@ contains
class(psb_s_base_sparse_mat), optional :: amold
class(psb_s_base_vect_type), optional :: vmold
class(psb_i_base_vect_type), optional :: imold
integer(psb_ipk_), optional :: nrl
integer(psb_ipk_), optional :: partition, nrl,iv(:)
! Local variables.
@@ -129,9 +132,13 @@ contains
type(psb_s_coo_sparse_mat) :: acoo
type(psb_s_csr_sparse_mat) :: acsr
real(psb_spk_) :: zt(nb),x,y,z
integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k
integer(psb_ipk_) :: m,n,nnz,nr,nt,glob_row,nlr,i,j,ii,ib,k, partition_
integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
integer(psb_ipk_) :: np, iam, nr, nt
! For 3D partition
integer(psb_ipk_) :: npx,npy,npz, npdims(3),iamx,iamy,iamz,mynx,myny,mynz
integer(psb_ipk_), allocatable :: bndx(:),bndy(:),bndz(:)
! Process grid
integer(psb_ipk_) :: np, iam
integer(psb_ipk_) :: icoeff
integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:)
real(psb_spk_), allocatable :: val(:)
@@ -161,6 +168,17 @@ contains
sqdeltah = deltah*deltah
deltah2 = 2.e0* deltah
if (present(partition)) then
if ((1<= partition).and.(partition <= 3)) then
partition_ = partition
else
write(*,*) 'Invalid partition choice ',partition,' defaulting to 3'
partition_ = 3
end if
else
partition_ = 3
end if
! initialize array descriptor and sparse matrix storage. provide an
! estimate of the number of non zeroes
@@ -169,29 +187,121 @@ contains
nnz = ((n*9)/(np))
if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n
if (present(nrl)) then
nr = nrl
else
!
! Using a simple BLOCK distribution.
!
nt = (m+np-1)/np
nr = max(0,min(nt,m-(iam*nt)))
end if
select case(partition_)
case(1)
! A BLOCK partition
if (present(nrl)) then
nr = nrl
else
!
! Using a simple BLOCK distribution.
!
nt = (m+np-1)/np
nr = max(0,min(nt,m-(iam*nt)))
end if
nt = nr
call psb_sum(ictxt,nt)
if (nt /= m) then
write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m
nt = nr
call psb_sum(ictxt,nt)
if (nt /= m) then
write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
return
end if
!
! First example of use of CDALL: specify for each process a number of
! contiguous rows
!
call psb_cdall(ictxt,desc_a,info,nl=nr)
myidx = desc_a%get_global_indices()
nlr = size(myidx)
case(2)
! A partition defined by the user through IV
if (present(iv)) then
if (size(iv) /= m) then
write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
return
end if
else
write(psb_err_unit,*) iam, 'Initialization error: IV not present'
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
return
end if
!
! Second example of use of CDALL: specify for each row the
! process that owns it
!
call psb_cdall(ictxt,desc_a,info,vg=iv)
myidx = desc_a%get_global_indices()
nlr = size(myidx)
case(3)
! A 3-dimensional partition
! A nifty MPI function will split the process list
npdims = 0
call mpi_dims_create(np,3,npdims,info)
npx = npdims(1)
npy = npdims(2)
npz = npdims(3)
allocate(bndx(0:npx),bndy(0:npy),bndz(0:npz))
! We can reuse idx2ijk for process indices as well.
call idx2ijk(iamx,iamy,iamz,iam,npx,npy,npz,base=0)
! Now let's split the 3D cube in hexahedra
call dist1Didx(bndx,idim,npx)
mynx = bndx(iamx+1)-bndx(iamx)
call dist1Didx(bndy,idim,npy)
myny = bndy(iamy+1)-bndy(iamy)
call dist1Didx(bndz,idim,npz)
mynz = bndz(iamz+1)-bndz(iamz)
! How many indices do I own?
nlr = mynx*myny*mynz
allocate(myidx(nlr))
! Now, let's generate the list of indices I own
nr = 0
do i=bndx(iamx),bndx(iamx+1)-1
do j=bndy(iamy),bndx(iamy+1)-1
do k=bndz(iamz),bndz(iamz+1)-1
nr = nr + 1
call ijk2idx(myidx(nr),i,j,k,idim,idim,idim)
end do
end do
end do
if (nr /= nlr) then
write(psb_err_unit,*) iam,iamx,iamy,iamz, 'Initialization error: NR vs NLR ',&
& nr,nlr,mynx,myny,mynz
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
end if
!
! Third example of use of CDALL: specify for each process
! the set of global indices it owns.
!
call psb_cdall(ictxt,desc_a,info,vl=myidx)
case default
write(psb_err_unit,*) iam, 'Initialization error: should not get here'
info = -1
call psb_barrier(ictxt)
call psb_abort(ictxt)
return
end if
return
end select
call psb_barrier(ictxt)
t0 = psb_wtime()
call psb_cdall(ictxt,desc_a,info,nl=nr)
if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz)
! define rhs from boundary conditions; also build initial guess
if (info == psb_success_) call psb_geall(xv,desc_a,info)
@@ -219,8 +329,6 @@ contains
goto 9999
endif
myidx = desc_a%get_global_indices()
nlr = size(myidx)
! loop over rows belonging to current process in a block
! distribution.
@@ -235,18 +343,8 @@ contains
! local matrix pointer
glob_row=myidx(i)
! compute gridpoint coordinates
if (mod(glob_row,(idim*idim)) == 0) then
ix = glob_row/(idim*idim)
else
ix = glob_row/(idim*idim)+1
endif
if (mod((glob_row-(ix-1)*idim*idim),idim) == 0) then
iy = (glob_row-(ix-1)*idim*idim)/idim
else
iy = (glob_row-(ix-1)*idim*idim)/idim+1
endif
iz = glob_row-(ix-1)*idim*idim-(iy-1)*idim
! x, y, x coordinates
call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
! x, y, z coordinates
x = (ix-1)*deltah
y = (iy-1)*deltah
z = (iz-1)*deltah
test/pargen/tryidxijk.f90
+1 -106
View File
@@ -1,5 +1,5 @@
program tryyidxijk
use psb_base_mod, only : psb_ipk_
use psb_util_mod
implicit none
integer(psb_ipk_) :: nx,ny,nz, base, i, j, k, idx,npx,npy,npz
@@ -41,109 +41,4 @@ program tryyidxijk
write(*,*) 'SZZ:',v(1:npz)-v(0:npz-1)
contains
!
! 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 idx2ijk(i,j,k,idx,nx,ny,nz,base)
use psb_base_mod, only : psb_ipk_
implicit none
integer(psb_ipk_), intent(out) :: i,j,k
integer(psb_ipk_), intent(in) :: idx,nx,ny,nz
integer(psb_ipk_), intent(in), optional :: base
integer(psb_ipk_) :: base_, idx_
if (present(base)) then
base_ = base
else
base_ = 1
end if
idx_ = idx - base_
k = mod(idx_,nz) + base_
j = mod(idx_/nz,ny) + base_
i = mod(idx_/(ny*nz),nx) + base_
end subroutine idx2ijk
!
! Given a triple (I,J,K) and the domain size (NX,NY,NZ)
! compute the global index IDX
! 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 ijk2idx(idx,i,j,k,nx,ny,nz,base)
use psb_base_mod, only : psb_ipk_
implicit none
integer(psb_ipk_), intent(out) :: idx
integer(psb_ipk_), intent(in) :: i,j,k,nx,ny,nz
integer(psb_ipk_), intent(in), optional :: base
integer(psb_ipk_) :: base_
if (present(base)) then
base_ = base
else
base_ = 1
end if
idx = ((i-base_)*nz*ny + (j-base_)*nz + k - base_) + base_
end subroutine ijk2idx
!
! dist1Didx
! Given an index space [base : N-(1-base)] and
! a set of NP processes, split the index base as
! evenly as possible, i.e. difference in size
! between any two processes is either 0 or 1,
! then return the boundaries in a vector
! such that
! V(P) : first index owned by process P
! V(P+1) : first index owned by process P+1
!
subroutine dist1Didx(v,n,np,base)
use psb_base_mod, only : psb_ipk_
implicit none
integer(psb_ipk_), intent(out) :: v(:)
integer(psb_ipk_), intent(in) :: n, np
integer(psb_ipk_), intent(in), optional :: base
!
integer(psb_ipk_) :: base_, nb, i
if (present(base)) then
base_ = base
else
base_ = 1
end if
nb = n/np
do i=1,mod(n,np)
v(i) = nb + 1
end do
do i=mod(n,np)+1,np
v(i) = nb
end do
v(2:np+1) = v(1:np)
v(1) = base_
do i=2,np+1
v(i) = v(i) + v(i-1)
end do
end subroutine dist1Didx
end program tryyidxijk
util/Makefile
+1 -1
View File
@@ -7,7 +7,7 @@ MODDIR=../modules
HERE=.
BASEOBJS= psb_blockpart_mod.o psb_metispart_mod.o \
BASEOBJS= psb_blockpart_mod.o psb_metispart_mod.o psb_partidx_mod.o \
psb_hbio_mod.o psb_mmio_mod.o psb_mat_dist_mod.o \
psb_s_mat_dist_mod.o psb_d_mat_dist_mod.o psb_c_mat_dist_mod.o psb_z_mat_dist_mod.o \
psb_renum_mod.o psb_gps_mod.o
util/psb_partidx_mod.F90
+259
View File
@@ -0,0 +1,259 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
!
! Purpose:
! Proide functions to handle a distribution of a general
! rectangular 2/3/n-dimensional domain onto
! a rectangular 2/3/n-dimensional grid of processes
!
! See test/pargen/psb_X_pdeNd for examples of usage
!
module psb_partidx_mod
use psb_base_mod, only : psb_ipk_
interface idx2ijk
module procedure idx2ijk3d, idx2ijkv, idx2ijk2d
end interface idx2ijk
interface ijk2idx
module procedure ijk2idx3d, ijk2idxv, ijk2idx2d
end interface ijk2idx
contains
!
! 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 idx2ijk3d(i,j,k,idx,nx,ny,nz,base)
use psb_base_mod, only : psb_ipk_
implicit none
integer(psb_ipk_), intent(out) :: i,j,k
integer(psb_ipk_), intent(in) :: idx,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 idx2ijk3d
subroutine idx2ijk2d(i,j,idx,nx,ny,base)
use psb_base_mod, only : psb_ipk_
implicit none
integer(psb_ipk_), intent(out) :: i,j
integer(psb_ipk_), intent(in) :: idx,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 idx2ijk2d
!
! 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 idx2ijkv(coords,idx,dims,base)
use psb_base_mod, only : psb_ipk_
implicit none
integer(psb_ipk_), intent(out) :: coords(:)
integer(psb_ipk_), intent(in) :: idx,dims(:)
integer(psb_ipk_), intent(in), optional :: base
integer(psb_ipk_) :: base_, idx_, 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 idx2ijkv
!
! Given a triple (I,J,K) and the domain size (NX,NY,NZ)
! compute the global index IDX
! 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 ijk2idxv(idx,coords,dims,base)
use psb_base_mod, only : psb_ipk_
implicit none
integer(psb_ipk_), intent(in) :: coords(:),dims(:)
integer(psb_ipk_), intent(out) :: idx
integer(psb_ipk_), intent(in), optional :: base
integer(psb_ipk_) :: base_, i, sz
if (present(base)) then
base_ = base
else
base_ = 1
end if
sz = size(coords)
if (sz /= size(dims)) then
write(0,*) 'Error: size mismatch ',size(coords),size(dims)
idx = 0
return
end if
idx = coords(1) - base_
do i=2, sz
idx = (idx * dims(i)) + coords(i) - base_
end do
idx = idx + base_
end subroutine ijk2idxv
!
! Given a triple (I,J,K) and the domain size (NX,NY,NZ)
! compute the global index IDX
! 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 ijk2idx3d(idx,i,j,k,nx,ny,nz,base)
use psb_base_mod, only : psb_ipk_
implicit none
integer(psb_ipk_), intent(out) :: idx
integer(psb_ipk_), intent(in) :: i,j,k,nx,ny,nz
integer(psb_ipk_), intent(in), optional :: base
! idx = ((i-base_)*nz*ny + (j-base_)*nz + k - base_) + base_
call ijk2idx(idx,[i,j,k],[nx,ny,nz],base)
end subroutine ijk2idx3d
subroutine ijk2idx2d(idx,i,j,nx,ny,base)
use psb_base_mod, only : psb_ipk_
implicit none
integer(psb_ipk_), intent(out) :: idx
integer(psb_ipk_), intent(in) :: i,j,nx,ny
integer(psb_ipk_), intent(in), optional :: base
! idx = ((i-base_)*ny + (j-base_) + base_
call ijk2idx(idx,[i,j],[nx,ny],base)
end subroutine ijk2idx2d
!
! dist1Didx
! Given an index space [base : N-(1-base)] and
! a set of NP processes, split the index base as
! evenly as possible, i.e. difference in size
! between any two processes is either 0 or 1,
! then return the boundaries in a vector
! such that
! V(P) : first index owned by process P
! V(P+1) : first index owned by process P+1
!
subroutine dist1Didx(v,n,np,base)
use psb_base_mod, only : psb_ipk_
implicit none
integer(psb_ipk_), intent(out) :: v(:)
integer(psb_ipk_), intent(in) :: n, np
integer(psb_ipk_), intent(in), optional :: base
!
integer(psb_ipk_) :: base_, nb, i
if (present(base)) then
base_ = base
else
base_ = 1
end if
nb = n/np
do i=1,mod(n,np)
v(i) = nb + 1
end do
do i=mod(n,np)+1,np
v(i) = nb
end do
v(2:np+1) = v(1:np)
v(1) = base_
do i=2,np+1
v(i) = v(i) + v(i-1)
end do
end subroutine dist1Didx
end module psb_partidx_mod
util/psb_util_mod.f90
+1 -2
View File
@@ -34,11 +34,10 @@
module psb_util_mod
use psb_blockpart_mod
use psb_metispart_mod
use psb_partidx_mod
use psb_hbio_mod
use psb_mmio_mod
use psb_mat_dist_mod
use psb_renum_mod
!!$ use psb_d_genpde_mod
!!$ use psb_s_genpde_mod
end module psb_util_mod