diff --git a/base/modules/psb_serial_mod.f90 b/base/modules/psb_serial_mod.f90 index 80f29212..8f2987b6 100644 --- a/base/modules/psb_serial_mod.f90 +++ b/base/modules/psb_serial_mod.f90 @@ -493,6 +493,40 @@ module psb_serial_mod integer, intent(in), optional :: lrw, nzin end subroutine psb_zspgetrow end interface + + + interface psb_gelp + ! 2-D version + subroutine psb_dgelp(trans,iperm,x,info) + real(kind(1.d0)), intent(inout) :: x(:,:) + integer, intent(in) :: iperm(:) + integer, intent(out) :: info + character, intent(in) :: trans + end subroutine psb_dgelp + ! 1-D version + subroutine psb_dgelpv(trans,iperm,x,info) + real(kind(1.d0)), intent(inout) :: x(:) + integer, intent(in) :: iperm(:) + integer, intent(out) :: info + character, intent(in) :: trans + end subroutine psb_dgelpv + ! 2-D version + subroutine psb_zgelp(trans,iperm,x,info) + complex(kind(1.d0)), intent(inout) :: x(:,:) + integer, intent(in) :: iperm(:) + integer, intent(out) :: info + character, intent(in) :: trans + end subroutine psb_zgelp + ! 1-D version + subroutine psb_zgelpv(trans,iperm,x,info) + complex(kind(1.d0)), intent(inout) :: x(:) + integer, intent(in) :: iperm(:) + integer, intent(out) :: info + character, intent(in) :: trans + end subroutine psb_zgelpv + end interface + + interface psb_msort module procedure imsort diff --git a/base/modules/psb_tools_mod.f90 b/base/modules/psb_tools_mod.f90 index 7f1c9803..904de70c 100644 --- a/base/modules/psb_tools_mod.f90 +++ b/base/modules/psb_tools_mod.f90 @@ -132,25 +132,27 @@ Module psb_tools_mod end interface interface psb_sphalo - Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) + Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& + & rowscale,colscale,outfmt,data) use psb_descriptor_type use psb_spmat_type Type(psb_dspmat_type),Intent(in) :: a Type(psb_dspmat_type),Intent(inout) :: blk Type(psb_desc_type),Intent(in),target :: desc_a integer, intent(out) :: info - logical, optional, intent(in) :: rwcnv,clcnv,cliprow + logical, optional, intent(in) :: rowcnv,colcnv,rowscale,colscale character(len=5), optional :: outfmt integer, intent(in), optional :: data end Subroutine psb_dsphalo - Subroutine psb_zsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) + Subroutine psb_zsphalo(a,desc_a,blk,info,rowcnv,colcnv,& + & rowscale,colscale,outfmt,data) use psb_descriptor_type use psb_spmat_type Type(psb_zspmat_type),Intent(in) :: a Type(psb_zspmat_type),Intent(inout) :: blk Type(psb_desc_type),Intent(in) :: desc_a integer, intent(out) :: info - logical, optional, intent(in) :: rwcnv,clcnv,cliprow + logical, optional, intent(in) :: rowcnv,colcnv,rowscale,colscale character(len=5), optional :: outfmt integer, intent(in), optional :: data end Subroutine psb_zsphalo @@ -227,42 +229,6 @@ Module psb_tools_mod end interface - interface psb_gelp - ! 2-D version - subroutine psb_dgelp(trans,iperm,x,desc_a,info) - use psb_descriptor_type - type(psb_desc_type), intent(in) :: desc_a - real(kind(1.d0)), intent(inout) :: x(:,:) - integer, intent(inout) :: iperm(:),info - character, intent(in) :: trans - end subroutine psb_dgelp - ! 1-D version - subroutine psb_dgelpv(trans,iperm,x,desc_a,info) - use psb_descriptor_type - type(psb_desc_type), intent(in) :: desc_a - real(kind(1.d0)), intent(inout) :: x(:) - integer, intent(inout) :: iperm(:), info - character, intent(in) :: trans - end subroutine psb_dgelpv - ! 2-D version - subroutine psb_zgelp(trans,iperm,x,desc_a,info) - use psb_descriptor_type - type(psb_desc_type), intent(in) :: desc_a - complex(kind(1.d0)), intent(inout) :: x(:,:) - integer, intent(inout) :: iperm(:),info - character, intent(in) :: trans - end subroutine psb_zgelp - ! 1-D version - subroutine psb_zgelpv(trans,iperm,x,desc_a,info) - use psb_descriptor_type - type(psb_desc_type), intent(in) :: desc_a - complex(kind(1.d0)), intent(inout) :: x(:) - integer, intent(inout) :: iperm(:), info - character, intent(in) :: trans - end subroutine psb_zgelpv - end interface - - interface psb_geins ! 2-D double precision version subroutine psb_dinsi(m,irw,val, x,desc_a,info,dupl) diff --git a/base/serial/Makefile b/base/serial/Makefile index 4c8f94af..86ff101f 100644 --- a/base/serial/Makefile +++ b/base/serial/Makefile @@ -12,7 +12,8 @@ FOBJS = psb_cest.o psb_dcoins.o psb_dcsmm.o psb_dcsmv.o \ psb_zfixcoo.o psb_zipcoo2csr.o psb_zipcsr2coo.o psb_zipcoo2csc.o \ psb_zcoins.o psb_zcsprt.o psb_zneigh.o psb_ztransp.o psb_ztransc.o\ psb_zrwextd.o psb_zsymbmm.o psb_znumbmm.o psb_zspscal.o psb_zspclip.o\ - psb_getifield.o psb_setifield.o psb_update_mod.o psb_getrow_mod.o + psb_getifield.o psb_setifield.o psb_update_mod.o psb_getrow_mod.o\ + psb_dgelp.o psb_zgelp.o LIBDIR=.. MODDIR=../modules diff --git a/base/tools/psb_dgelp.f90 b/base/serial/psb_dgelp.f90 similarity index 77% rename from base/tools/psb_dgelp.f90 rename to base/serial/psb_dgelp.f90 index ed32f990..a3bf37fb 100644 --- a/base/tools/psb_dgelp.f90 +++ b/base/serial/psb_dgelp.f90 @@ -30,9 +30,6 @@ !!$ ! File: psb_dgelp.f90 ! -! WARNING: This routine should be changed and moved to the serial part -! i.e. taking out the descriptor. -! ! ! Subroutine: psb_dgelp ! Apply a left permutation to a dense matrix @@ -42,24 +39,22 @@ ! iperm - integer. ! x - real, dimension(:,:). ! info - integer. Return code. -subroutine psb_dgelp(trans,iperm,x,desc_a,info) - use psb_descriptor_type - use psb_serial_mod +subroutine psb_dgelp(trans,iperm,x,info) + use psb_serial_mod, psb_protect_name => psb_dgelp use psb_const_mod - use psb_psblas_mod use psb_error_mod - use psb_penv_mod implicit none - type(psb_desc_type), intent(in) :: desc_a real(kind(1.d0)), intent(inout) :: x(:,:) - integer, intent(inout) :: iperm(:),info + integer, intent(in) :: iperm(:) + integer, intent(out) :: info character, intent(in) :: trans ! local variables integer :: ictxt,np, me,nrow,ncol - real(kind(1.d0)),pointer :: dtemp(:) + real(kind(1.d0)),allocatable :: dtemp(:) integer :: int_err(5), i1sz, i2sz, dectype, i, err_act + integer, allocatable :: itemp(:) real(kind(1.d0)),parameter :: one=1 logical, parameter :: debug=.false. @@ -88,47 +83,34 @@ subroutine psb_dgelp(trans,iperm,x,desc_a,info) info=0 call psb_erractionsave(err_act) - ictxt = psb_cd_get_context(desc_a) - dectype = psb_cd_get_dectype(desc_a) - nrow = psb_cd_get_local_rows(desc_a) - ncol = psb_cd_get_local_cols(desc_a) i1sz = size(x,dim=1) i2sz = size(x,dim=2) - call psb_info(ictxt, me, np) - - if (debug) write(*,*) 'asb start: ',np,me,& - & psb_cd_get_dectype(desc_a) - ! ....verify blacs grid correctness.. - if (np == -1) then - info = 2010 + if (debug) write(*,*) 'gelp: ',i1sz,i2sz + allocate(dtemp(i1sz),itemp(size(iperm)),stat=info) + if (info /= 0) then + info=2040 call psb_errpush(info,name) goto 9999 - else if (.not.psb_is_asb_desc(desc_a)) then - info = 3110 - call psb_errpush(info,name) - goto 9999 - endif - - - if (.not.isaperm(i1sz,iperm)) then + end if + itemp(:) = iperm(:) + + if (.not.isaperm(i1sz,itemp)) then info = 70 int_err(1) = 1 call psb_errpush(info,name,i_err=int_err) goto 9999 endif - if (debug) write(*,*) 'asb: ',i1sz,i2sz,nrow,ncol - allocate(dtemp(i1sz),stat=info) - call dgelp(trans,i1sz,i2sz,iperm,x,i1sz,dtemp,i1sz,info) + call dgelp(trans,i1sz,i2sz,itemp,x,i1sz,dtemp,i1sz,info) if(info /= 0) then info=4010 ch_err='dgelp' call psb_errpush(info,name,a_err=ch_err) end if - deallocate(dtemp) + deallocate(dtemp,itemp) call psb_erractionrestore(err_act) return @@ -178,9 +160,6 @@ end subroutine psb_dgelp !!$ !!$ ! -! WARNING: This routine should be changed and moved to the serial part -! i.e. taking out the descriptor. -! ! ! Subroutine: psb_dgelpv ! Apply a left permutation to a dense matrix @@ -190,24 +169,22 @@ end subroutine psb_dgelp ! iperm - integer. ! x - real, dimension(:). ! info - integer. Return code. -subroutine psb_dgelpv(trans,iperm,x,desc_a,info) - use psb_descriptor_type - use psb_serial_mod +subroutine psb_dgelpv(trans,iperm,x,info) + use psb_serial_mod, psb_protect_name => psb_dgelpv use psb_const_mod - use psb_psblas_mod use psb_error_mod - use psb_penv_mod implicit none - type(psb_desc_type), intent(in) :: desc_a real(kind(1.d0)), intent(inout) :: x(:) - integer, intent(inout) :: iperm(:), info + integer, intent(in) :: iperm(:) + integer, intent(out) :: info character, intent(in) :: trans ! local variables integer :: ictxt,np,me integer :: int_err(5), i1sz,nrow,ncol,dectype, err_act - real(kind(1.d0)),pointer :: dtemp(:) + real(kind(1.d0)),allocatable :: dtemp(:) + integer, allocatable :: itemp(:) real(kind(1.d0)),parameter :: one=1 logical, parameter :: debug=.false. @@ -238,43 +215,30 @@ subroutine psb_dgelpv(trans,iperm,x,desc_a,info) i1sz = size(x) - ictxt = psb_cd_get_context(desc_a) - dectype = psb_cd_get_dectype(desc_a) - nrow = psb_cd_get_local_rows(desc_a) - ncol = psb_cd_get_local_cols(desc_a) - - call psb_info(ictxt, me, np) - - ! ....verify blacs grid correctness.. - if (np == -1) then - info = 2010 + if (debug) write(*,*) 'gelp: ',i1sz + allocate(dtemp(i1sz),itemp(size(iperm)),stat=info) + if (info /= 0) then + info=2040 call psb_errpush(info,name) goto 9999 - else if (.not.psb_is_asb_desc(desc_a)) then - info = 3110 - call psb_errpush(info,name) - goto 9999 - endif - - if (debug) write(0,*) 'calling isaperm ',i1sz,size(iperm),trans - - if (.not.isaperm(i1sz,iperm)) then + end if + itemp(:) = iperm(:) + + if (.not.isaperm(i1sz,itemp)) then info = 70 int_err(1) = 1 call psb_errpush(info,name,i_err=int_err) goto 9999 endif - allocate(dtemp(i1sz),stat=info) - - call dgelp(trans,i1sz,1,iperm,x,i1sz,dtemp,i1sz,info) + call dgelp(trans,i1sz,1,itemp,x,i1sz,dtemp,i1sz,info) if(info /= 0) then info=4010 ch_err='dgelp' call psb_errpush(info,name,a_err=ch_err) end if - deallocate(dtemp) + deallocate(dtemp,itemp) call psb_erractionrestore(err_act) return diff --git a/base/tools/psb_zgelp.f90 b/base/serial/psb_zgelp.f90 similarity index 79% rename from base/tools/psb_zgelp.f90 rename to base/serial/psb_zgelp.f90 index da5a3c76..62418918 100644 --- a/base/tools/psb_zgelp.f90 +++ b/base/serial/psb_zgelp.f90 @@ -30,9 +30,6 @@ !!$ ! File: psb_zgelp.f90 ! -! WARNING: This routine should be changed and moved to the serial part -! i.e. taking out the descriptor. -! ! ! Subroutine: psb_zgelp ! Apply a left permutation to a dense matrix @@ -42,23 +39,21 @@ ! iperm - integer. ! x - real, dimension(:,:). ! info - integer. Return code. -subroutine psb_zgelp(trans,iperm,x,desc_a,info) - use psb_descriptor_type - use psb_serial_mod +subroutine psb_zgelp(trans,iperm,x,info) + use psb_serial_mod, psb_protect_name => psb_zgelp use psb_const_mod - use psb_psblas_mod use psb_error_mod - use psb_penv_mod implicit none - type(psb_desc_type), intent(in) :: desc_a complex(kind(1.d0)), intent(inout) :: x(:,:) - integer, intent(inout) :: iperm(:),info + integer, intent(in) :: iperm(:) + integer, intent(out) :: info character, intent(in) :: trans ! local variables integer :: ictxt,np,me,nrow,ncol - complex(kind(1.d0)),pointer :: dtemp(:) + complex(kind(1.d0)),allocatable :: dtemp(:) + integer, allocatable :: itemp(:) integer :: int_err(5), i1sz, i2sz, i, err_act real(kind(1.d0)),parameter :: one=1 logical, parameter :: debug=.false. @@ -88,46 +83,34 @@ subroutine psb_zgelp(trans,iperm,x,desc_a,info) info=0 call psb_erractionsave(err_act) - ictxt = psb_cd_get_context(desc_a) - nrow = psb_cd_get_local_rows(desc_a) - ncol = psb_cd_get_local_cols(desc_a) i1sz = size(x,dim=1) i2sz = size(x,dim=2) - call psb_info(ictxt, me, np) - - if (debug) write(*,*) 'asb start: ',np,me,& - & psb_cd_get_dectype(desc_a) - ! ....verify blacs grid correctness.. - if (np == -1) then - info = 2010 + if (debug) write(*,*) 'gelp: ',i1sz,i2sz + allocate(dtemp(i1sz),itemp(size(iperm)),stat=info) + if (info /= 0) then + info=2040 call psb_errpush(info,name) goto 9999 - else if (.not.psb_is_asb_desc(desc_a)) then - info = 3110 - call psb_errpush(info,name) - goto 9999 - endif - - - if (.not.isaperm(i1sz,iperm)) then + end if + itemp(:) = iperm(:) + + if (.not.isaperm(i1sz,itemp)) then info = 70 int_err(1) = 1 call psb_errpush(info,name,i_err=int_err) goto 9999 endif - if (debug) write(*,*) 'asb: ',i1sz,i2sz,nrow,ncol - allocate(dtemp(i1sz),stat=info) - call zgelp(trans,i1sz,i2sz,iperm,x,i1sz,dtemp,i1sz,info) + call zgelp(trans,i1sz,i2sz,itemp,x,i1sz,dtemp,i1sz,info) if(info.ne.0) then info=4010 ch_err='zgelp' call psb_errpush(info,name,a_err=ch_err) end if - deallocate(dtemp) + deallocate(dtemp,itemp) call psb_erractionrestore(err_act) return @@ -189,24 +172,22 @@ end subroutine psb_zgelp ! iperm - integer. ! x - real, dimension(:). ! info - integer. Return code. -subroutine psb_zgelpv(trans,iperm,x,desc_a,info) - use psb_descriptor_type - use psb_serial_mod +subroutine psb_zgelpv(trans,iperm,x,info) + use psb_serial_mod, psb_protect_name => psb_zgelpv use psb_const_mod - use psb_psblas_mod use psb_error_mod - use psb_penv_mod implicit none - type(psb_desc_type), intent(in) :: desc_a complex(kind(1.d0)), intent(inout) :: x(:) - integer, intent(inout) :: iperm(:), info + integer, intent(in) :: iperm(:) + integer, intent(out) :: info character, intent(in) :: trans ! local variables integer :: ictxt,np,me integer :: int_err(5), i1sz,nrow,ncol, i, err_act - complex(kind(1.d0)),pointer :: dtemp(:) + complex(kind(1.d0)),allocatable :: dtemp(:) + integer, allocatable :: itemp(:) real(kind(1.d0)),parameter :: one=1 logical, parameter :: debug=.false. @@ -237,35 +218,25 @@ subroutine psb_zgelpv(trans,iperm,x,desc_a,info) i1sz = size(x) - ictxt = psb_cd_get_context(desc_a) - nrow = psb_cd_get_local_rows(desc_a) - ncol = psb_cd_get_local_cols(desc_a) - - call psb_info(ictxt, me, np) - ! ....verify blacs grid correctness.. - if (np == -1) then - info = 2010 - call psb_errpush(info,name) - goto 9999 - else if (.not.psb_is_asb_desc(desc_a)) then - info = 3110 + if (debug) write(*,*) 'gelp: ',i1sz + allocate(dtemp(i1sz),itemp(size(iperm)),stat=info) + if (info /= 0) then + info=2040 call psb_errpush(info,name) goto 9999 - endif - - if (debug) write(0,*) 'calling isaperm ',i1sz,size(iperm),trans - - if (.not.isaperm(i1sz,iperm)) then + end if + itemp(:) = iperm(:) + + if (.not.isaperm(i1sz,itemp)) then info = 70 int_err(1) = 1 call psb_errpush(info,name,i_err=int_err) goto 9999 endif - allocate(dtemp(i1sz),stat=info) - call zgelp(trans,i1sz,1,iperm,x,i1sz,dtemp,i1sz,info) + call zgelp(trans,i1sz,1,itemp,x,i1sz,dtemp,i1sz,info) if(info.ne.0) then info=4010 ch_err='zgelp' diff --git a/base/tools/Makefile b/base/tools/Makefile index f7819d50..799d7595 100644 --- a/base/tools/Makefile +++ b/base/tools/Makefile @@ -1,7 +1,7 @@ include ../../Make.inc FOBJS = psb_dallc.o psb_dasb.o psb_dcsrp.o psb_cdprt.o \ - psb_dfree.o psb_dgelp.o psb_dins.o \ + psb_dfree.o psb_dins.o \ psb_cdals.o psb_cdalv.o psb_cd_inloc.o psb_cdcpy.o \ psb_cdfree.o psb_cdins.o \ psb_cdren.o psb_cdrep.o psb_cdtransfer.o psb_get_overlap.o\ @@ -11,7 +11,7 @@ FOBJS = psb_dallc.o psb_dasb.o psb_dcsrp.o psb_cdprt.o \ psb_ifree.o psb_iins.o psb_loc_to_glob.o\ psb_zallc.o psb_zasb.o psb_zfree.o psb_zins.o \ psb_zspalloc.o psb_zspasb.o psb_zspfree.o\ - psb_zspins.o psb_zsprn.o psb_zcdovr.o psb_zgelp.o + psb_zspins.o psb_zsprn.o psb_zcdovr.o MPFOBJS = psb_dsphalo.o psb_zsphalo.o psb_icdasb.o psb_dcdovr.o psb_zcdovr.o diff --git a/base/tools/psb_cd_inloc.f90 b/base/tools/psb_cd_inloc.f90 index d61bc15d..397784e9 100644 --- a/base/tools/psb_cd_inloc.f90 +++ b/base/tools/psb_cd_inloc.f90 @@ -31,13 +31,13 @@ ! File: psb_cdalv.f90 ! ! Subroutine: psb_cdalv -! Allocate descriptor -! and checks correctness of PARTS subroutine +! Allocate descriptor with a local vector V containing the list +! of indices that are assigned to the current process. The global size +! is equal to the largest index found on any process. ! ! Parameters: -! m - integer. The number of rows. ! v - integer, dimension(:). The array containg the partitioning scheme. -! ictxt - integer. The communication context. +! ictxt - integer. The communication context. ! desc_a - type(). The communication descriptor. ! info - integer. Eventually returns an error code subroutine psb_cd_inloc(v, ictxt, desc_a, info) diff --git a/base/tools/psb_cdals.f90 b/base/tools/psb_cdals.f90 index 2da3be41..a78eef14 100644 --- a/base/tools/psb_cdals.f90 +++ b/base/tools/psb_cdals.f90 @@ -37,7 +37,8 @@ ! Parameters: ! m - integer. The number of rows. ! n - integer. The number of columns. -! parts - external subroutine. The routine that contains the partitioning scheme. +! parts - external subroutine. The routine that contains the +! partitioning scheme. ! ictxt - integer. The communication context. ! desc_a - type(). The communication descriptor. ! info - integer. Error code (if any). diff --git a/base/tools/psb_cdalv.f90 b/base/tools/psb_cdalv.f90 index c22378ae..f6ad2fa8 100644 --- a/base/tools/psb_cdalv.f90 +++ b/base/tools/psb_cdalv.f90 @@ -31,16 +31,17 @@ ! File: psb_cdalv.f90 ! ! Subroutine: psb_cdalv -! Allocate descriptor -! and checks correctness of PARTS subroutine +! Allocate descriptor by means of a global map vector V: index I +! is assigned to process V(I). It is assumed that V is identical +! on all calling processes. +! ! ! Parameters: -! m - integer. The number of rows. ! v - integer, dimension(:). The array containg the partitioning scheme. -! ictxt - integer. The communication context. +! ictxt - integer. The communication context. ! desc_a - type(). The communication descriptor. -! info - integer. Eventually returns an error code -! flag - integer. ??? +! info - integer. Return code +! flag - integer. Are V's contents 0- or 1-based? subroutine psb_cdalv(v, ictxt, desc_a, info, flag) use psb_descriptor_type use psb_serial_mod diff --git a/base/tools/psb_cdcpy.f90 b/base/tools/psb_cdcpy.f90 index 9b0efcb2..ca3c7a36 100644 --- a/base/tools/psb_cdcpy.f90 +++ b/base/tools/psb_cdcpy.f90 @@ -36,7 +36,7 @@ ! Parameters: ! desc_in - type(). The communication descriptor to be cloned. ! desc_out - type(). The output communication descriptor. -! info - integer. Eventually returns an error code. +! info - integer. Return code. subroutine psb_cdcpy(desc_in, desc_out, info) use psb_descriptor_type diff --git a/base/tools/psb_cdins.f90 b/base/tools/psb_cdins.f90 index fac22403..ceb6a291 100644 --- a/base/tools/psb_cdins.f90 +++ b/base/tools/psb_cdins.f90 @@ -38,7 +38,10 @@ ! ia - integer,dimension(:). The row indices of the points. ! ja - integer,dimension(:). The column indices of the points. ! desc_a - type(). The communication descriptor to be freed. -! info - integer. Eventually returns an error code. +! info - integer. Return code. +! ila - integer,dimension(:). The row indices in local numbering +! jla - integer,dimension(:). The col indices in local numbering +! subroutine psb_cdins(nz,ia,ja,desc_a,info,ila,jla) use psb_descriptor_type diff --git a/base/tools/psb_cdprt.f90 b/base/tools/psb_cdprt.f90 index 0c58a33b..d895ee8c 100644 --- a/base/tools/psb_cdprt.f90 +++ b/base/tools/psb_cdprt.f90 @@ -38,6 +38,7 @@ ! desc_p - type(). The communication descriptor to be printed. ! glob - logical(otpional). Wheter to print out global or local data. ! short - logical(optional). Used to choose a verbose output. +! subroutine psb_cdprt(iout,desc_p,glob,short) use psb_const_mod use psb_descriptor_type diff --git a/base/tools/psb_cdren.f90 b/base/tools/psb_cdren.f90 index cd0eafec..b25e8d74 100644 --- a/base/tools/psb_cdren.f90 +++ b/base/tools/psb_cdren.f90 @@ -29,14 +29,18 @@ !!$ !!$ ! File: psb_cdren.f90 -! +! +! WARNING: this routine is almost certainly obsolete. Must be reviewed. +! ! Subroutine: psb_cdren ! Updates a communication descriptor according to a renumbering scheme. ! ! Parameters: -! trans - character. Whether iperm or its transpose should be applied. +! trans - character. Whether iperm or its transpose +! should be applied. ! iperm - integer,dimension(:). The renumbering scheme. -! desc_a - type(). The communication descriptor to be updated. +! desc_a - type(). The communication descriptor +! to be updated. ! info - integer. Eventually returns an error code. ! subroutine psb_cdren(trans,iperm,desc_a,info) diff --git a/base/tools/psb_cdtransfer.f90 b/base/tools/psb_cdtransfer.f90 index dead1f7a..c950b2aa 100644 --- a/base/tools/psb_cdtransfer.f90 +++ b/base/tools/psb_cdtransfer.f90 @@ -34,9 +34,10 @@ ! Transfers data and allocation from in to out (just like MOVE_ALLOC). ! ! Parameters: -! desc_in - type(). The communication descriptor to be transferred. +! desc_in - type(). The communication descriptor to be +! transferred. ! desc_out - type(). The output communication descriptor. -! info - integer. Eventually returns an error code. +! info - integer. Return code. subroutine psb_cdtransfer(desc_in, desc_out, info) use psb_descriptor_type @@ -75,18 +76,31 @@ subroutine psb_cdtransfer(desc_in, desc_out, info) endif call psb_transfer( desc_in%matrix_data , desc_out%matrix_data , info) - if (info == 0) call psb_transfer( desc_in%halo_index , desc_out%halo_index , info) - if (info == 0) call psb_transfer( desc_in%bnd_elem , desc_out%bnd_elem , info) - if (info == 0) call psb_transfer( desc_in%ovrlap_elem , desc_out%ovrlap_elem , info) - if (info == 0) call psb_transfer( desc_in%ovrlap_index, desc_out%ovrlap_index , info) - if (info == 0) call psb_transfer( desc_in%ext_index , desc_out%ext_index , info) - if (info == 0) call psb_transfer( desc_in%loc_to_glob , desc_out%loc_to_glob , info) - if (info == 0) call psb_transfer( desc_in%glob_to_loc , desc_out%glob_to_loc , info) - if (info == 0) call psb_transfer( desc_in%lprm , desc_out%lprm , info) - if (info == 0) call psb_transfer( desc_in%idx_space , desc_out%idx_space , info) - if (info == 0) call psb_transfer( desc_in%hashv , desc_out%hashv , info) - if (info == 0) call psb_transfer( desc_in%glb_lc , desc_out%glb_lc , info) - if (info == 0) call psb_transfer( desc_in%ptree , desc_out%ptree , info) + if (info == 0) & + & call psb_transfer( desc_in%halo_index , desc_out%halo_index , info) + if (info == 0) & + & call psb_transfer( desc_in%bnd_elem , desc_out%bnd_elem , info) + if (info == 0) & + & call psb_transfer( desc_in%ovrlap_elem , desc_out%ovrlap_elem , info) + if (info == 0) & + & call psb_transfer( desc_in%ovrlap_index, desc_out%ovrlap_index , info) + if (info == 0) & + & call psb_transfer( desc_in%ext_index , desc_out%ext_index , info) + if (info == 0) & + & call psb_transfer( desc_in%loc_to_glob , desc_out%loc_to_glob , info) + if (info == 0) & + & call psb_transfer( desc_in%glob_to_loc , desc_out%glob_to_loc , info) + if (info == 0) & + & call psb_transfer( desc_in%lprm , desc_out%lprm , info) + if (info == 0) & + & call psb_transfer( desc_in%idx_space , desc_out%idx_space , info) + if (info == 0) & + & call psb_transfer( desc_in%hashv , desc_out%hashv , info) + if (info == 0) & + & call psb_transfer( desc_in%glb_lc , desc_out%glb_lc , info) + if (info == 0) & + & call psb_transfer( desc_in%ptree , desc_out%ptree , info) + if (info /= 0) then info = 4010 call psb_errpush(info,name) diff --git a/base/tools/psb_dcdovr.F90 b/base/tools/psb_dcdovr.F90 index 1e9c7ce0..be12a616 100644 --- a/base/tools/psb_dcdovr.F90 +++ b/base/tools/psb_dcdovr.F90 @@ -33,8 +33,7 @@ ! Subroutine: psb_cdovr ! This routine takes a matrix A with its descriptor, and builds the ! auxiliary descriptor corresponding to the number of overlap levels -! specified on input. It really is just a size estimation/allocation -! front end for . +! specified on input. ! ! Parameters: ! a - type(). The input sparse matrix. @@ -42,7 +41,22 @@ ! novr - integer. The number of overlap levels. ! desc_ov - type(). The auxiliary output communication ! descriptor. -! info - integer. Eventually returns an error code. +! info - integer. Return code. +! extype - integer. Choice of type of overlap: +! psb_ovt_xhal_: build a descriptor with an extended +! stencil, i.e. enlarge the existing +! halo by novr additional layers. +! psb_ovt_asov_: build a descriptor suitable +! for Additive Schwarz preconditioner. +! This last choice implies that: +! a. The novr halo layers are added to +! the overlap; +! b. The novr layers are also copied to +! the ext_ structure to provide +! the mapping between the base +! descriptor and the overlapped one. +! c. The (novr+1) layer becomes the new +! halo. ! Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) @@ -670,7 +684,7 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) ! Build an overlapped descriptor for Additive Schwarz ! with overlap enlargement; we need the overlap, ! the (new) halo and the mapping into the new index space. - ! Here we need: 1. (orig_ovr + t_ovr_idx -> ovrlap + ! Here we need: 1. (orig_ovr + t_ovr_idx) -> ovrlap ! 2. (tmp_halo) -> ext ! 3. (t_halo_in) -> halo ! 4. n_row(ov) current. diff --git a/base/tools/psb_dcsrp.f90 b/base/tools/psb_dcsrp.f90 index 92872650..292586f0 100644 --- a/base/tools/psb_dcsrp.f90 +++ b/base/tools/psb_dcsrp.f90 @@ -30,13 +30,18 @@ !!$ ! File: psb_dcsrp.f90 ! +! WARNING: This routine should be changed and moved to the serial part +! i.e. taking out the descriptor. +! ! Subroutine: psb_dcsrp ! Apply a right permutation to a sparse matrix, i.e. permute the column ! indices. ! ! Parameters: -! trans - character. Whether iperm or its transpose should be applied -! iperm - integer, pointer, dimension(:). A permutation vector; its size must be either N_ROW or N_COL +! trans - character. Whether iperm or its transpose +! should be applied +! iperm - integer, dimension(:) A permutation vector; its size +! must be either N_ROW or N_COL ! a - type(). The communication descriptor. ! info - integer. Eventually returns an error code diff --git a/base/tools/psb_dfree.f90 b/base/tools/psb_dfree.f90 index f1c4d605..f1cd3340 100644 --- a/base/tools/psb_dfree.f90 +++ b/base/tools/psb_dfree.f90 @@ -34,9 +34,9 @@ ! frees a dense matrix structure ! ! Parameters: -! x - real, pointer, dimension(:,:). The dense matrix to be freed. +! x - real, allocatable, dimension(:,:). The dense matrix to be freed. ! desc_a - type(). The communication descriptor. -! info - integer. Eventually returns an error code +! info - integer. Return code subroutine psb_dfree(x, desc_a, info) !...free dense matrix structure... use psb_const_mod @@ -109,9 +109,9 @@ end subroutine psb_dfree ! frees a dense matrix structure ! ! Parameters: -! x - real, pointer, dimension(:). The dense matrix to be freed. +! x - real, allocatable, dimension(:). The dense matrix to be freed. ! desc_a - type(). The communication descriptor. -! info - integer. Eventually returns an error code +! info - integer. Return code subroutine psb_dfreev(x, desc_a, info) !...free dense matrix structure... use psb_const_mod diff --git a/base/tools/psb_dsphalo.F90 b/base/tools/psb_dsphalo.F90 index 38ed63d2..a66f5f83 100644 --- a/base/tools/psb_dsphalo.F90 +++ b/base/tools/psb_dsphalo.F90 @@ -30,16 +30,32 @@ !!$ ! File: psb_dsphalo.f90 ! -!***************************************************************************** -!* * -!* This routine does the retrieval of remote matrix rows. * -!* Note that retrieval is done through GTBLK, therefore it should work * -!* for any format. * -!* * -!* * -!* * -!***************************************************************************** -Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) +! Subroutine: psb_dsphalo +! This routine does the retrieval of remote matrix rows. +! Note that retrieval is done through GTBLK, therefore it should work +! for any matrix format in A; as for the output, default is CSR. +! +! +! Parameters: +! a - type(psb_dspmat_type) The local part of input matrix A +! desc_a - type(). The communication descriptor. +! blck - type(psb_dspmat_type) The local part of output matrix BLCK +! info - integer. Return code +! rowcnv - logical Should row/col indices be converted +! colcnv - logical to/from global numbering when sent/received? +! default is .TRUE. +! rowscale - logical Should row/col indices on output be remapped +! colscale - logical from MIN:MAX to 1:(MAX-MIN+1) ? +! default is .FALSE. +! (commmon use is ROWSCALE=.TRUE., COLSCALE=.FALSE.) +! data - integer Which index list in desc_a should be used to retrieve +! rows, default psb_comm_halo_ (i.e.: use halo_index) +! other value psb_comm_ext_, no longer accepting +! psb_comm_ovrl_, perhaps should be reinstated in the future. +! +! +Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& + & rowscale,colscale,outfmt,data) use psb_const_mod use psb_serial_mod @@ -60,20 +76,20 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) Type(psb_dspmat_type),Intent(inout) :: blk Type(psb_desc_type),Intent(in), target :: desc_a integer, intent(out) :: info - logical, optional, intent(in) :: rwcnv,clcnv,cliprow + logical, optional, intent(in) :: rowcnv,colcnv,rowscale,colscale character(len=5), optional :: outfmt integer, intent(in), optional :: data ! ...local scalars.... Integer :: np,me,counter,proc,i, & & n_el_send,k,n_el_recv,ictxt, idx, r, tot_elem,& & n_elem, j, ipx,mat_recv, iszs, iszr,idxs,idxr,nz,& - & nrmin,data_,ngtz + & irmin,icmin,irmax,icmax,data_,ngtz Integer :: l1, icomm, err_act Integer, allocatable :: sdid(:,:), brvindx(:),rvid(:,:), & & rvsz(:), bsdindx(:),sdsz(:), iasnd(:), jasnd(:) real(kind(1.d0)), allocatable :: valsnd(:) integer, pointer :: idxv(:) - logical :: rwcnv_,clcnv_,cliprow_ + logical :: rowcnv_,colcnv_,rowscale_,colscale_ character(len=5) :: outfmt_ Logical,Parameter :: debug=.false., debugprt=.false. real(kind(1.d0)) :: t1,t2,t3,t4,t5,t6,t7,t8,t9 @@ -85,20 +101,25 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) call psb_erractionsave(err_act) if(debug) write(0,*)'Inside DSPHALO' - if (present(rwcnv)) then - rwcnv_ = rwcnv + if (present(rowcnv)) then + rowcnv_ = rowcnv + else + rowcnv_ = .true. + endif + if (present(colcnv)) then + colcnv_ = colcnv else - rwcnv_ = .true. + colcnv_ = .true. endif - if (present(clcnv)) then - clcnv_ = clcnv + if (present(rowscale)) then + rowscale_ = rowscale else - clcnv_ = .true. + rowscale_ = .false. endif - if (present(cliprow)) then - cliprow_ = cliprow + if (present(colscale)) then + colscale_ = colscale else - cliprow_ = .false. + colscale_ = .false. endif if (present(data)) then data_ = data @@ -133,12 +154,12 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) case(psb_comm_halo_) idxv => desc_a%halo_index - case(psb_comm_ovr_) - idxv => desc_a%ovrlap_index - case(psb_comm_ext_) idxv => desc_a%ext_index +!!$ case(psb_comm_ovr_) +!!$ idxv => desc_a%ovrlap_index +!!$ ! Do not accept OVRLAP_INDEX any longer. case default call psb_errpush(4010,name,a_err='wrong Data selector') goto 9999 @@ -261,8 +282,8 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) Enddo nz = tot_elem - if (rwcnv_) call psb_loc_to_glob(iasnd(1:nz),desc_a,info,iact='I') - if (clcnv_) call psb_loc_to_glob(jasnd(1:nz),desc_a,info,iact='I') + if (rowcnv_) call psb_loc_to_glob(iasnd(1:nz),desc_a,info,iact='I') + if (colcnv_) call psb_loc_to_glob(jasnd(1:nz),desc_a,info,iact='I') if (info /= 0) then info=4010 ch_err='psb_loc_to_glob' @@ -290,8 +311,8 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) ! ! Convert into local numbering ! - if (rwcnv_) call psb_glob_to_loc(blk%ia1(1:iszr),desc_a,info,iact='I',owned=cliprow_) - if (clcnv_) call psb_glob_to_loc(blk%ia2(1:iszr),desc_a,info,iact='I') + if (rowcnv_) call psb_glob_to_loc(blk%ia1(1:iszr),desc_a,info,iact='I') + if (colcnv_) call psb_glob_to_loc(blk%ia2(1:iszr),desc_a,info,iact='I') if (info /= 0) then info=4010 @@ -309,23 +330,42 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) end if l1 = 0 blk%m=0 - nrmin=max(0,a%m) + ! + irmin = huge(irmin) + icmin = huge(icmin) + irmax = 0 + icmax = 0 Do i=1,iszr -!!$ write(0,*) work5(i),work6(i) r=(blk%ia1(i)) k=(blk%ia2(i)) - If ((r>nrmin).and.(k>0)) Then + ! Just in case some of the conversions were out-of-range + If ((r>0).and.(k>0)) Then l1=l1+1 blk%aspk(l1) = blk%aspk(i) - blk%ia1(l1) = r + blk%ia1(l1) = r blk%ia2(l1) = k - blk%k = max(blk%k,k) - blk%m = max(blk%m,r) + irmin = min(irmin,r) + irmax = max(irmax,r) + icmin = min(icmin,k) + icmax = max(icmax,k) End If Enddo - blk%fida='COO' - blk%infoa(psb_nnz_)=l1 - blk%m = blk%m - a%m + if (rowscale_) then + blk%m = irmax-irmin+1 + blk%ia1(1:l1) = blk%ia1(1:l1) - irmin + 1 + else + blk%m = irmax + end if + if (colscale_) then + blk%k = icmax-icmin+1 + blk%ia2(1:l1) = blk%ia2(1:l1) - icmin + 1 + else + blk%k = icmax + end if + + blk%fida = 'COO' + blk%infoa(psb_nnz_) = l1 + if (debugprt) then open(50+me) call psb_csprt(50+me,blk,head='% SPHALO border .') @@ -336,32 +376,16 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) if(debug) Write(0,*)me,'End first loop',counter,l1,blk%m - ! - ! Combined sort & conversion to CSR. - ! - if(debug) write(0,*) me,'Calling ipcoo2csr from dsphalo ',blk%m,blk%k,l1,blk%ia2(2) - - select case(outfmt_) - case ('CSR') - call psb_ipcoo2csr(blk,info,rwshr=.true.) - if (info /= 0) then - info=4010 - ch_err='psb_ipcoo2csr' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - case('COO') - ! Do nothing! - case default - write(0,*) 'Error in DSPHALO : invalid outfmt "',outfmt_,'"' + ! Do we expect any duplicates to appear???? + call psb_spcnv(blk,info,afmt=outfmt_,dupl=psb_dupl_add_) + if (info /= 0) then info=4010 - ch_err='Bad outfmt' + ch_err='psb_spcnv' call psb_errpush(info,name,a_err=ch_err) goto 9999 - end select - t5 = psb_wtime() - + end if + t5 = psb_wtime() !!$ write(0,'(i3,1x,a,4(1x,i14))') me,'DSPHALO sizes:',iszr,iszs !!$ write(0,'(i3,1x,a,4(1x,g14.5))') me,'DSPHALO timings:',t6-t2,t7-t6,t8-t7,t3-t8 diff --git a/base/tools/psb_zcdovr.F90 b/base/tools/psb_zcdovr.F90 index a48e7823..5ee59750 100644 --- a/base/tools/psb_zcdovr.F90 +++ b/base/tools/psb_zcdovr.F90 @@ -33,15 +33,30 @@ ! Subroutine: psb_zcdovr ! This routine takes a matrix A with its descriptor, and builds the ! auxiliary descriptor corresponding to the number of overlap levels -! specified on input. It really is just a size estimation/allocation -! front end for . +! specified on input. ! ! Parameters: ! a - type(). The input sparse matrix. ! desc_a - type(). The input communication descriptor. -! norv - integer. The number of overlap levels. -! desc_ov - type(). The auxiliary output communication descriptor. -! info - integer. Eventually returns an error code. +! novr - integer. The number of overlap levels. +! desc_ov - type(). The auxiliary output communication +! descriptor. +! info - integer. Return code. +! extype - integer. Choice of type of overlap: +! psb_ovt_xhal_: build a descriptor with an extended +! stencil, i.e. enlarge the existing +! halo by novr additional layers. +! psb_ovt_asov_: build a descriptor suitable +! for Additive Schwarz preconditioner. +! This last choice implies that: +! a. The novr halo layers are added to +! the overlap; +! b. The novr layers are also copied to +! the ext_ structure to provide +! the mapping between the base +! descriptor and the overlapped one. +! c. The (novr+1) layer becomes the new +! halo. ! Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) @@ -668,7 +683,7 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) ! Build an overlapped descriptor for Additive Schwarz ! with overlap enlargement; we need the overlap, ! the (new) halo and the mapping into the new index space. - ! Here we need: 1. (orig_ovr + t_ovr_idx -> ovrlap + ! Here we need: 1. (orig_ovr + t_ovr_idx) -> ovrlap ! 2. (tmp_halo) -> ext ! 3. (t_halo_in) -> halo ! 4. n_row(ov) current. diff --git a/base/tools/psb_zcsrp.f90 b/base/tools/psb_zcsrp.f90 index 115c125e..f2c5d192 100644 --- a/base/tools/psb_zcsrp.f90 +++ b/base/tools/psb_zcsrp.f90 @@ -30,13 +30,18 @@ !!$ ! File: psb_zcsrp.f90 ! +! WARNING: This routine should be changed and moved to the serial part +! i.e. taking out the descriptor. +! ! Subroutine: psb_zcsrp ! Apply a right permutation to a sparse matrix, i.e. permute the column ! indices. ! ! Parameters: -! trans - character. Whether iperm or its transpose should be applied -! iperm - integer, pointer, dimension(:). A permutation vector; its size must be either N_ROW or N_COL +! trans - character. Whether iperm or its transpose +! should be applied +! iperm - integer, dimension(:) A permutation vector; its size +! must be either N_ROW or N_COL ! a - type(). The communication descriptor. ! info - integer. Eventually returns an error code diff --git a/base/tools/psb_zfree.f90 b/base/tools/psb_zfree.f90 index cda58bb6..eb637423 100644 --- a/base/tools/psb_zfree.f90 +++ b/base/tools/psb_zfree.f90 @@ -34,9 +34,9 @@ ! frees a dense matrix structure ! ! Parameters: -! x - real, pointer, dimension(:,:). The dense matrix to be freed. +! x - real, allocatable, dimension(:,:). The dense matrix to be freed. ! desc_a - type(). The communication descriptor. -! info - integer. Eventually returns an error code +! info - integer. Return code subroutine psb_zfree(x, desc_a, info) !...free dense matrix structure... use psb_const_mod @@ -109,9 +109,9 @@ end subroutine psb_zfree ! frees a dense matrix structure ! ! Parameters: -! x - real, pointer, dimension(:). The dense matrix to be freed. +! x - real, allocatable, dimension(:). The dense matrix to be freed. ! desc_a - type(). The communication descriptor. -! info - integer. Eventually returns an error code +! info - integer. Return code subroutine psb_zfreev(x, desc_a, info) !...free dense matrix structure... use psb_const_mod diff --git a/base/tools/psb_zsphalo.F90 b/base/tools/psb_zsphalo.F90 index 219a92e7..e699ca73 100644 --- a/base/tools/psb_zsphalo.F90 +++ b/base/tools/psb_zsphalo.F90 @@ -30,16 +30,32 @@ !!$ ! File: psb_zsphalo.f90 ! -!***************************************************************************** -!* * -!* This routine does the retrieval of remote matrix rows. * -!* Note that retrieval is done through GTBLK, therefore it should work * -!* for any format. * -!* * -!* * -!* * -!***************************************************************************** -Subroutine psb_zsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) +! Subroutine: psb_zsphalo +! This routine does the retrieval of remote matrix rows. +! Note that retrieval is done through GTBLK, therefore it should work +! for any matrix format in A; as for the output, default is CSR. +! +! +! Parameters: +! a - type(psb_zspmat_type) The local part of input matrix A +! desc_a - type(). The communication descriptor. +! blck - type(psb_zspmat_type) The local part of output matrix BLCK +! info - integer. Return code +! rowcnv - logical Should row/col indices be converted +! colcnv - logical to/from global numbering when sent/received? +! default is .TRUE. +! rowscale - logical Should row/col indices on output be remapped +! colscale - logical from MIN:MAX to 1:(MAX-MIN+1) ? +! default is .FALSE. +! (commmon use is ROWSCALE=.TRUE., COLSCALE=.FALSE.) +! data - integer Which index list in desc_a should be used to retrieve +! rows, default psb_comm_halo_ (i.e.: use halo_index) +! other value psb_comm_ext_, no longer accepting +! psb_comm_ovrl_, perhaps should be reinstated in the future. +! +! +Subroutine psb_zsphalo(a,desc_a,blk,info,rowcnv,colcnv,& + & rowscale,colscale,outfmt,data) use psb_const_mod use psb_serial_mod @@ -60,20 +76,20 @@ Subroutine psb_zsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) Type(psb_zspmat_type),Intent(inout) :: blk Type(psb_desc_type),Intent(in), target :: desc_a integer, intent(out) :: info - logical, optional, intent(in) :: rwcnv,clcnv,cliprow + logical, optional, intent(in) :: rowcnv,colcnv,rowscale,colscale character(len=5), optional :: outfmt integer, intent(in), optional :: data ! ...local scalars.... Integer :: np,me,counter,proc,i, & & n_el_send,k,n_el_recv,ictxt, idx, r, tot_elem,& & n_elem, j, ipx,mat_recv, iszs, iszr,idxs,idxr,nz,& - & nrmin,data_,ngtz + & irmin,icmin,irmax,icmax,data_,ngtz Integer :: l1, icomm, err_act Integer, allocatable :: sdid(:,:), brvindx(:),rvid(:,:), & & rvsz(:), bsdindx(:),sdsz(:), iasnd(:), jasnd(:) complex(kind(1.d0)), allocatable :: valsnd(:) integer, pointer :: idxv(:) - logical :: rwcnv_,clcnv_,cliprow_ + logical :: rowcnv_,colcnv_,rowscale_,colscale_ character(len=5) :: outfmt_ Logical,Parameter :: debug=.false., debugprt=.false. real(kind(1.d0)) :: t1,t2,t3,t4,t5,t6,t7,t8,t9 @@ -85,20 +101,25 @@ Subroutine psb_zsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) call psb_erractionsave(err_act) if(debug) write(0,*)'Inside DSPHALO' - if (present(rwcnv)) then - rwcnv_ = rwcnv + if (present(rowcnv)) then + rowcnv_ = rowcnv + else + rowcnv_ = .true. + endif + if (present(colcnv)) then + colcnv_ = colcnv else - rwcnv_ = .true. + colcnv_ = .true. endif - if (present(clcnv)) then - clcnv_ = clcnv + if (present(rowscale)) then + rowscale_ = rowscale else - clcnv_ = .true. + rowscale_ = .false. endif - if (present(cliprow)) then - cliprow_ = cliprow + if (present(colscale)) then + colscale_ = colscale else - cliprow_ = .false. + colscale_ = .false. endif if (present(data)) then data_ = data @@ -133,12 +154,12 @@ Subroutine psb_zsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) case(psb_comm_halo_) idxv => desc_a%halo_index - case(psb_comm_ovr_) - idxv => desc_a%ovrlap_index - case(psb_comm_ext_) idxv => desc_a%ext_index +!!$ case(psb_comm_ovr_) +!!$ idxv => desc_a%ovrlap_index +!!$ ! Do not accept OVRLAP_INDEX any longer. case default call psb_errpush(4010,name,a_err='wrong Data selector') goto 9999 @@ -261,8 +282,8 @@ Subroutine psb_zsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) Enddo nz = tot_elem - if (rwcnv_) call psb_loc_to_glob(iasnd(1:nz),desc_a,info,iact='I') - if (clcnv_) call psb_loc_to_glob(jasnd(1:nz),desc_a,info,iact='I') + if (rowcnv_) call psb_loc_to_glob(iasnd(1:nz),desc_a,info,iact='I') + if (colcnv_) call psb_loc_to_glob(jasnd(1:nz),desc_a,info,iact='I') if (info /= 0) then info=4010 ch_err='psb_loc_to_glob' @@ -290,8 +311,8 @@ Subroutine psb_zsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) ! ! Convert into local numbering ! - if (rwcnv_) call psb_glob_to_loc(blk%ia1(1:iszr),desc_a,info,iact='I',owned=cliprow_) - if (clcnv_) call psb_glob_to_loc(blk%ia2(1:iszr),desc_a,info,iact='I') + if (rowcnv_) call psb_glob_to_loc(blk%ia1(1:iszr),desc_a,info,iact='I') + if (colcnv_) call psb_glob_to_loc(blk%ia2(1:iszr),desc_a,info,iact='I') if (info /= 0) then info=4010 @@ -309,23 +330,42 @@ Subroutine psb_zsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) end if l1 = 0 blk%m=0 - nrmin=max(0,a%m) + ! + irmin = huge(irmin) + icmin = huge(icmin) + irmax = 0 + icmax = 0 Do i=1,iszr -!!$ write(0,*) work5(i),work6(i) r=(blk%ia1(i)) k=(blk%ia2(i)) - If ((r>nrmin).and.(k>0)) Then + ! Just in case some of the conversions were out-of-range + If ((r>0).and.(k>0)) Then l1=l1+1 blk%aspk(l1) = blk%aspk(i) - blk%ia1(l1) = r + blk%ia1(l1) = r blk%ia2(l1) = k - blk%k = max(blk%k,k) - blk%m = max(blk%m,r) + irmin = min(irmin,r) + irmax = max(irmax,r) + icmin = min(icmin,k) + icmax = max(icmax,k) End If Enddo - blk%fida='COO' - blk%infoa(psb_nnz_)=l1 - blk%m = blk%m - a%m + if (rowscale_) then + blk%m = irmax-irmin+1 + blk%ia1(1:l1) = blk%ia1(1:l1) - irmin + 1 + else + blk%m = irmax + end if + if (colscale_) then + blk%k = icmax-icmin+1 + blk%ia2(1:l1) = blk%ia2(1:l1) - icmin + 1 + else + blk%k = icmax + end if + + blk%fida = 'COO' + blk%infoa(psb_nnz_) = l1 + if (debugprt) then open(50+me) call psb_csprt(50+me,blk,head='% SPHALO border .') @@ -336,32 +376,16 @@ Subroutine psb_zsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) if(debug) Write(0,*)me,'End first loop',counter,l1,blk%m - ! - ! Combined sort & conversion to CSR. - ! - if(debug) write(0,*) me,'Calling ipcoo2csr from dsphalo ',blk%m,blk%k,l1,blk%ia2(2) - - select case(outfmt_) - case ('CSR') - call psb_ipcoo2csr(blk,info,rwshr=.true.) - if (info /= 0) then - info=4010 - ch_err='psb_ipcoo2csr' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - case('COO') - ! Do nothing! - case default - write(0,*) 'Error in DSPHALO : invalid outfmt "',outfmt_,'"' + ! Do we expect any duplicates to appear???? + call psb_spcnv(blk,info,afmt=outfmt_,dupl=psb_dupl_add_) + if (info /= 0) then info=4010 - ch_err='Bad outfmt' + ch_err='psb_spcnv' call psb_errpush(info,name,a_err=ch_err) goto 9999 - end select - t5 = psb_wtime() - + end if + t5 = psb_wtime() !!$ write(0,'(i3,1x,a,4(1x,i14))') me,'DSPHALO sizes:',iszr,iszs !!$ write(0,'(i3,1x,a,4(1x,g14.5))') me,'DSPHALO timings:',t6-t2,t7-t6,t8-t7,t3-t8