From afecc0a1ede53b5d7ace909173de30570d77c44a Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Sat, 3 Dec 2011 11:22:34 +0000 Subject: [PATCH] psblas-testpre: base/serial/impl/psb_c_coo_impl.f90 base/serial/impl/psb_c_mat_impl.F90 base/serial/impl/psb_d_coo_impl.f90 base/serial/impl/psb_d_mat_impl.F90 base/serial/impl/psb_s_coo_impl.f90 base/serial/impl/psb_s_mat_impl.F90 base/serial/impl/psb_z_coo_impl.f90 base/serial/impl/psb_z_mat_impl.F90 base/tools/psb_csphalo.F90 base/tools/psb_dsphalo.F90 base/tools/psb_ssphalo.F90 base/tools/psb_zsphalo.F90 Reverted fix due to bug in GNU fortran. --- base/serial/impl/psb_c_coo_impl.f90 | 29 ++- base/serial/impl/psb_c_mat_impl.F90 | 28 +-- base/serial/impl/psb_d_coo_impl.f90 | 29 ++- base/serial/impl/psb_d_mat_impl.F90 | 28 +-- base/serial/impl/psb_s_coo_impl.f90 | 29 ++- base/serial/impl/psb_s_mat_impl.F90 | 28 +-- base/serial/impl/psb_z_coo_impl.f90 | 29 ++- base/serial/impl/psb_z_mat_impl.F90 | 28 +-- base/tools/psb_csphalo.F90 | 358 ++++++++++++++-------------- base/tools/psb_dsphalo.F90 | 358 ++++++++++++++-------------- base/tools/psb_ssphalo.F90 | 358 ++++++++++++++-------------- base/tools/psb_zsphalo.F90 | 358 ++++++++++++++-------------- 12 files changed, 796 insertions(+), 864 deletions(-) diff --git a/base/serial/impl/psb_c_coo_impl.f90 b/base/serial/impl/psb_c_coo_impl.f90 index 9976ed77..142f5b61 100644 --- a/base/serial/impl/psb_c_coo_impl.f90 +++ b/base/serial/impl/psb_c_coo_impl.f90 @@ -2164,6 +2164,10 @@ contains end subroutine psb_c_coo_csgetptn +! +! NZ is the number of non-zeros on output. +! The output is guaranteed to be sorted +! subroutine psb_c_coo_csgetrow(imin,imax,a,nz,ia,ja,val,info,& & jmin,jmax,iren,append,nzin,rscale,cscale) ! Output is always in COO format @@ -2271,6 +2275,7 @@ contains use psb_error_mod use psb_realloc_mod use psb_sort_mod + use psb_ip_reord_mod implicit none class(psb_c_coo_sparse_mat), intent(in) :: a @@ -2369,25 +2374,25 @@ contains if (present(iren)) then do i=ip,jp if ((jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then - nzin_ = nzin_ + 1 nz = nz + 1 - val(nzin_) = a%val(i) - ia(nzin_) = iren(a%ia(i)) - ja(nzin_) = iren(a%ja(i)) + val(nzin_+nz) = a%val(i) + ia(nzin_+nz) = iren(a%ia(i)) + ja(nzin_+nz) = iren(a%ja(i)) end if enddo + call psb_c_fix_coo_inner(nzin_+nz,psb_dupl_add_,ia,ja,val,nz,info) + nz = nz - nzin_ else do i=ip,jp if ((jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then - nzin_ = nzin_ + 1 nz = nz + 1 - val(nzin_) = a%val(i) - ia(nzin_) = a%ia(i) - ja(nzin_) = a%ja(i) + val(nzin_+nz) = a%val(i) + ia(nzin_+nz) = a%ia(i) + ja(nzin_+nz) = a%ja(i) end if enddo end if - else + else nz = 0 end if @@ -2438,9 +2443,9 @@ contains ja(nzin_+k) = (a%ja(i)) endif enddo - nzin_=nzin_+k - end if - nz = k + end if + call psb_c_fix_coo_inner(nzin_+k,psb_dupl_add_,ia,ja,val,nz,info) + nz = nz - nzin_ end if end subroutine coo_getrow diff --git a/base/serial/impl/psb_c_mat_impl.F90 b/base/serial/impl/psb_c_mat_impl.F90 index 2078fb82..9bfb1008 100644 --- a/base/serial/impl/psb_c_mat_impl.F90 +++ b/base/serial/impl/psb_c_mat_impl.F90 @@ -867,7 +867,7 @@ subroutine psb_c_csgetblk(imin,imax,a,b,info,& Integer :: err_act character(len=20) :: name='csget' logical, parameter :: debug=.false. - class(psb_c_base_sparse_mat), allocatable :: acoo + type(psb_c_coo_sparse_mat), allocatable :: acoo info = psb_success_ @@ -878,17 +878,11 @@ subroutine psb_c_csgetblk(imin,imax,a,b,info,& goto 9999 endif - allocate(psb_c_coo_sparse_mat :: acoo,stat=info) + allocate(acoo,stat=info) if (info == psb_success_) then - select type (acoo) - type is (psb_c_coo_sparse_mat) - call a%a%csget(imin,imax,acoo,info,& - & jmin,jmax,iren,append,rscale,cscale) - class default - ! This is impossible - info = psb_err_internal_error_ - end select + call a%a%csget(imin,imax,acoo,info,& + & jmin,jmax,iren,append,rscale,cscale) else info = psb_err_alloc_dealloc_ end if @@ -929,7 +923,7 @@ subroutine psb_c_csclip(a,b,info,& Integer :: err_act character(len=20) :: name='csclip' logical, parameter :: debug=.false. - class(psb_c_base_sparse_mat), allocatable :: acoo + type(psb_c_coo_sparse_mat), allocatable :: acoo info = psb_success_ call psb_erractionsave(err_act) @@ -939,17 +933,11 @@ subroutine psb_c_csclip(a,b,info,& goto 9999 endif - allocate(psb_c_coo_sparse_mat :: acoo,stat=info) + allocate(acoo,stat=info) if (info == psb_success_) then - select type (acoo) - type is (psb_c_coo_sparse_mat) - call a%a%csclip(acoo,info,& - & imin,imax,jmin,jmax,rscale,cscale) - class default - ! This is impossible - info = psb_err_internal_error_ - end select + call a%a%csclip(acoo,info,& + & imin,imax,jmin,jmax,rscale,cscale) else info = psb_err_alloc_dealloc_ end if diff --git a/base/serial/impl/psb_d_coo_impl.f90 b/base/serial/impl/psb_d_coo_impl.f90 index a5dd2b6b..fa7b0e8d 100644 --- a/base/serial/impl/psb_d_coo_impl.f90 +++ b/base/serial/impl/psb_d_coo_impl.f90 @@ -2164,6 +2164,10 @@ contains end subroutine psb_d_coo_csgetptn +! +! NZ is the number of non-zeros on output. +! The output is guaranteed to be sorted +! subroutine psb_d_coo_csgetrow(imin,imax,a,nz,ia,ja,val,info,& & jmin,jmax,iren,append,nzin,rscale,cscale) ! Output is always in COO format @@ -2271,6 +2275,7 @@ contains use psb_error_mod use psb_realloc_mod use psb_sort_mod + use psb_ip_reord_mod implicit none class(psb_d_coo_sparse_mat), intent(in) :: a @@ -2369,25 +2374,25 @@ contains if (present(iren)) then do i=ip,jp if ((jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then - nzin_ = nzin_ + 1 nz = nz + 1 - val(nzin_) = a%val(i) - ia(nzin_) = iren(a%ia(i)) - ja(nzin_) = iren(a%ja(i)) + val(nzin_+nz) = a%val(i) + ia(nzin_+nz) = iren(a%ia(i)) + ja(nzin_+nz) = iren(a%ja(i)) end if enddo + call psb_d_fix_coo_inner(nzin_+nz,psb_dupl_add_,ia,ja,val,nz,info) + nz = nz - nzin_ else do i=ip,jp if ((jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then - nzin_ = nzin_ + 1 nz = nz + 1 - val(nzin_) = a%val(i) - ia(nzin_) = a%ia(i) - ja(nzin_) = a%ja(i) + val(nzin_+nz) = a%val(i) + ia(nzin_+nz) = a%ia(i) + ja(nzin_+nz) = a%ja(i) end if enddo end if - else + else nz = 0 end if @@ -2438,9 +2443,9 @@ contains ja(nzin_+k) = (a%ja(i)) endif enddo - nzin_=nzin_+k - end if - nz = k + end if + call psb_d_fix_coo_inner(nzin_+k,psb_dupl_add_,ia,ja,val,nz,info) + nz = nz - nzin_ end if end subroutine coo_getrow diff --git a/base/serial/impl/psb_d_mat_impl.F90 b/base/serial/impl/psb_d_mat_impl.F90 index bc45e56b..fc7a5b61 100644 --- a/base/serial/impl/psb_d_mat_impl.F90 +++ b/base/serial/impl/psb_d_mat_impl.F90 @@ -867,7 +867,7 @@ subroutine psb_d_csgetblk(imin,imax,a,b,info,& Integer :: err_act character(len=20) :: name='csget' logical, parameter :: debug=.false. - class(psb_d_base_sparse_mat), allocatable :: acoo + type(psb_d_coo_sparse_mat), allocatable :: acoo info = psb_success_ @@ -878,17 +878,11 @@ subroutine psb_d_csgetblk(imin,imax,a,b,info,& goto 9999 endif - allocate(psb_d_coo_sparse_mat :: acoo,stat=info) + allocate(acoo,stat=info) if (info == psb_success_) then - select type (acoo) - type is (psb_d_coo_sparse_mat) - call a%a%csget(imin,imax,acoo,info,& - & jmin,jmax,iren,append,rscale,cscale) - class default - ! This is impossible - info = psb_err_internal_error_ - end select + call a%a%csget(imin,imax,acoo,info,& + & jmin,jmax,iren,append,rscale,cscale) else info = psb_err_alloc_dealloc_ end if @@ -929,7 +923,7 @@ subroutine psb_d_csclip(a,b,info,& Integer :: err_act character(len=20) :: name='csclip' logical, parameter :: debug=.false. - class(psb_d_base_sparse_mat), allocatable :: acoo + type(psb_d_coo_sparse_mat), allocatable :: acoo info = psb_success_ call psb_erractionsave(err_act) @@ -939,17 +933,11 @@ subroutine psb_d_csclip(a,b,info,& goto 9999 endif - allocate(psb_d_coo_sparse_mat :: acoo,stat=info) + allocate(acoo,stat=info) if (info == psb_success_) then - select type (acoo) - type is (psb_d_coo_sparse_mat) - call a%a%csclip(acoo,info,& - & imin,imax,jmin,jmax,rscale,cscale) - class default - ! This is impossible - info = psb_err_internal_error_ - end select + call a%a%csclip(acoo,info,& + & imin,imax,jmin,jmax,rscale,cscale) else info = psb_err_alloc_dealloc_ end if diff --git a/base/serial/impl/psb_s_coo_impl.f90 b/base/serial/impl/psb_s_coo_impl.f90 index a89fd662..a2e6403b 100644 --- a/base/serial/impl/psb_s_coo_impl.f90 +++ b/base/serial/impl/psb_s_coo_impl.f90 @@ -2164,6 +2164,10 @@ contains end subroutine psb_s_coo_csgetptn +! +! NZ is the number of non-zeros on output. +! The output is guaranteed to be sorted +! subroutine psb_s_coo_csgetrow(imin,imax,a,nz,ia,ja,val,info,& & jmin,jmax,iren,append,nzin,rscale,cscale) ! Output is always in COO format @@ -2271,6 +2275,7 @@ contains use psb_error_mod use psb_realloc_mod use psb_sort_mod + use psb_ip_reord_mod implicit none class(psb_s_coo_sparse_mat), intent(in) :: a @@ -2369,25 +2374,25 @@ contains if (present(iren)) then do i=ip,jp if ((jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then - nzin_ = nzin_ + 1 nz = nz + 1 - val(nzin_) = a%val(i) - ia(nzin_) = iren(a%ia(i)) - ja(nzin_) = iren(a%ja(i)) + val(nzin_+nz) = a%val(i) + ia(nzin_+nz) = iren(a%ia(i)) + ja(nzin_+nz) = iren(a%ja(i)) end if enddo + call psb_s_fix_coo_inner(nzin_+nz,psb_dupl_add_,ia,ja,val,nz,info) + nz = nz - nzin_ else do i=ip,jp if ((jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then - nzin_ = nzin_ + 1 nz = nz + 1 - val(nzin_) = a%val(i) - ia(nzin_) = a%ia(i) - ja(nzin_) = a%ja(i) + val(nzin_+nz) = a%val(i) + ia(nzin_+nz) = a%ia(i) + ja(nzin_+nz) = a%ja(i) end if enddo end if - else + else nz = 0 end if @@ -2438,9 +2443,9 @@ contains ja(nzin_+k) = (a%ja(i)) endif enddo - nzin_=nzin_+k - end if - nz = k + end if + call psb_s_fix_coo_inner(nzin_+k,psb_dupl_add_,ia,ja,val,nz,info) + nz = nz - nzin_ end if end subroutine coo_getrow diff --git a/base/serial/impl/psb_s_mat_impl.F90 b/base/serial/impl/psb_s_mat_impl.F90 index de6e5679..bfd59727 100644 --- a/base/serial/impl/psb_s_mat_impl.F90 +++ b/base/serial/impl/psb_s_mat_impl.F90 @@ -867,7 +867,7 @@ subroutine psb_s_csgetblk(imin,imax,a,b,info,& Integer :: err_act character(len=20) :: name='csget' logical, parameter :: debug=.false. - class(psb_s_base_sparse_mat), allocatable :: acoo + type(psb_s_coo_sparse_mat), allocatable :: acoo info = psb_success_ @@ -878,17 +878,11 @@ subroutine psb_s_csgetblk(imin,imax,a,b,info,& goto 9999 endif - allocate(psb_s_coo_sparse_mat :: acoo,stat=info) + allocate(acoo,stat=info) if (info == psb_success_) then - select type (acoo) - type is (psb_s_coo_sparse_mat) - call a%a%csget(imin,imax,acoo,info,& - & jmin,jmax,iren,append,rscale,cscale) - class default - ! This is impossible - info = psb_err_internal_error_ - end select + call a%a%csget(imin,imax,acoo,info,& + & jmin,jmax,iren,append,rscale,cscale) else info = psb_err_alloc_dealloc_ end if @@ -929,7 +923,7 @@ subroutine psb_s_csclip(a,b,info,& Integer :: err_act character(len=20) :: name='csclip' logical, parameter :: debug=.false. - class(psb_s_base_sparse_mat), allocatable :: acoo + type(psb_s_coo_sparse_mat), allocatable :: acoo info = psb_success_ call psb_erractionsave(err_act) @@ -939,17 +933,11 @@ subroutine psb_s_csclip(a,b,info,& goto 9999 endif - allocate(psb_s_coo_sparse_mat :: acoo,stat=info) + allocate(acoo,stat=info) if (info == psb_success_) then - select type (acoo) - type is (psb_s_coo_sparse_mat) - call a%a%csclip(acoo,info,& - & imin,imax,jmin,jmax,rscale,cscale) - class default - ! This is impossible - info = psb_err_internal_error_ - end select + call a%a%csclip(acoo,info,& + & imin,imax,jmin,jmax,rscale,cscale) else info = psb_err_alloc_dealloc_ end if diff --git a/base/serial/impl/psb_z_coo_impl.f90 b/base/serial/impl/psb_z_coo_impl.f90 index e4a7a6eb..542410a1 100644 --- a/base/serial/impl/psb_z_coo_impl.f90 +++ b/base/serial/impl/psb_z_coo_impl.f90 @@ -2164,6 +2164,10 @@ contains end subroutine psb_z_coo_csgetptn +! +! NZ is the number of non-zeros on output. +! The output is guaranteed to be sorted +! subroutine psb_z_coo_csgetrow(imin,imax,a,nz,ia,ja,val,info,& & jmin,jmax,iren,append,nzin,rscale,cscale) ! Output is always in COO format @@ -2271,6 +2275,7 @@ contains use psb_error_mod use psb_realloc_mod use psb_sort_mod + use psb_ip_reord_mod implicit none class(psb_z_coo_sparse_mat), intent(in) :: a @@ -2369,25 +2374,25 @@ contains if (present(iren)) then do i=ip,jp if ((jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then - nzin_ = nzin_ + 1 nz = nz + 1 - val(nzin_) = a%val(i) - ia(nzin_) = iren(a%ia(i)) - ja(nzin_) = iren(a%ja(i)) + val(nzin_+nz) = a%val(i) + ia(nzin_+nz) = iren(a%ia(i)) + ja(nzin_+nz) = iren(a%ja(i)) end if enddo + call psb_z_fix_coo_inner(nzin_+nz,psb_dupl_add_,ia,ja,val,nz,info) + nz = nz - nzin_ else do i=ip,jp if ((jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then - nzin_ = nzin_ + 1 nz = nz + 1 - val(nzin_) = a%val(i) - ia(nzin_) = a%ia(i) - ja(nzin_) = a%ja(i) + val(nzin_+nz) = a%val(i) + ia(nzin_+nz) = a%ia(i) + ja(nzin_+nz) = a%ja(i) end if enddo end if - else + else nz = 0 end if @@ -2438,9 +2443,9 @@ contains ja(nzin_+k) = (a%ja(i)) endif enddo - nzin_=nzin_+k - end if - nz = k + end if + call psb_z_fix_coo_inner(nzin_+k,psb_dupl_add_,ia,ja,val,nz,info) + nz = nz - nzin_ end if end subroutine coo_getrow diff --git a/base/serial/impl/psb_z_mat_impl.F90 b/base/serial/impl/psb_z_mat_impl.F90 index 4fd791cd..c8423502 100644 --- a/base/serial/impl/psb_z_mat_impl.F90 +++ b/base/serial/impl/psb_z_mat_impl.F90 @@ -867,7 +867,7 @@ subroutine psb_z_csgetblk(imin,imax,a,b,info,& Integer :: err_act character(len=20) :: name='csget' logical, parameter :: debug=.false. - class(psb_z_base_sparse_mat), allocatable :: acoo + type(psb_z_coo_sparse_mat), allocatable :: acoo info = psb_success_ @@ -878,17 +878,11 @@ subroutine psb_z_csgetblk(imin,imax,a,b,info,& goto 9999 endif - allocate(psb_z_coo_sparse_mat :: acoo,stat=info) + allocate(acoo,stat=info) if (info == psb_success_) then - select type (acoo) - type is (psb_z_coo_sparse_mat) - call a%a%csget(imin,imax,acoo,info,& - & jmin,jmax,iren,append,rscale,cscale) - class default - ! This is impossible - info = psb_err_internal_error_ - end select + call a%a%csget(imin,imax,acoo,info,& + & jmin,jmax,iren,append,rscale,cscale) else info = psb_err_alloc_dealloc_ end if @@ -929,7 +923,7 @@ subroutine psb_z_csclip(a,b,info,& Integer :: err_act character(len=20) :: name='csclip' logical, parameter :: debug=.false. - class(psb_z_base_sparse_mat), allocatable :: acoo + type(psb_z_coo_sparse_mat), allocatable :: acoo info = psb_success_ call psb_erractionsave(err_act) @@ -939,17 +933,11 @@ subroutine psb_z_csclip(a,b,info,& goto 9999 endif - allocate(psb_z_coo_sparse_mat :: acoo,stat=info) + allocate(acoo,stat=info) if (info == psb_success_) then - select type (acoo) - type is (psb_z_coo_sparse_mat) - call a%a%csclip(acoo,info,& - & imin,imax,jmin,jmax,rscale,cscale) - class default - ! This is impossible - info = psb_err_internal_error_ - end select + call a%a%csclip(acoo,info,& + & imin,imax,jmin,jmax,rscale,cscale) else info = psb_err_alloc_dealloc_ end if diff --git a/base/tools/psb_csphalo.F90 b/base/tools/psb_csphalo.F90 index 1335da95..7471cc64 100644 --- a/base/tools/psb_csphalo.F90 +++ b/base/tools/psb_csphalo.F90 @@ -83,7 +83,7 @@ Subroutine psb_csphalo(a,desc_a,blk,info,rowcnv,colcnv,& Integer, allocatable :: sdid(:,:), brvindx(:),rvid(:,:), & & rvsz(:), bsdindx(:),sdsz(:), iasnd(:), jasnd(:) complex(psb_spk_), allocatable :: valsnd(:) - class(psb_c_base_sparse_mat), allocatable :: acoo + type(psb_c_coo_sparse_mat), allocatable :: acoo integer, pointer :: idxv(:) logical :: rowcnv_,colcnv_,rowscale_,colscale_ character(len=5) :: outfmt_ @@ -164,200 +164,190 @@ Subroutine psb_csphalo(a,desc_a,blk,info,rowcnv,colcnv,& end select - allocate(psb_c_coo_sparse_mat :: acoo,stat=info) + allocate(acoo,stat=info) if (info /= psb_success_) then info = psb_err_alloc_dealloc_ call psb_errpush(info,name) goto 9999 end if - select type (acoo) - type is (psb_c_coo_sparse_mat) - - - l1 = 0 - - sdsz(:)=0 - rvsz(:)=0 - ipx = 1 - brvindx(ipx) = 0 - bsdindx(ipx) = 0 - counter=1 - idx = 0 - idxs = 0 - idxr = 0 - call acoo%allocate(0,a%get_ncols(),info) - ! For all rows in the halo descriptor, extract and send/receive. - Do - proc=idxv(counter) - if (proc == -1) exit - n_el_recv = idxv(counter+psb_n_elem_recv_) - counter = counter+n_el_recv - n_el_send = idxv(counter+psb_n_elem_send_) - tot_elem = 0 - Do j=0,n_el_send-1 - idx = idxv(counter+psb_elem_send_+j) - n_elem = a%get_nz_row(idx) - tot_elem = tot_elem+n_elem - Enddo - sdsz(proc+1) = tot_elem - call acoo%set_nrows(acoo%get_nrows() + n_el_recv) - counter = counter+n_el_send+3 - Enddo - call mpi_alltoall(sdsz,1,mpi_integer,rvsz,1,mpi_integer,icomm,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='mpi_alltoall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - idxs = 0 - idxr = 0 - counter = 1 - Do - proc=idxv(counter) - if (proc == -1) exit - n_el_recv = idxv(counter+psb_n_elem_recv_) - counter = counter+n_el_recv - n_el_send = idxv(counter+psb_n_elem_send_) - - bsdindx(proc+1) = idxs - idxs = idxs + sdsz(proc+1) - brvindx(proc+1) = idxr - idxr = idxr + rvsz(proc+1) - counter = counter+n_el_send+3 + l1 = 0 + + sdsz(:)=0 + rvsz(:)=0 + ipx = 1 + brvindx(ipx) = 0 + bsdindx(ipx) = 0 + counter=1 + idx = 0 + idxs = 0 + idxr = 0 + call acoo%allocate(0,a%get_ncols(),info) + ! For all rows in the halo descriptor, extract and send/receive. + Do + proc=idxv(counter) + if (proc == -1) exit + n_el_recv = idxv(counter+psb_n_elem_recv_) + counter = counter+n_el_recv + n_el_send = idxv(counter+psb_n_elem_send_) + tot_elem = 0 + Do j=0,n_el_send-1 + idx = idxv(counter+psb_elem_send_+j) + n_elem = a%get_nz_row(idx) + tot_elem = tot_elem+n_elem Enddo + sdsz(proc+1) = tot_elem + call acoo%set_nrows(acoo%get_nrows() + n_el_recv) + counter = counter+n_el_send+3 + Enddo - iszr=sum(rvsz) - call acoo%reallocate(max(iszr,1)) - if (debug_level >= psb_debug_outer_)& - & write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),& - & ' Send:',sdsz(:),' Receive:',rvsz(:) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_sp_reall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - mat_recv = iszr - iszs=sum(sdsz) - call psb_ensure_size(max(iszs,1),iasnd,info) - if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info) - if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info) - - l1 = 0 - ipx = 1 - counter=1 - idx = 0 - - tot_elem=0 - Do - proc=idxv(counter) - if (proc == -1) exit - n_el_recv=idxv(counter+psb_n_elem_recv_) - counter=counter+n_el_recv - n_el_send=idxv(counter+psb_n_elem_send_) - - Do j=0,n_el_send-1 - idx = idxv(counter+psb_elem_send_+j) - n_elem = a%get_nz_row(idx) - call a%csget(idx,idx,ngtz,iasnd,jasnd,valsnd,info,& - & append=.true.,nzin=tot_elem) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_sp_getrow' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - tot_elem=tot_elem+n_elem - Enddo - ipx = ipx + 1 - counter = counter+n_el_send+3 - Enddo - nz = tot_elem - - 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 /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_loc_to_glob' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - - call mpi_alltoallv(valsnd,sdsz,bsdindx,mpi_complex,& - & acoo%val,rvsz,brvindx,mpi_complex,icomm,info) - call mpi_alltoallv(iasnd,sdsz,bsdindx,mpi_integer,& - & acoo%ia,rvsz,brvindx,mpi_integer,icomm,info) - call mpi_alltoallv(jasnd,sdsz,bsdindx,mpi_integer,& - & acoo%ja,rvsz,brvindx,mpi_integer,icomm,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='mpi_alltoallv' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - ! - ! Convert into local numbering - ! - if (rowcnv_) call psb_glob_to_loc(acoo%ia(1:iszr),desc_a,info,iact='I') - if (colcnv_) call psb_glob_to_loc(acoo%ja(1:iszr),desc_a,info,iact='I') - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psbglob_to_loc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - l1 = 0 - call acoo%set_nrows(0) - ! - irmin = huge(irmin) - icmin = huge(icmin) - irmax = 0 - icmax = 0 - Do i=1,iszr - r=(acoo%ia(i)) - k=(acoo%ja(i)) - ! Just in case some of the conversions were out-of-range - If ((r>0).and.(k>0)) Then - l1=l1+1 - acoo%val(l1) = acoo%val(i) - acoo%ia(l1) = r - acoo%ja(l1) = k - irmin = min(irmin,r) - irmax = max(irmax,r) - icmin = min(icmin,k) - icmax = max(icmax,k) - End If + call mpi_alltoall(sdsz,1,mpi_integer,rvsz,1,mpi_integer,icomm,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='mpi_alltoall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + idxs = 0 + idxr = 0 + counter = 1 + Do + proc=idxv(counter) + if (proc == -1) exit + n_el_recv = idxv(counter+psb_n_elem_recv_) + counter = counter+n_el_recv + n_el_send = idxv(counter+psb_n_elem_send_) + + bsdindx(proc+1) = idxs + idxs = idxs + sdsz(proc+1) + brvindx(proc+1) = idxr + idxr = idxr + rvsz(proc+1) + counter = counter+n_el_send+3 + Enddo + + iszr=sum(rvsz) + call acoo%reallocate(max(iszr,1)) + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),& + & ' Send:',sdsz(:),' Receive:',rvsz(:) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_reall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + mat_recv = iszr + iszs=sum(sdsz) + call psb_ensure_size(max(iszs,1),iasnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info) + + l1 = 0 + ipx = 1 + counter=1 + idx = 0 + + tot_elem=0 + Do + proc=idxv(counter) + if (proc == -1) exit + n_el_recv=idxv(counter+psb_n_elem_recv_) + counter=counter+n_el_recv + n_el_send=idxv(counter+psb_n_elem_send_) + + Do j=0,n_el_send-1 + idx = idxv(counter+psb_elem_send_+j) + n_elem = a%get_nz_row(idx) + call a%csget(idx,idx,ngtz,iasnd,jasnd,valsnd,info,& + & append=.true.,nzin=tot_elem) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_getrow' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + tot_elem=tot_elem+n_elem Enddo - if (rowscale_) then - call acoo%set_nrows(max(irmax-irmin+1,0)) - acoo%ia(1:l1) = acoo%ia(1:l1) - irmin + 1 - else - call acoo%set_nrows(irmax) - end if - if (colscale_) then - call acoo%set_ncols(max(icmax-icmin+1,0)) - acoo%ja(1:l1) = acoo%ja(1:l1) - icmin + 1 - else - call acoo%set_ncols(icmax) - end if - - call acoo%set_nzeros(l1) - - - class default - ! This is impossible - info = psb_err_internal_error_ - call psb_Errpush(info,name) + ipx = ipx + 1 + counter = counter+n_el_send+3 + Enddo + nz = tot_elem + + 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 /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_loc_to_glob' + call psb_errpush(info,name,a_err=ch_err) goto 9999 - end select + end if + + + call mpi_alltoallv(valsnd,sdsz,bsdindx,mpi_complex,& + & acoo%val,rvsz,brvindx,mpi_complex,icomm,info) + call mpi_alltoallv(iasnd,sdsz,bsdindx,mpi_integer,& + & acoo%ia,rvsz,brvindx,mpi_integer,icomm,info) + call mpi_alltoallv(jasnd,sdsz,bsdindx,mpi_integer,& + & acoo%ja,rvsz,brvindx,mpi_integer,icomm,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='mpi_alltoallv' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! + ! Convert into local numbering + ! + if (rowcnv_) call psb_glob_to_loc(acoo%ia(1:iszr),desc_a,info,iact='I') + if (colcnv_) call psb_glob_to_loc(acoo%ja(1:iszr),desc_a,info,iact='I') + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psbglob_to_loc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + l1 = 0 + call acoo%set_nrows(0) + ! + irmin = huge(irmin) + icmin = huge(icmin) + irmax = 0 + icmax = 0 + Do i=1,iszr + r=(acoo%ia(i)) + k=(acoo%ja(i)) + ! Just in case some of the conversions were out-of-range + If ((r>0).and.(k>0)) Then + l1=l1+1 + acoo%val(l1) = acoo%val(i) + acoo%ia(l1) = r + acoo%ja(l1) = k + irmin = min(irmin,r) + irmax = max(irmax,r) + icmin = min(icmin,k) + icmax = max(icmax,k) + End If + Enddo + if (rowscale_) then + call acoo%set_nrows(max(irmax-irmin+1,0)) + acoo%ia(1:l1) = acoo%ia(1:l1) - irmin + 1 + else + call acoo%set_nrows(irmax) + end if + if (colscale_) then + call acoo%set_ncols(max(icmax-icmin+1,0)) + acoo%ja(1:l1) = acoo%ja(1:l1) - icmin + 1 + else + call acoo%set_ncols(icmax) + end if + + call acoo%set_nzeros(l1) + if (debug_level >= psb_debug_outer_)& & write(debug_unit,*) me,' ',trim(name),& diff --git a/base/tools/psb_dsphalo.F90 b/base/tools/psb_dsphalo.F90 index 2896649d..568e4655 100644 --- a/base/tools/psb_dsphalo.F90 +++ b/base/tools/psb_dsphalo.F90 @@ -83,7 +83,7 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& Integer, allocatable :: sdid(:,:), brvindx(:),rvid(:,:), & & rvsz(:), bsdindx(:),sdsz(:), iasnd(:), jasnd(:) real(psb_dpk_), allocatable :: valsnd(:) - class(psb_d_base_sparse_mat), allocatable :: acoo + type(psb_d_coo_sparse_mat), allocatable :: acoo integer, pointer :: idxv(:) logical :: rowcnv_,colcnv_,rowscale_,colscale_ character(len=5) :: outfmt_ @@ -164,200 +164,190 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& end select - allocate(psb_d_coo_sparse_mat :: acoo,stat=info) + allocate(acoo,stat=info) if (info /= psb_success_) then info = psb_err_alloc_dealloc_ call psb_errpush(info,name) goto 9999 end if - select type (acoo) - type is (psb_d_coo_sparse_mat) - - - l1 = 0 - - sdsz(:)=0 - rvsz(:)=0 - ipx = 1 - brvindx(ipx) = 0 - bsdindx(ipx) = 0 - counter=1 - idx = 0 - idxs = 0 - idxr = 0 - call acoo%allocate(0,a%get_ncols(),info) - ! For all rows in the halo descriptor, extract and send/receive. - Do - proc=idxv(counter) - if (proc == -1) exit - n_el_recv = idxv(counter+psb_n_elem_recv_) - counter = counter+n_el_recv - n_el_send = idxv(counter+psb_n_elem_send_) - tot_elem = 0 - Do j=0,n_el_send-1 - idx = idxv(counter+psb_elem_send_+j) - n_elem = a%get_nz_row(idx) - tot_elem = tot_elem+n_elem - Enddo - sdsz(proc+1) = tot_elem - call acoo%set_nrows(acoo%get_nrows() + n_el_recv) - counter = counter+n_el_send+3 - Enddo - call mpi_alltoall(sdsz,1,mpi_integer,rvsz,1,mpi_integer,icomm,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='mpi_alltoall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - idxs = 0 - idxr = 0 - counter = 1 - Do - proc=idxv(counter) - if (proc == -1) exit - n_el_recv = idxv(counter+psb_n_elem_recv_) - counter = counter+n_el_recv - n_el_send = idxv(counter+psb_n_elem_send_) - - bsdindx(proc+1) = idxs - idxs = idxs + sdsz(proc+1) - brvindx(proc+1) = idxr - idxr = idxr + rvsz(proc+1) - counter = counter+n_el_send+3 + l1 = 0 + + sdsz(:)=0 + rvsz(:)=0 + ipx = 1 + brvindx(ipx) = 0 + bsdindx(ipx) = 0 + counter=1 + idx = 0 + idxs = 0 + idxr = 0 + call acoo%allocate(0,a%get_ncols(),info) + ! For all rows in the halo descriptor, extract and send/receive. + Do + proc=idxv(counter) + if (proc == -1) exit + n_el_recv = idxv(counter+psb_n_elem_recv_) + counter = counter+n_el_recv + n_el_send = idxv(counter+psb_n_elem_send_) + tot_elem = 0 + Do j=0,n_el_send-1 + idx = idxv(counter+psb_elem_send_+j) + n_elem = a%get_nz_row(idx) + tot_elem = tot_elem+n_elem Enddo + sdsz(proc+1) = tot_elem + call acoo%set_nrows(acoo%get_nrows() + n_el_recv) + counter = counter+n_el_send+3 + Enddo - iszr=sum(rvsz) - call acoo%reallocate(max(iszr,1)) - if (debug_level >= psb_debug_outer_)& - & write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),& - & ' Send:',sdsz(:),' Receive:',rvsz(:) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_sp_reall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - mat_recv = iszr - iszs=sum(sdsz) - call psb_ensure_size(max(iszs,1),iasnd,info) - if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info) - if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info) - - l1 = 0 - ipx = 1 - counter=1 - idx = 0 - - tot_elem=0 - Do - proc=idxv(counter) - if (proc == -1) exit - n_el_recv=idxv(counter+psb_n_elem_recv_) - counter=counter+n_el_recv - n_el_send=idxv(counter+psb_n_elem_send_) - - Do j=0,n_el_send-1 - idx = idxv(counter+psb_elem_send_+j) - n_elem = a%get_nz_row(idx) - call a%csget(idx,idx,ngtz,iasnd,jasnd,valsnd,info,& - & append=.true.,nzin=tot_elem) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_sp_getrow' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - tot_elem=tot_elem+n_elem - Enddo - ipx = ipx + 1 - counter = counter+n_el_send+3 - Enddo - nz = tot_elem - - 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 /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_loc_to_glob' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - - call mpi_alltoallv(valsnd,sdsz,bsdindx,mpi_double_precision,& - & acoo%val,rvsz,brvindx,mpi_double_precision,icomm,info) - call mpi_alltoallv(iasnd,sdsz,bsdindx,mpi_integer,& - & acoo%ia,rvsz,brvindx,mpi_integer,icomm,info) - call mpi_alltoallv(jasnd,sdsz,bsdindx,mpi_integer,& - & acoo%ja,rvsz,brvindx,mpi_integer,icomm,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='mpi_alltoallv' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - ! - ! Convert into local numbering - ! - if (rowcnv_) call psb_glob_to_loc(acoo%ia(1:iszr),desc_a,info,iact='I') - if (colcnv_) call psb_glob_to_loc(acoo%ja(1:iszr),desc_a,info,iact='I') - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psbglob_to_loc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - l1 = 0 - call acoo%set_nrows(0) - ! - irmin = huge(irmin) - icmin = huge(icmin) - irmax = 0 - icmax = 0 - Do i=1,iszr - r=(acoo%ia(i)) - k=(acoo%ja(i)) - ! Just in case some of the conversions were out-of-range - If ((r>0).and.(k>0)) Then - l1=l1+1 - acoo%val(l1) = acoo%val(i) - acoo%ia(l1) = r - acoo%ja(l1) = k - irmin = min(irmin,r) - irmax = max(irmax,r) - icmin = min(icmin,k) - icmax = max(icmax,k) - End If + call mpi_alltoall(sdsz,1,mpi_integer,rvsz,1,mpi_integer,icomm,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='mpi_alltoall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + idxs = 0 + idxr = 0 + counter = 1 + Do + proc=idxv(counter) + if (proc == -1) exit + n_el_recv = idxv(counter+psb_n_elem_recv_) + counter = counter+n_el_recv + n_el_send = idxv(counter+psb_n_elem_send_) + + bsdindx(proc+1) = idxs + idxs = idxs + sdsz(proc+1) + brvindx(proc+1) = idxr + idxr = idxr + rvsz(proc+1) + counter = counter+n_el_send+3 + Enddo + + iszr=sum(rvsz) + call acoo%reallocate(max(iszr,1)) + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),& + & ' Send:',sdsz(:),' Receive:',rvsz(:) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_reall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + mat_recv = iszr + iszs=sum(sdsz) + call psb_ensure_size(max(iszs,1),iasnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info) + + l1 = 0 + ipx = 1 + counter=1 + idx = 0 + + tot_elem=0 + Do + proc=idxv(counter) + if (proc == -1) exit + n_el_recv=idxv(counter+psb_n_elem_recv_) + counter=counter+n_el_recv + n_el_send=idxv(counter+psb_n_elem_send_) + + Do j=0,n_el_send-1 + idx = idxv(counter+psb_elem_send_+j) + n_elem = a%get_nz_row(idx) + call a%csget(idx,idx,ngtz,iasnd,jasnd,valsnd,info,& + & append=.true.,nzin=tot_elem) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_getrow' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + tot_elem=tot_elem+n_elem Enddo - if (rowscale_) then - call acoo%set_nrows(max(irmax-irmin+1,0)) - acoo%ia(1:l1) = acoo%ia(1:l1) - irmin + 1 - else - call acoo%set_nrows(irmax) - end if - if (colscale_) then - call acoo%set_ncols(max(icmax-icmin+1,0)) - acoo%ja(1:l1) = acoo%ja(1:l1) - icmin + 1 - else - call acoo%set_ncols(icmax) - end if - - call acoo%set_nzeros(l1) - - - class default - ! This is impossible - info = psb_err_internal_error_ - call psb_Errpush(info,name) + ipx = ipx + 1 + counter = counter+n_el_send+3 + Enddo + nz = tot_elem + + 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 /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_loc_to_glob' + call psb_errpush(info,name,a_err=ch_err) goto 9999 - end select + end if + + + call mpi_alltoallv(valsnd,sdsz,bsdindx,mpi_double_precision,& + & acoo%val,rvsz,brvindx,mpi_double_precision,icomm,info) + call mpi_alltoallv(iasnd,sdsz,bsdindx,mpi_integer,& + & acoo%ia,rvsz,brvindx,mpi_integer,icomm,info) + call mpi_alltoallv(jasnd,sdsz,bsdindx,mpi_integer,& + & acoo%ja,rvsz,brvindx,mpi_integer,icomm,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='mpi_alltoallv' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! + ! Convert into local numbering + ! + if (rowcnv_) call psb_glob_to_loc(acoo%ia(1:iszr),desc_a,info,iact='I') + if (colcnv_) call psb_glob_to_loc(acoo%ja(1:iszr),desc_a,info,iact='I') + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psbglob_to_loc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + l1 = 0 + call acoo%set_nrows(0) + ! + irmin = huge(irmin) + icmin = huge(icmin) + irmax = 0 + icmax = 0 + Do i=1,iszr + r=(acoo%ia(i)) + k=(acoo%ja(i)) + ! Just in case some of the conversions were out-of-range + If ((r>0).and.(k>0)) Then + l1=l1+1 + acoo%val(l1) = acoo%val(i) + acoo%ia(l1) = r + acoo%ja(l1) = k + irmin = min(irmin,r) + irmax = max(irmax,r) + icmin = min(icmin,k) + icmax = max(icmax,k) + End If + Enddo + if (rowscale_) then + call acoo%set_nrows(max(irmax-irmin+1,0)) + acoo%ia(1:l1) = acoo%ia(1:l1) - irmin + 1 + else + call acoo%set_nrows(irmax) + end if + if (colscale_) then + call acoo%set_ncols(max(icmax-icmin+1,0)) + acoo%ja(1:l1) = acoo%ja(1:l1) - icmin + 1 + else + call acoo%set_ncols(icmax) + end if + + call acoo%set_nzeros(l1) + if (debug_level >= psb_debug_outer_)& & write(debug_unit,*) me,' ',trim(name),& diff --git a/base/tools/psb_ssphalo.F90 b/base/tools/psb_ssphalo.F90 index 40bd1195..ecaba85d 100644 --- a/base/tools/psb_ssphalo.F90 +++ b/base/tools/psb_ssphalo.F90 @@ -83,7 +83,7 @@ Subroutine psb_ssphalo(a,desc_a,blk,info,rowcnv,colcnv,& Integer, allocatable :: sdid(:,:), brvindx(:),rvid(:,:), & & rvsz(:), bsdindx(:),sdsz(:), iasnd(:), jasnd(:) real(psb_spk_), allocatable :: valsnd(:) - class(psb_s_base_sparse_mat), allocatable :: acoo + type(psb_s_coo_sparse_mat), allocatable :: acoo integer, pointer :: idxv(:) logical :: rowcnv_,colcnv_,rowscale_,colscale_ character(len=5) :: outfmt_ @@ -164,200 +164,190 @@ Subroutine psb_ssphalo(a,desc_a,blk,info,rowcnv,colcnv,& end select - allocate(psb_s_coo_sparse_mat :: acoo,stat=info) + allocate(acoo,stat=info) if (info /= psb_success_) then info = psb_err_alloc_dealloc_ call psb_errpush(info,name) goto 9999 end if - select type (acoo) - type is (psb_s_coo_sparse_mat) - - - l1 = 0 - - sdsz(:)=0 - rvsz(:)=0 - ipx = 1 - brvindx(ipx) = 0 - bsdindx(ipx) = 0 - counter=1 - idx = 0 - idxs = 0 - idxr = 0 - call acoo%allocate(0,a%get_ncols(),info) - ! For all rows in the halo descriptor, extract and send/receive. - Do - proc=idxv(counter) - if (proc == -1) exit - n_el_recv = idxv(counter+psb_n_elem_recv_) - counter = counter+n_el_recv - n_el_send = idxv(counter+psb_n_elem_send_) - tot_elem = 0 - Do j=0,n_el_send-1 - idx = idxv(counter+psb_elem_send_+j) - n_elem = a%get_nz_row(idx) - tot_elem = tot_elem+n_elem - Enddo - sdsz(proc+1) = tot_elem - call acoo%set_nrows(acoo%get_nrows() + n_el_recv) - counter = counter+n_el_send+3 - Enddo - call mpi_alltoall(sdsz,1,mpi_integer,rvsz,1,mpi_integer,icomm,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='mpi_alltoall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - idxs = 0 - idxr = 0 - counter = 1 - Do - proc=idxv(counter) - if (proc == -1) exit - n_el_recv = idxv(counter+psb_n_elem_recv_) - counter = counter+n_el_recv - n_el_send = idxv(counter+psb_n_elem_send_) - - bsdindx(proc+1) = idxs - idxs = idxs + sdsz(proc+1) - brvindx(proc+1) = idxr - idxr = idxr + rvsz(proc+1) - counter = counter+n_el_send+3 + l1 = 0 + + sdsz(:)=0 + rvsz(:)=0 + ipx = 1 + brvindx(ipx) = 0 + bsdindx(ipx) = 0 + counter=1 + idx = 0 + idxs = 0 + idxr = 0 + call acoo%allocate(0,a%get_ncols(),info) + ! For all rows in the halo descriptor, extract and send/receive. + Do + proc=idxv(counter) + if (proc == -1) exit + n_el_recv = idxv(counter+psb_n_elem_recv_) + counter = counter+n_el_recv + n_el_send = idxv(counter+psb_n_elem_send_) + tot_elem = 0 + Do j=0,n_el_send-1 + idx = idxv(counter+psb_elem_send_+j) + n_elem = a%get_nz_row(idx) + tot_elem = tot_elem+n_elem Enddo + sdsz(proc+1) = tot_elem + call acoo%set_nrows(acoo%get_nrows() + n_el_recv) + counter = counter+n_el_send+3 + Enddo - iszr=sum(rvsz) - call acoo%reallocate(max(iszr,1)) - if (debug_level >= psb_debug_outer_)& - & write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),& - & ' Send:',sdsz(:),' Receive:',rvsz(:) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_sp_reall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - mat_recv = iszr - iszs=sum(sdsz) - call psb_ensure_size(max(iszs,1),iasnd,info) - if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info) - if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info) - - l1 = 0 - ipx = 1 - counter=1 - idx = 0 - - tot_elem=0 - Do - proc=idxv(counter) - if (proc == -1) exit - n_el_recv=idxv(counter+psb_n_elem_recv_) - counter=counter+n_el_recv - n_el_send=idxv(counter+psb_n_elem_send_) - - Do j=0,n_el_send-1 - idx = idxv(counter+psb_elem_send_+j) - n_elem = a%get_nz_row(idx) - call a%csget(idx,idx,ngtz,iasnd,jasnd,valsnd,info,& - & append=.true.,nzin=tot_elem) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_sp_getrow' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - tot_elem=tot_elem+n_elem - Enddo - ipx = ipx + 1 - counter = counter+n_el_send+3 - Enddo - nz = tot_elem - - 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 /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_loc_to_glob' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - - call mpi_alltoallv(valsnd,sdsz,bsdindx,mpi_real,& - & acoo%val,rvsz,brvindx,mpi_real,icomm,info) - call mpi_alltoallv(iasnd,sdsz,bsdindx,mpi_integer,& - & acoo%ia,rvsz,brvindx,mpi_integer,icomm,info) - call mpi_alltoallv(jasnd,sdsz,bsdindx,mpi_integer,& - & acoo%ja,rvsz,brvindx,mpi_integer,icomm,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='mpi_alltoallv' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - ! - ! Convert into local numbering - ! - if (rowcnv_) call psb_glob_to_loc(acoo%ia(1:iszr),desc_a,info,iact='I') - if (colcnv_) call psb_glob_to_loc(acoo%ja(1:iszr),desc_a,info,iact='I') - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psbglob_to_loc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - l1 = 0 - call acoo%set_nrows(0) - ! - irmin = huge(irmin) - icmin = huge(icmin) - irmax = 0 - icmax = 0 - Do i=1,iszr - r=(acoo%ia(i)) - k=(acoo%ja(i)) - ! Just in case some of the conversions were out-of-range - If ((r>0).and.(k>0)) Then - l1=l1+1 - acoo%val(l1) = acoo%val(i) - acoo%ia(l1) = r - acoo%ja(l1) = k - irmin = min(irmin,r) - irmax = max(irmax,r) - icmin = min(icmin,k) - icmax = max(icmax,k) - End If + call mpi_alltoall(sdsz,1,mpi_integer,rvsz,1,mpi_integer,icomm,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='mpi_alltoall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + idxs = 0 + idxr = 0 + counter = 1 + Do + proc=idxv(counter) + if (proc == -1) exit + n_el_recv = idxv(counter+psb_n_elem_recv_) + counter = counter+n_el_recv + n_el_send = idxv(counter+psb_n_elem_send_) + + bsdindx(proc+1) = idxs + idxs = idxs + sdsz(proc+1) + brvindx(proc+1) = idxr + idxr = idxr + rvsz(proc+1) + counter = counter+n_el_send+3 + Enddo + + iszr=sum(rvsz) + call acoo%reallocate(max(iszr,1)) + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),& + & ' Send:',sdsz(:),' Receive:',rvsz(:) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_reall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + mat_recv = iszr + iszs=sum(sdsz) + call psb_ensure_size(max(iszs,1),iasnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info) + + l1 = 0 + ipx = 1 + counter=1 + idx = 0 + + tot_elem=0 + Do + proc=idxv(counter) + if (proc == -1) exit + n_el_recv=idxv(counter+psb_n_elem_recv_) + counter=counter+n_el_recv + n_el_send=idxv(counter+psb_n_elem_send_) + + Do j=0,n_el_send-1 + idx = idxv(counter+psb_elem_send_+j) + n_elem = a%get_nz_row(idx) + call a%csget(idx,idx,ngtz,iasnd,jasnd,valsnd,info,& + & append=.true.,nzin=tot_elem) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_getrow' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + tot_elem=tot_elem+n_elem Enddo - if (rowscale_) then - call acoo%set_nrows(max(irmax-irmin+1,0)) - acoo%ia(1:l1) = acoo%ia(1:l1) - irmin + 1 - else - call acoo%set_nrows(irmax) - end if - if (colscale_) then - call acoo%set_ncols(max(icmax-icmin+1,0)) - acoo%ja(1:l1) = acoo%ja(1:l1) - icmin + 1 - else - call acoo%set_ncols(icmax) - end if - - call acoo%set_nzeros(l1) - - - class default - ! This is impossible - info = psb_err_internal_error_ - call psb_Errpush(info,name) + ipx = ipx + 1 + counter = counter+n_el_send+3 + Enddo + nz = tot_elem + + 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 /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_loc_to_glob' + call psb_errpush(info,name,a_err=ch_err) goto 9999 - end select + end if + + + call mpi_alltoallv(valsnd,sdsz,bsdindx,mpi_real,& + & acoo%val,rvsz,brvindx,mpi_real,icomm,info) + call mpi_alltoallv(iasnd,sdsz,bsdindx,mpi_integer,& + & acoo%ia,rvsz,brvindx,mpi_integer,icomm,info) + call mpi_alltoallv(jasnd,sdsz,bsdindx,mpi_integer,& + & acoo%ja,rvsz,brvindx,mpi_integer,icomm,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='mpi_alltoallv' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! + ! Convert into local numbering + ! + if (rowcnv_) call psb_glob_to_loc(acoo%ia(1:iszr),desc_a,info,iact='I') + if (colcnv_) call psb_glob_to_loc(acoo%ja(1:iszr),desc_a,info,iact='I') + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psbglob_to_loc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + l1 = 0 + call acoo%set_nrows(0) + ! + irmin = huge(irmin) + icmin = huge(icmin) + irmax = 0 + icmax = 0 + Do i=1,iszr + r=(acoo%ia(i)) + k=(acoo%ja(i)) + ! Just in case some of the conversions were out-of-range + If ((r>0).and.(k>0)) Then + l1=l1+1 + acoo%val(l1) = acoo%val(i) + acoo%ia(l1) = r + acoo%ja(l1) = k + irmin = min(irmin,r) + irmax = max(irmax,r) + icmin = min(icmin,k) + icmax = max(icmax,k) + End If + Enddo + if (rowscale_) then + call acoo%set_nrows(max(irmax-irmin+1,0)) + acoo%ia(1:l1) = acoo%ia(1:l1) - irmin + 1 + else + call acoo%set_nrows(irmax) + end if + if (colscale_) then + call acoo%set_ncols(max(icmax-icmin+1,0)) + acoo%ja(1:l1) = acoo%ja(1:l1) - icmin + 1 + else + call acoo%set_ncols(icmax) + end if + + call acoo%set_nzeros(l1) + if (debug_level >= psb_debug_outer_)& & write(debug_unit,*) me,' ',trim(name),& diff --git a/base/tools/psb_zsphalo.F90 b/base/tools/psb_zsphalo.F90 index 30437952..38274115 100644 --- a/base/tools/psb_zsphalo.F90 +++ b/base/tools/psb_zsphalo.F90 @@ -83,7 +83,7 @@ Subroutine psb_zsphalo(a,desc_a,blk,info,rowcnv,colcnv,& Integer, allocatable :: sdid(:,:), brvindx(:),rvid(:,:), & & rvsz(:), bsdindx(:),sdsz(:), iasnd(:), jasnd(:) complex(psb_dpk_), allocatable :: valsnd(:) - class(psb_z_base_sparse_mat), allocatable :: acoo + type(psb_z_coo_sparse_mat), allocatable :: acoo integer, pointer :: idxv(:) logical :: rowcnv_,colcnv_,rowscale_,colscale_ character(len=5) :: outfmt_ @@ -164,200 +164,190 @@ Subroutine psb_zsphalo(a,desc_a,blk,info,rowcnv,colcnv,& end select - allocate(psb_z_coo_sparse_mat :: acoo,stat=info) + allocate(acoo,stat=info) if (info /= psb_success_) then info = psb_err_alloc_dealloc_ call psb_errpush(info,name) goto 9999 end if - select type (acoo) - type is (psb_z_coo_sparse_mat) - - - l1 = 0 - - sdsz(:)=0 - rvsz(:)=0 - ipx = 1 - brvindx(ipx) = 0 - bsdindx(ipx) = 0 - counter=1 - idx = 0 - idxs = 0 - idxr = 0 - call acoo%allocate(0,a%get_ncols(),info) - ! For all rows in the halo descriptor, extract and send/receive. - Do - proc=idxv(counter) - if (proc == -1) exit - n_el_recv = idxv(counter+psb_n_elem_recv_) - counter = counter+n_el_recv - n_el_send = idxv(counter+psb_n_elem_send_) - tot_elem = 0 - Do j=0,n_el_send-1 - idx = idxv(counter+psb_elem_send_+j) - n_elem = a%get_nz_row(idx) - tot_elem = tot_elem+n_elem - Enddo - sdsz(proc+1) = tot_elem - call acoo%set_nrows(acoo%get_nrows() + n_el_recv) - counter = counter+n_el_send+3 - Enddo - call mpi_alltoall(sdsz,1,mpi_integer,rvsz,1,mpi_integer,icomm,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='mpi_alltoall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - idxs = 0 - idxr = 0 - counter = 1 - Do - proc=idxv(counter) - if (proc == -1) exit - n_el_recv = idxv(counter+psb_n_elem_recv_) - counter = counter+n_el_recv - n_el_send = idxv(counter+psb_n_elem_send_) - - bsdindx(proc+1) = idxs - idxs = idxs + sdsz(proc+1) - brvindx(proc+1) = idxr - idxr = idxr + rvsz(proc+1) - counter = counter+n_el_send+3 + l1 = 0 + + sdsz(:)=0 + rvsz(:)=0 + ipx = 1 + brvindx(ipx) = 0 + bsdindx(ipx) = 0 + counter=1 + idx = 0 + idxs = 0 + idxr = 0 + call acoo%allocate(0,a%get_ncols(),info) + ! For all rows in the halo descriptor, extract and send/receive. + Do + proc=idxv(counter) + if (proc == -1) exit + n_el_recv = idxv(counter+psb_n_elem_recv_) + counter = counter+n_el_recv + n_el_send = idxv(counter+psb_n_elem_send_) + tot_elem = 0 + Do j=0,n_el_send-1 + idx = idxv(counter+psb_elem_send_+j) + n_elem = a%get_nz_row(idx) + tot_elem = tot_elem+n_elem Enddo + sdsz(proc+1) = tot_elem + call acoo%set_nrows(acoo%get_nrows() + n_el_recv) + counter = counter+n_el_send+3 + Enddo - iszr=sum(rvsz) - call acoo%reallocate(max(iszr,1)) - if (debug_level >= psb_debug_outer_)& - & write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),& - & ' Send:',sdsz(:),' Receive:',rvsz(:) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_sp_reall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - mat_recv = iszr - iszs=sum(sdsz) - call psb_ensure_size(max(iszs,1),iasnd,info) - if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info) - if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info) - - l1 = 0 - ipx = 1 - counter=1 - idx = 0 - - tot_elem=0 - Do - proc=idxv(counter) - if (proc == -1) exit - n_el_recv=idxv(counter+psb_n_elem_recv_) - counter=counter+n_el_recv - n_el_send=idxv(counter+psb_n_elem_send_) - - Do j=0,n_el_send-1 - idx = idxv(counter+psb_elem_send_+j) - n_elem = a%get_nz_row(idx) - call a%csget(idx,idx,ngtz,iasnd,jasnd,valsnd,info,& - & append=.true.,nzin=tot_elem) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_sp_getrow' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - tot_elem=tot_elem+n_elem - Enddo - ipx = ipx + 1 - counter = counter+n_el_send+3 - Enddo - nz = tot_elem - - 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 /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_loc_to_glob' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - - call mpi_alltoallv(valsnd,sdsz,bsdindx,mpi_double_complex,& - & acoo%val,rvsz,brvindx,mpi_double_complex,icomm,info) - call mpi_alltoallv(iasnd,sdsz,bsdindx,mpi_integer,& - & acoo%ia,rvsz,brvindx,mpi_integer,icomm,info) - call mpi_alltoallv(jasnd,sdsz,bsdindx,mpi_integer,& - & acoo%ja,rvsz,brvindx,mpi_integer,icomm,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='mpi_alltoallv' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - ! - ! Convert into local numbering - ! - if (rowcnv_) call psb_glob_to_loc(acoo%ia(1:iszr),desc_a,info,iact='I') - if (colcnv_) call psb_glob_to_loc(acoo%ja(1:iszr),desc_a,info,iact='I') - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psbglob_to_loc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - l1 = 0 - call acoo%set_nrows(0) - ! - irmin = huge(irmin) - icmin = huge(icmin) - irmax = 0 - icmax = 0 - Do i=1,iszr - r=(acoo%ia(i)) - k=(acoo%ja(i)) - ! Just in case some of the conversions were out-of-range - If ((r>0).and.(k>0)) Then - l1=l1+1 - acoo%val(l1) = acoo%val(i) - acoo%ia(l1) = r - acoo%ja(l1) = k - irmin = min(irmin,r) - irmax = max(irmax,r) - icmin = min(icmin,k) - icmax = max(icmax,k) - End If + call mpi_alltoall(sdsz,1,mpi_integer,rvsz,1,mpi_integer,icomm,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='mpi_alltoall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + idxs = 0 + idxr = 0 + counter = 1 + Do + proc=idxv(counter) + if (proc == -1) exit + n_el_recv = idxv(counter+psb_n_elem_recv_) + counter = counter+n_el_recv + n_el_send = idxv(counter+psb_n_elem_send_) + + bsdindx(proc+1) = idxs + idxs = idxs + sdsz(proc+1) + brvindx(proc+1) = idxr + idxr = idxr + rvsz(proc+1) + counter = counter+n_el_send+3 + Enddo + + iszr=sum(rvsz) + call acoo%reallocate(max(iszr,1)) + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),& + & ' Send:',sdsz(:),' Receive:',rvsz(:) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_reall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + mat_recv = iszr + iszs=sum(sdsz) + call psb_ensure_size(max(iszs,1),iasnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info) + + l1 = 0 + ipx = 1 + counter=1 + idx = 0 + + tot_elem=0 + Do + proc=idxv(counter) + if (proc == -1) exit + n_el_recv=idxv(counter+psb_n_elem_recv_) + counter=counter+n_el_recv + n_el_send=idxv(counter+psb_n_elem_send_) + + Do j=0,n_el_send-1 + idx = idxv(counter+psb_elem_send_+j) + n_elem = a%get_nz_row(idx) + call a%csget(idx,idx,ngtz,iasnd,jasnd,valsnd,info,& + & append=.true.,nzin=tot_elem) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_getrow' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + tot_elem=tot_elem+n_elem Enddo - if (rowscale_) then - call acoo%set_nrows(max(irmax-irmin+1,0)) - acoo%ia(1:l1) = acoo%ia(1:l1) - irmin + 1 - else - call acoo%set_nrows(irmax) - end if - if (colscale_) then - call acoo%set_ncols(max(icmax-icmin+1,0)) - acoo%ja(1:l1) = acoo%ja(1:l1) - icmin + 1 - else - call acoo%set_ncols(icmax) - end if - - call acoo%set_nzeros(l1) - - - class default - ! This is impossible - info = psb_err_internal_error_ - call psb_Errpush(info,name) + ipx = ipx + 1 + counter = counter+n_el_send+3 + Enddo + nz = tot_elem + + 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 /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_loc_to_glob' + call psb_errpush(info,name,a_err=ch_err) goto 9999 - end select + end if + + + call mpi_alltoallv(valsnd,sdsz,bsdindx,mpi_double_complex,& + & acoo%val,rvsz,brvindx,mpi_double_complex,icomm,info) + call mpi_alltoallv(iasnd,sdsz,bsdindx,mpi_integer,& + & acoo%ia,rvsz,brvindx,mpi_integer,icomm,info) + call mpi_alltoallv(jasnd,sdsz,bsdindx,mpi_integer,& + & acoo%ja,rvsz,brvindx,mpi_integer,icomm,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='mpi_alltoallv' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! + ! Convert into local numbering + ! + if (rowcnv_) call psb_glob_to_loc(acoo%ia(1:iszr),desc_a,info,iact='I') + if (colcnv_) call psb_glob_to_loc(acoo%ja(1:iszr),desc_a,info,iact='I') + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psbglob_to_loc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + l1 = 0 + call acoo%set_nrows(0) + ! + irmin = huge(irmin) + icmin = huge(icmin) + irmax = 0 + icmax = 0 + Do i=1,iszr + r=(acoo%ia(i)) + k=(acoo%ja(i)) + ! Just in case some of the conversions were out-of-range + If ((r>0).and.(k>0)) Then + l1=l1+1 + acoo%val(l1) = acoo%val(i) + acoo%ia(l1) = r + acoo%ja(l1) = k + irmin = min(irmin,r) + irmax = max(irmax,r) + icmin = min(icmin,k) + icmax = max(icmax,k) + End If + Enddo + if (rowscale_) then + call acoo%set_nrows(max(irmax-irmin+1,0)) + acoo%ia(1:l1) = acoo%ia(1:l1) - irmin + 1 + else + call acoo%set_nrows(irmax) + end if + if (colscale_) then + call acoo%set_ncols(max(icmax-icmin+1,0)) + acoo%ja(1:l1) = acoo%ja(1:l1) - icmin + 1 + else + call acoo%set_ncols(icmax) + end if + + call acoo%set_nzeros(l1) + if (debug_level >= psb_debug_outer_)& & write(debug_unit,*) me,' ',trim(name),&