diff --git a/base/comm/Makefile b/base/comm/Makefile index 5b6b5786..3a587281 100644 --- a/base/comm/Makefile +++ b/base/comm/Makefile @@ -1,9 +1,9 @@ include ../../Make.inc OBJS = psb_dgather.o psb_dhalo.o psb_dovrl.o \ - psb_ihalo.o psb_zgather.o psb_zhalo.o psb_zovrl.o + psb_igather.o psb_ihalo.o psb_zgather.o psb_zhalo.o psb_zovrl.o -MPFOBJS = psb_dscatter.o psb_zscatter.o +MPFOBJS = psb_dscatter.o psb_zscatter.o psb_iscatter.o LIBDIR = .. MODDIR = ../modules INCDIRS = -I $(LIBDIR) -I $(MODDIR) -I . diff --git a/base/comm/psb_dgather.f90 b/base/comm/psb_dgather.f90 index 863172eb..ffac9e34 100644 --- a/base/comm/psb_dgather.f90 +++ b/base/comm/psb_dgather.f90 @@ -39,20 +39,13 @@ ! locx - real,dimension(:,:). The local piece of the distributed ! matrix to be gathered. ! desc_a - type(). The communication descriptor. -! info - integer. Eventually returns an error code. +! info - integer. Error code. ! iroot - integer. The process that has to own the -! global matrix. If -1 all +! global matrix. If -1 all ! the processes will have a copy. -! iiglobx - integer(optional). The starting row of the global matrix. -! ijglobx - integer(optional). The starting column of the global matrix. -! iilocx - integer(optional). The starting row of the local piece -! of matrix. -! ijlocx - integer(optional). The starting column of the local piece -! of matrix. -! ik - integer(optional). The number of columns to gather. +! Default: -1. ! -subroutine psb_dgatherm(globx, locx, desc_a, info, iroot,& - & iiglobx, ijglobx, iilocx,ijlocx,ik) +subroutine psb_dgatherm(globx, locx, desc_a, info, iroot) use psb_descriptor_type use psb_check_mod use psb_error_mod @@ -63,7 +56,7 @@ subroutine psb_dgatherm(globx, locx, desc_a, info, iroot,& real(kind(1.d0)), intent(out) :: globx(:,:) type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info - integer, intent(in), optional :: iroot, iiglobx, ijglobx, iilocx, ijlocx, ik + integer, intent(in), optional :: iroot ! locals @@ -104,29 +97,10 @@ subroutine psb_dgatherm(globx, locx, desc_a, info, iroot,& iiroot = root endif - if (present(iiglobx)) then - iglobx = iiglobx - else - iglobx = 1 - end if - - if (present(ijglobx)) then - jglobx = ijglobx - else - jglobx = 1 - end if - - if (present(iilocx)) then - ilocx = iilocx - else - ilocx = 1 - end if - - if (present(ijlocx)) then - jlocx = ijlocx - else - jlocx = 1 - end if + iglobx = 1 + jglobx = 1 + ilocx = 1 + jlocx = 1 lda_globx = size(globx,1) lda_locx = size(locx, 1) @@ -137,16 +111,7 @@ subroutine psb_dgatherm(globx, locx, desc_a, info, iroot,& lock=size(locx,2)-jlocx+1 globk=size(globx,2)-jglobx+1 maxk=min(lock,globk) - - if(present(ik)) then - if(ik.gt.maxk) then - k=maxk - else - k=ik - end if - else - k = maxk - end if + k = maxk call psb_bcast(ictxt,k,root=iiroot) @@ -245,16 +210,12 @@ end subroutine psb_dgatherm ! locx - real,dimension(:). The local piece of the ditributed ! vector to be gathered. ! desc_a - type(). The communication descriptor. -! info - integer. Eventually returns an error code. +! info - integer. Error code. ! iroot - integer. The process that has to own the -! global vector. If -1 all +! global matrix. If -1 all ! the processes will have a copy. -! iiglobx - integer(optional). The starting row of the global vector. -! iilocx - integer(optional). The starting row of the local piece -! of vector. ! -subroutine psb_dgatherv(globx, locx, desc_a, info, iroot,& - & iiglobx, iilocx) +subroutine psb_dgatherv(globx, locx, desc_a, info, iroot) use psb_descriptor_type use psb_check_mod use psb_error_mod @@ -265,7 +226,7 @@ subroutine psb_dgatherv(globx, locx, desc_a, info, iroot,& real(kind(1.d0)), intent(out) :: globx(:) type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info - integer, intent(in), optional :: iroot, iiglobx, iilocx + integer, intent(in), optional :: iroot ! locals @@ -303,18 +264,9 @@ subroutine psb_dgatherv(globx, locx, desc_a, info, iroot,& end if jglobx=1 - if (present(iiglobx)) then - iglobx = iiglobx - else - iglobx = 1 - end if - + iglobx = 1 jlocx=1 - if (present(iilocx)) then - ilocx = iilocx - else - ilocx = 1 - end if + ilocx = 1 lda_globx = size(globx) lda_locx = size(locx) diff --git a/base/comm/psb_dscatter.f90 b/base/comm/psb_dscatter.f90 index b344f88b..2076dd9f 100644 --- a/base/comm/psb_dscatter.f90 +++ b/base/comm/psb_dscatter.f90 @@ -38,17 +38,11 @@ ! globx - real,dimension(:,:). The global matrix to scatter. ! locx - real,dimension(:,:). The local piece of the ditributed matrix. ! desc_a - type(). The communication descriptor. -! info - integer. Eventually returns an error code. +! info - integer. Error code. ! iroot - integer(optional). The process that owns the global matrix. If -1 all -! the processes have a copy. -! iiglobx - integer(optional). The starting row of the global matrix. -! ijglobx - integer(optional). The starting column of the global matrix. -! iilocx - integer(optional). The starting row of the local piece of matrix. -! ijlocx - integer(optional). The starting column of the local piece of matrix. -! ik - integer(optional). The number of columns to gather. +! the processes have a copy. Default -1. ! -subroutine psb_dscatterm(globx, locx, desc_a, info, iroot,& - & iiglobx, ijglobx, iilocx,ijlocx,ik) +subroutine psb_dscatterm(globx, locx, desc_a, info, iroot) use psb_descriptor_type use psb_check_mod @@ -61,8 +55,7 @@ subroutine psb_dscatterm(globx, locx, desc_a, info, iroot,& real(kind(1.d0)), intent(in) :: globx(:,:) type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info - integer, intent(in), optional :: iroot,iiglobx,& - & ijglobx,iilocx,ijlocx,ik + integer, intent(in), optional :: iroot ! locals @@ -70,8 +63,8 @@ subroutine psb_dscatterm(globx, locx, desc_a, info, iroot,& & err_act, m, n, i, j, idx, nrow, iiroot, iglobx, jglobx,& & ilocx, jlocx, lda_locx, lda_globx, lock, globk, icomm, k, maxk, root, ilx,& & jlx, myrank, rootrank, c, pos - real(kind(1.d0)),pointer :: scatterv(:) - integer, pointer :: displ(:), l_t_g_all(:), all_dim(:) + real(kind(1.d0)), allocatable :: scatterv(:) + integer, allocatable :: displ(:), l_t_g_all(:), all_dim(:) character(len=20) :: name, ch_err name='psb_scatterm' @@ -104,30 +97,10 @@ subroutine psb_dscatterm(globx, locx, desc_a, info, iroot,& iiroot=0 endif - if (present(iiglobx)) then - iglobx = iiglobx - else - iglobx = 1 - end if - - if (present(ijglobx)) then - jglobx = ijglobx - else - jglobx = 1 - end if - - if (present(iilocx)) then - ilocx = iilocx - else - ilocx = 1 - end if - - if (present(ijlocx)) then - jlocx = ijlocx - else - jlocx = 1 - end if - + iglobx = 1 + jglobx = 1 + ilocx = 1 + jlocx = 1 lda_globx = size(globx,1) lda_locx = size(locx, 1) @@ -137,17 +110,7 @@ subroutine psb_dscatterm(globx, locx, desc_a, info, iroot,& lock=size(locx,2)-jlocx+1 globk=size(globx,2)-jglobx+1 maxk=min(lock,globk) - - if(present(ik)) then - if(ik.gt.maxk) then - k=maxk - else - k=ik - end if - else - k = maxk - end if - + k = maxk call psb_get_mpicomm(ictxt,icomm) call psb_get_rank(myrank,ictxt,me) @@ -303,7 +266,7 @@ end subroutine psb_dscatterm ! globx - real,dimension(:). The global vector to scatter. ! locx - real,dimension(:). The local piece of the ditributed vector. ! desc_a - type(). The communication descriptor. -! info - integer. Eventually returns an error code. +! info - integer. Error code. ! iroot - integer(optional). The process that owns the global vector. If -1 all ! the processes have a copy. ! @@ -327,8 +290,8 @@ subroutine psb_dscatterv(globx, locx, desc_a, info, iroot) & err_act, m, n, i, j, idx, nrow, iiroot, iglobx, jglobx,& & ilocx, jlocx, lda_locx, lda_globx, root, k, icomm, myrank,& & rootrank, pos, ilx, jlx - real(kind(1.d0)),pointer :: scatterv(:) - integer, pointer :: displ(:), l_t_g_all(:), all_dim(:) + real(kind(1.d0)), allocatable :: scatterv(:) + integer, allocatable :: displ(:), l_t_g_all(:), all_dim(:) character(len=20) :: name, ch_err name='psb_scatterv' diff --git a/base/comm/psb_igather.f90 b/base/comm/psb_igather.f90 new file mode 100644 index 00000000..b975e4f1 --- /dev/null +++ b/base/comm/psb_igather.f90 @@ -0,0 +1,326 @@ +!!$ +!!$ Parallel Sparse BLAS v2.0 +!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ +!!$ 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. +!!$ +!!$ +! File: psb_igather.f90 +! +! Subroutine: psb_igatherm +! This subroutine gathers pieces of a distributed dense matrix into a local one. +! +! Parameters: +! globx - integer,dimension(:,:). The local matrix into which gather +! the distributed pieces. +! locx - integer,dimension(:,:). The local piece of the distributed +! matrix to be gathered. +! desc_a - type(). The communication descriptor. +! info - integer. Error code. +! iroot - integer. The process that has to own the +! global matrix. If -1 all +! the processes will have a copy. +! Default: -1. +! +subroutine psb_igatherm(globx, locx, desc_a, info, iroot) + use psb_descriptor_type + use psb_check_mod + use psb_error_mod + use psb_penv_mod + implicit none + + integer, intent(in) :: locx(:,:) + integer, intent(out) :: globx(:,:) + type(psb_desc_type), intent(in) :: desc_a + integer, intent(out) :: info + integer, intent(in), optional :: iroot + + + ! locals + integer :: int_err(5), ictxt, np, me,& + & err_act, n, root, iiroot, ilocx, iglobx, jlocx,& + & jglobx, lda_locx, lda_globx, m, lock, globk, maxk, k, jlx, ilx, i, j, idx + character(len=20) :: name, ch_err + + name='psb_igatherm' + if(psb_get_errstatus().ne.0) return + info=0 + call psb_erractionsave(err_act) + + ictxt=psb_cd_get_context(desc_a) + + ! check on blacs grid + call psb_info(ictxt, me, np) + if (np == -1) then + info = 2010 + call psb_errpush(info,name) + goto 9999 + endif + + if (present(iroot)) then + root = iroot + if((root.lt.-1).or.(root.gt.np)) then + info=30 + int_err(1:2)=(/5,root/) + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + else + root = -1 + end if + if (root==-1) then + iiroot=0 + else + iiroot = root + endif + + iglobx = 1 + jglobx = 1 + ilocx = 1 + jlocx = 1 + + lda_globx = size(globx,1) + lda_locx = size(locx, 1) + + m = psb_cd_get_global_rows(desc_a) + n = psb_cd_get_global_cols(desc_a) + + lock=size(locx,2)-jlocx+1 + globk=size(globx,2)-jglobx+1 + maxk=min(lock,globk) + k = maxk + + call psb_bcast(ictxt,k,root=iiroot) + + ! there should be a global check on k here!!! + + call psb_chkglobvect(m,n,size(globx,1),iglobx,jglobx,desc_a,info) + call psb_chkvect(m,n,size(locx,1),ilocx,jlocx,desc_a,info,ilx,jlx) + if(info.ne.0) then + info=4010 + ch_err='psb_chk(glob)vect' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if ((ilx.ne.1).or.(iglobx.ne.1)) then + info=3040 + call psb_errpush(info,name) + goto 9999 + end if + + globx(:,:)=0 + + do j=1,k + do i=1,psb_cd_get_local_rows(desc_a) + idx = desc_a%loc_to_glob(i) + globx(idx,jglobx+j-1) = locx(i,jlx+j-1) + end do + ! adjust overlapped elements + i=1 + do while (desc_a%ovrlap_elem(i).ne.-1) + idx=desc_a%ovrlap_elem(i+psb_ovrlp_elem_) + idx=desc_a%loc_to_glob(idx) + globx(idx,jglobx+j-1) = & + & globx(idx,jglobx+j-1)/desc_a%ovrlap_elem(i+psb_n_dom_ovr_) + i=i+2 + end do + end do + + call psb_sum(ictxt,globx(1:m,jglobx:jglobx+k-1),root=root) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act.eq.psb_act_abort_) then + call psb_error(ictxt) + return + end if + return + +end subroutine psb_igatherm + + + + + + +!!$ +!!$ Parallel Sparse BLAS v2.0 +!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ +!!$ 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. +!!$ +!!$ +! Subroutine: psb_igatherv +! This subroutine gathers pieces of a distributed dense vector into a local one. +! +! Parameters: +! globx - integer,dimension(:). The local vector into which gather the +! distributed pieces. +! locx - integer,dimension(:). The local piece of the ditributed +! vector to be gathered. +! desc_a - type(). The communication descriptor. +! info - integer. Error code. +! iroot - integer. The process that has to own the +! global matrix. If -1 all +! the processes will have a copy. +! +subroutine psb_igatherv(globx, locx, desc_a, info, iroot) + use psb_descriptor_type + use psb_check_mod + use psb_error_mod + use psb_penv_mod + implicit none + + integer, intent(in) :: locx(:) + integer, intent(out) :: globx(:) + type(psb_desc_type), intent(in) :: desc_a + integer, intent(out) :: info + integer, intent(in), optional :: iroot + + + ! locals + integer :: int_err(5), ictxt, np, me, & + & err_act, n, root, ilocx, iglobx, jlocx,& + & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx + + character(len=20) :: name, ch_err + + name='psb_igatherv' + if(psb_get_errstatus().ne.0) return + info=0 + call psb_erractionsave(err_act) + + ictxt=psb_cd_get_context(desc_a) + + ! check on blacs grid + call psb_info(ictxt, me, np) + if (np == -1) then + info = 2010 + call psb_errpush(info,name) + goto 9999 + endif + + if (present(iroot)) then + root = iroot + if((root.lt.-1).or.(root.gt.np)) then + info=30 + int_err(1:2)=(/5,root/) + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + else + root = -1 + end if + + jglobx=1 + iglobx = 1 + jlocx=1 + ilocx = 1 + + lda_globx = size(globx) + lda_locx = size(locx) + + m = psb_cd_get_global_rows(desc_a) + n = psb_cd_get_global_cols(desc_a) + + k = 1 + + + ! there should be a global check on k here!!! + + call psb_chkglobvect(m,n,size(globx),iglobx,jglobx,desc_a,info) + call psb_chkvect(m,n,size(locx),ilocx,jlocx,desc_a,info,ilx,jlx) + if(info.ne.0) then + info=4010 + ch_err='psb_chk(glob)vect' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if ((ilx.ne.1).or.(iglobx.ne.1)) then + info=3040 + call psb_errpush(info,name) + goto 9999 + end if + + globx(:)=0 + + do i=1,psb_cd_get_local_rows(desc_a) + idx = desc_a%loc_to_glob(i) + globx(idx) = locx(i) + end do + ! adjust overlapped elements + i=1 + do while (desc_a%ovrlap_elem(i).ne.-1) + idx=desc_a%ovrlap_elem(i+psb_ovrlp_elem_) + idx=desc_a%loc_to_glob(idx) + globx(idx) = globx(idx)/desc_a%ovrlap_elem(i+psb_n_dom_ovr_) + i=i+2 + end do + + call psb_sum(ictxt,globx(1:m),root=root) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act.eq.psb_act_abort_) then + call psb_error(ictxt) + return + end if + return + +end subroutine psb_igatherv diff --git a/base/comm/psb_iscatter.f90 b/base/comm/psb_iscatter.f90 new file mode 100644 index 00000000..e80fcfc8 --- /dev/null +++ b/base/comm/psb_iscatter.f90 @@ -0,0 +1,416 @@ +!!$ +!!$ Parallel Sparse BLAS v2.0 +!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ +!!$ 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. +!!$ +!!$ +! File: psb_iscatter.f90 +! +! Subroutine: psb_iscatterm +! This subroutine scatters a global matrix locally owned by one process +! into pieces that are local to alle the processes. +! +! Parameters: +! globx - integer,dimension(:,:). The global matrix to scatter. +! locx - integer,dimension(:,:). The local piece of the ditributed matrix. +! desc_a - type(). The communication descriptor. +! info - integer. Error code. +! iroot - integer(optional). The process that owns the global matrix. If -1 all +! the processes have a copy. +! +subroutine psb_iscatterm(globx, locx, desc_a, info, iroot) + + use psb_descriptor_type + use psb_check_mod + use psb_error_mod + use mpi + use psb_penv_mod + implicit none + + integer, intent(out) :: locx(:,:) + integer, intent(in) :: globx(:,:) + type(psb_desc_type), intent(in) :: desc_a + integer, intent(out) :: info + integer, intent(in), optional :: iroot + + ! locals + integer :: int_err(5), ictxt, np, me,& + & err_act, m, n, i, j, idx, nrow, iiroot, iglobx, jglobx,& + & ilocx, jlocx, lda_locx, lda_globx, lock, globk, icomm, k, maxk, root, ilx,& + & jlx, myrank, rootrank, c, pos + integer, allocatable :: scatterv(:) + integer, allocatable :: displ(:), l_t_g_all(:), all_dim(:) + character(len=20) :: name, ch_err + + name='psb_scatterm' + if(psb_get_errstatus() /= 0) return + info=0 + call psb_erractionsave(err_act) + + ictxt=psb_cd_get_context(desc_a) + + ! check on blacs grid + call psb_info(ictxt, me, np) + if (np == -1) then + info = 2010 + call psb_errpush(info,name) + goto 9999 + endif + + if (present(iroot)) then + root = iroot + if((root.lt.-1).or.(root.gt.np)) then + info=30 + int_err(1:2)=(/5,root/) + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + else + root = -1 + end if + if (root==-1) then + iiroot=0 + endif + + iglobx = 1 + jglobx = 1 + ilocx = 1 + jlocx = 1 + lda_globx = size(globx,1) + lda_locx = size(locx, 1) + + m = psb_cd_get_global_rows(desc_a) + n = psb_cd_get_global_cols(desc_a) + + lock=size(locx,2)-jlocx+1 + globk=size(globx,2)-jglobx+1 + maxk=min(lock,globk) + k = maxk + call psb_get_mpicomm(ictxt,icomm) + call psb_get_rank(myrank,ictxt,me) + + + lda_globx = size(globx) + lda_locx = size(locx) + + m = psb_cd_get_global_rows(desc_a) + n = psb_cd_get_global_cols(desc_a) + + if (me == iiroot) then + call igebs2d(ictxt, 'all', ' ', 1, 1, k, 1) + else + call igebr2d(ictxt, 'all', ' ', 1, 1, k, 1, iiroot, 0) + end if + + ! there should be a global check on k here!!! + + call psb_chkglobvect(m,n,size(globx),iglobx,jglobx,desc_a,info) + if (info == 0) call psb_chkvect(m,n,size(locx),ilocx,jlocx,desc_a,info,ilx,jlx) + if(info /= 0) then + info=4010 + ch_err='psb_chk(glob)vect' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if ((ilx /= 1).or.(iglobx /= 1)) then + info=3040 + call psb_errpush(info,name) + goto 9999 + end if + + nrow=psb_cd_get_local_rows(desc_a) + + if(root == -1) then + ! extract my chunk + do j=1,k + do i=1, nrow + idx=desc_a%loc_to_glob(i) + locx(i,jlocx+j-1)=globx(idx,jglobx+j-1) + end do + end do + else + call psb_get_rank(rootrank,ictxt,root) + end if + + ! root has to gather size information + allocate(displ(np),all_dim(np),stat=info) + if(info /= 0) then + info=4010 + ch_err='Allocate' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call mpi_gather(nrow,1,mpi_integer,all_dim,& + & np,mpi_integer,rootrank,icomm,info) + + displ(1)=1 + displ(2:)=all_dim(1:np-1)+1 + + ! root has to gather loc_glob from each process + if(me == root) then + allocate(l_t_g_all(sum(all_dim)),scatterv(sum(all_dim)),stat=info) + if(info /= 0) then + info=4010 + ch_err='Allocate' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + end if + + call mpi_gatherv(desc_a%loc_to_glob,nrow,& + & mpi_integer,l_t_g_all,all_dim,& + & displ,mpi_integer,rootrank,icomm,info) + + + do c=1, k + ! prepare vector to scatter + if(me == root) then + do i=1,np + pos=displ(i) + do j=1, all_dim(i) + idx=l_t_g_all(pos+j-1) + scatterv(pos+j-1)=globx(idx,jglobx+c-1) + end do + end do + end if + + ! scatter !!! + call mpi_scatterv(scatterv,all_dim,displ,& + & mpi_integer,locx(1,jlocx+c-1),nrow,& + & mpi_integer,rootrank,icomm,info) + + end do + + deallocate(all_dim, l_t_g_all, displ, scatterv) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return + +end subroutine psb_iscatterm + + + + +!!$ +!!$ Parallel Sparse BLAS v2.0 +!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ +!!$ 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. +!!$ +!!$ + +! Subroutine: psb_iscatterv +! This subroutine scatters a global vector locally owned by one process +! into pieces that are local to alle the processes. +! +! Parameters: +! globx - integer,dimension(:). The global vector to scatter. +! locx - integer,dimension(:). The local piece of the ditributed vector. +! desc_a - type(). The communication descriptor. +! info - integer. Error code. +! iroot - integer(optional). The process that owns the global vector. If -1 all +! the processes have a copy. +! +subroutine psb_iscatterv(globx, locx, desc_a, info, iroot) + use psb_descriptor_type + use psb_check_mod + use psb_error_mod + use mpi + use psb_penv_mod + implicit none + + integer, intent(out) :: locx(:) + integer, intent(in) :: globx(:) + type(psb_desc_type), intent(in) :: desc_a + integer, intent(out) :: info + integer, intent(in), optional :: iroot + + + ! locals + integer :: int_err(5), ictxt, np, me, & + & err_act, m, n, i, j, idx, nrow, iiroot, iglobx, jglobx,& + & ilocx, jlocx, lda_locx, lda_globx, root, k, icomm, myrank,& + & rootrank, pos, ilx, jlx + integer, allocatable :: scatterv(:) + integer, allocatable :: displ(:), l_t_g_all(:), all_dim(:) + character(len=20) :: name, ch_err + + name='psb_scatterv' + if(psb_get_errstatus() /= 0) return + info=0 + call psb_erractionsave(err_act) + + ictxt=psb_cd_get_context(desc_a) + + ! check on blacs grid + call psb_info(ictxt, me, np) + if (np == -1) then + info = 2010 + call psb_errpush(info,name) + goto 9999 + endif + + if (present(iroot)) then + root = iroot + if((root.lt.-1).or.(root.gt.np)) then + info=30 + int_err(1:2)=(/5,root/) + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + else + root = -1 + end if + + call psb_get_mpicomm(ictxt,icomm) + call psb_get_rank(myrank,ictxt,me) + + lda_globx = size(globx) + lda_locx = size(locx) + + m = psb_cd_get_global_rows(desc_a) + n = psb_cd_get_global_cols(desc_a) + + k = 1 + + if (me == iiroot) then + call igebs2d(ictxt, 'all', ' ', 1, 1, k, 1) + else + call igebr2d(ictxt, 'all', ' ', 1, 1, k, 1, iiroot, 0) + end if + + ! there should be a global check on k here!!! + + call psb_chkglobvect(m,n,size(globx),iglobx,jglobx,desc_a,info) + if (info == 0) & + & call psb_chkvect(m,n,size(locx),ilocx,jlocx,desc_a,info,ilx,jlx) + if(info /= 0) then + info=4010 + ch_err='psb_chk(glob)vect' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if ((ilx /= 1).or.(iglobx /= 1)) then + info=3040 + call psb_errpush(info,name) + goto 9999 + end if + + nrow=psb_cd_get_local_rows(desc_a) + + if(root == -1) then + ! extract my chunk + do i=1, nrow + idx=desc_a%loc_to_glob(i) + locx(i)=globx(idx) + end do + else + call psb_get_rank(rootrank,ictxt,root) + end if + + ! root has to gather size information + allocate(displ(np),all_dim(np)) + call mpi_gather(nrow,1,mpi_integer,all_dim,& + & np,mpi_integer,rootrank,icomm,info) + + displ(1)=1 + displ(2:)=all_dim(1:np-1)+1 + + ! root has to gather loc_glob from each process + if(me == root) then + allocate(l_t_g_all(sum(all_dim)),scatterv(sum(all_dim))) + end if + + call mpi_gatherv(desc_a%loc_to_glob,nrow,& + & mpi_integer,l_t_g_all,all_dim,& + & displ,mpi_integer,rootrank,icomm,info) + + ! prepare vector to scatter + if(me == root) then + do i=1,np + pos=displ(i) + do j=1, all_dim(i) + idx=l_t_g_all(pos+j-1) + scatterv(pos+j-1)=globx(idx) + end do + end do + end if + + call mpi_scatterv(scatterv,all_dim,displ,& + & mpi_integer,locx,nrow,& + & mpi_integer,rootrank,icomm,info) + + deallocate(all_dim, l_t_g_all, displ, scatterv) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return + +end subroutine psb_iscatterv diff --git a/base/comm/psb_zgather.f90 b/base/comm/psb_zgather.f90 index 238496d4..a460a8bd 100644 --- a/base/comm/psb_zgather.f90 +++ b/base/comm/psb_zgather.f90 @@ -34,25 +34,17 @@ ! This subroutine gathers pieces of a distributed dense matrix into a local one. ! ! Parameters: -! globx - real,dimension(:,:). The local matrix into which gather +! globx - cplx,dimension(:,:). The local matrix into which gather ! the distributed pieces. -! locx - real,dimension(:,:). The local piece of the distributed +! locx - cplx,dimension(:,:). The local piece of the distributed ! matrix to be gathered. ! desc_a - type(). The communication descriptor. -! info - integer. Eventually returns an error code. +! info - integer. Error code. ! iroot - integer. The process that has to own the -! global matrix. If -1 all +! global matrix. If -1 all ! the processes will have a copy. -! iiglobx - integer(optional). The starting row of the global matrix. -! ijglobx - integer(optional). The starting column of the global matrix. -! iilocx - integer(optional). The starting row of the local piece -! of matrix. -! ijlocx - integer(optional). The starting column of the local piece -! of matrix. -! ik - integer(optional). The number of columns to gather. ! -subroutine psb_zgatherm(globx, locx, desc_a, info, iroot,& - & iiglobx, ijglobx, iilocx,ijlocx,ik) +subroutine psb_zgatherm(globx, locx, desc_a, info, iroot) use psb_descriptor_type use psb_check_mod use psb_error_mod @@ -63,7 +55,7 @@ subroutine psb_zgatherm(globx, locx, desc_a, info, iroot,& complex(kind(1.d0)), intent(out) :: globx(:,:) type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info - integer, intent(in), optional :: iroot, iiglobx, ijglobx, iilocx, ijlocx, ik + integer, intent(in), optional :: iroot ! locals @@ -105,29 +97,10 @@ subroutine psb_zgatherm(globx, locx, desc_a, info, iroot,& iiroot = root endif - if (present(iiglobx)) then - iglobx = iiglobx - else - iglobx = 1 - end if - - if (present(ijglobx)) then - jglobx = ijglobx - else - jglobx = 1 - end if - - if (present(iilocx)) then - ilocx = iilocx - else - ilocx = 1 - end if - - if (present(ijlocx)) then - jlocx = ijlocx - else - jlocx = 1 - end if + iglobx = 1 + jglobx = 1 + ilocx = 1 + jlocx = 1 lda_globx = size(globx,1) lda_locx = size(locx, 1) @@ -139,15 +112,7 @@ subroutine psb_zgatherm(globx, locx, desc_a, info, iroot,& globk=size(globx,2)-jglobx+1 maxk=min(lock,globk) - if(present(ik)) then - if(ik.gt.maxk) then - k=maxk - else - k=ik - end if - else - k = maxk - end if + k = maxk call psb_bcast(ictxt,k,root=iiroot) @@ -164,9 +129,9 @@ subroutine psb_zgatherm(globx, locx, desc_a, info, iroot,& end if if ((ilx.ne.1).or.(iglobx.ne.1)) then - info=3040 - call psb_errpush(info,name) - goto 9999 + info=3040 + call psb_errpush(info,name) + goto 9999 end if globx(:,:)=0.d0 @@ -242,21 +207,18 @@ end subroutine psb_zgatherm ! This subroutine gathers pieces of a distributed dense vector into a local one. ! ! Parameters: -! globx - real,dimension(:). The local vector into which gather +! globx - cplx,dimension(:). The local vector into which gather ! the distributed pieces. -! locx - real,dimension(:). The local piece of the distributed +! locx - cplx,dimension(:). The local piece of the distributed ! vector to be gathered. ! desc_a - type(). The communication descriptor. -! info - integer. Eventually returns an error code. +! info - integer. Error code. ! iroot - integer. The process that has to own the -! global vector. If -1 all +! global matrix. If -1 all ! the processes will have a copy. -! iiglobx - integer(optional). The starting row of the global vector. -! iilocx - integer(optional). The starting row of the local piece -! of vector. +! default: -1 ! -subroutine psb_zgatherv(globx, locx, desc_a, info, iroot,& - & iiglobx, iilocx) +subroutine psb_zgatherv(globx, locx, desc_a, info, iroot) use psb_descriptor_type use psb_check_mod use psb_error_mod @@ -267,7 +229,7 @@ subroutine psb_zgatherv(globx, locx, desc_a, info, iroot,& complex(kind(1.d0)), intent(out) :: globx(:) type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info - integer, intent(in), optional :: iroot, iiglobx, iilocx + integer, intent(in), optional :: iroot ! locals @@ -305,18 +267,9 @@ subroutine psb_zgatherv(globx, locx, desc_a, info, iroot,& end if jglobx=1 - if (present(iiglobx)) then - iglobx = iiglobx - else - iglobx = 1 - end if - + iglobx = 1 jlocx=1 - if (present(iilocx)) then - ilocx = iilocx - else - ilocx = 1 - end if + ilocx = 1 lda_globx = size(globx) lda_locx = size(locx) diff --git a/base/comm/psb_zscatter.f90 b/base/comm/psb_zscatter.f90 index 3c677834..0d46ca9a 100644 --- a/base/comm/psb_zscatter.f90 +++ b/base/comm/psb_zscatter.f90 @@ -38,17 +38,10 @@ ! globx - complex,dimension(:,:). The global matrix to scatter. ! locx - complex,dimension(:,:). The local piece of the ditributed matrix. ! desc_a - type(). The communication descriptor. -! info - integer. Eventually returns an error code. +! info - integer. Error code. ! iroot - integer(optional). The process that owns the global matrix. If -1 all -! the processes have a copy. -! iiglobx - integer(optional). The starting row of the global matrix. -! ijglobx - integer(optional). The starting column of the global matrix. -! iilocx - integer(optional). The starting row of the local piece of matrix. -! ijlocx - integer(optional). The starting column of the local piece of matrix. -! ik - integer(optional). The number of columns to gather. -! -subroutine psb_zscatterm(globx, locx, desc_a, info, iroot,& - & iiglobx, ijglobx, iilocx,ijlocx,ik) +! the processes have a copy. Default -1 +subroutine psb_zscatterm(globx, locx, desc_a, info, iroot) use psb_descriptor_type use psb_check_mod @@ -61,8 +54,7 @@ subroutine psb_zscatterm(globx, locx, desc_a, info, iroot,& complex(kind(1.d0)), intent(in) :: globx(:,:) type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info - integer, intent(in), optional :: iroot,iiglobx,& - & ijglobx,iilocx,ijlocx,ik + integer, intent(in), optional :: iroot ! locals @@ -70,8 +62,8 @@ subroutine psb_zscatterm(globx, locx, desc_a, info, iroot,& & err_act, m, n, i, j, idx, nrow, iiroot, iglobx, jglobx,& & ilocx, jlocx, lda_locx, lda_globx, lock, globk, icomm, k, maxk, root, ilx,& & jlx, myrank, rootrank, c, pos - complex(kind(1.d0)),pointer :: scatterv(:) - integer, pointer :: displ(:), l_t_g_all(:), all_dim(:) + complex(kind(1.d0)),allocatable :: scatterv(:) + integer, allocatable :: displ(:), l_t_g_all(:), all_dim(:) character(len=20) :: name, ch_err name='psb_scatterm' @@ -104,30 +96,12 @@ subroutine psb_zscatterm(globx, locx, desc_a, info, iroot,& iiroot=0 endif - if (present(iiglobx)) then - iglobx = iiglobx - else - iglobx = 1 - end if - - if (present(ijglobx)) then - jglobx = ijglobx - else - jglobx = 1 - end if - - if (present(iilocx)) then - ilocx = iilocx - else - ilocx = 1 - end if - - if (present(ijlocx)) then - jlocx = ijlocx - else - jlocx = 1 - end if + iglobx = 1 + jglobx = 1 + ilocx = 1 + jlocx = 1 + lda_globx = size(globx,1) lda_locx = size(locx, 1) @@ -137,16 +111,7 @@ subroutine psb_zscatterm(globx, locx, desc_a, info, iroot,& lock=size(locx,2)-jlocx+1 globk=size(globx,2)-jglobx+1 maxk=min(lock,globk) - - if(present(ik)) then - if(ik.gt.maxk) then - k=maxk - else - k=ik - end if - else - k = maxk - end if + k = maxk call psb_get_mpicomm(ictxt,icomm) call psb_get_rank(myrank,ictxt,me) @@ -327,8 +292,8 @@ subroutine psb_zscatterv(globx, locx, desc_a, info, iroot) & err_act, m, n, i, j, idx, nrow, iiroot, iglobx, jglobx,& & ilocx, jlocx, lda_locx, lda_globx, root, k, icomm, myrank,& & rootrank, c, pos, ilx, jlx - complex(kind(1.d0)),pointer :: scatterv(:) - integer, pointer :: displ(:), l_t_g_all(:), all_dim(:) + complex(kind(1.d0)), allocatable :: scatterv(:) + integer, allocatable :: displ(:), l_t_g_all(:), all_dim(:) character(len=20) :: name, ch_err name='psb_scatterv'