From cbc36fdd481d8fc9dc7f38d71f236139999207f2 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Thu, 4 May 2006 12:00:46 +0000 Subject: [PATCH] *** empty log message *** --- src/modules/psb_serial_mod.f90 | 10 +- src/serial/psb_dcoins.f90 | 474 +++++++++++++++------ src/serial/psb_update_mod.f90 | 729 +++++++++++++++++++++++++-------- src/serial/psb_zcoins.f90 | 475 ++++++++++++++------- src/tools/psb_dspins.f90 | 6 +- src/tools/psb_zspins.f90 | 6 +- 6 files changed, 1244 insertions(+), 456 deletions(-) diff --git a/src/modules/psb_serial_mod.f90 b/src/modules/psb_serial_mod.f90 index ead82596..421adcb9 100644 --- a/src/modules/psb_serial_mod.f90 +++ b/src/modules/psb_serial_mod.f90 @@ -310,22 +310,24 @@ module psb_serial_mod end interface interface psb_coins - subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild) + subroutine psb_dcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild) use psb_spmat_type integer, intent(in) :: nz, imin,imax,jmin,jmax - integer, intent(in) :: ia(:),ja(:),gtl(:) + integer, intent(in) :: ia(:),ja(:) real(kind(1.d0)), intent(in) :: val(:) type(psb_dspmat_type), intent(inout) :: a integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) logical, optional, intent(in) :: rebuild end subroutine psb_dcoins - subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild) + subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild) use psb_spmat_type integer, intent(in) :: nz, imin,imax,jmin,jmax - integer, intent(in) :: ia(:),ja(:),gtl(:) + integer, intent(in) :: ia(:),ja(:) complex(kind(1.d0)), intent(in) :: val(:) type(psb_zspmat_type), intent(inout) :: a integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) logical, optional, intent(in) :: rebuild end subroutine psb_zcoins end interface diff --git a/src/serial/psb_dcoins.f90 b/src/serial/psb_dcoins.f90 index 562d5a5c..23c09a80 100644 --- a/src/serial/psb_dcoins.f90 +++ b/src/serial/psb_dcoins.f90 @@ -31,7 +31,7 @@ ! File: psbdcoins.f90 ! Subroutine: ! Parameters: -subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild) +subroutine psb_dcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild) use psb_spmat_type use psb_const_mod @@ -43,10 +43,11 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild) implicit none integer, intent(in) :: nz, imin,imax,jmin,jmax - integer, intent(in) :: ia(:),ja(:),gtl(:) + integer, intent(in) :: ia(:),ja(:) real(kind(1.d0)), intent(in) :: val(:) type(psb_dspmat_type), intent(inout) :: a integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) logical, intent(in), optional :: rebuild character(len=5) :: ufida @@ -95,153 +96,312 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild) end if call touppers(a%fida,ufida) - ng = size(gtl) spstate = psb_sp_getifld(psb_state_,a,info) - select case(spstate) - case(psb_spmat_bld_) - if ((ufida /= 'COO').and.(ufida/='COI')) then - info = 134 - ch_err(1:3)=ufida(1:3) - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - call psb_spinfo(psb_nztotreq_,a,nza,info) - call psb_spinfo(psb_nzsizereq_,a,isza,info) - if(info /= izero) then - info=4010 - ch_err='psb_spinfo' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - endif + if (present(gtl)) then - if ((nza+nz)>isza) then - call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info) + ng = size(gtl) + select case(spstate) + case(psb_spmat_bld_) + if ((ufida /= 'COO').and.(ufida/='COI')) then + info = 134 + ch_err(1:3)=ufida(1:3) + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_spinfo(psb_nztotreq_,a,nza,info) + call psb_spinfo(psb_nzsizereq_,a,isza,info) if(info /= izero) then info=4010 - ch_err='psb_sp_reall' + ch_err='psb_spinfo' call psb_errpush(info,name,a_err=ch_err) goto 9999 endif - endif - call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,gtl,ng,& - & imin,imax,jmin,jmax,info) - if(info /= izero) then - info=4010 - ch_err='psb_inner_ins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - endif - if (debug) then - if ((nza - a%infoa(psb_nnz_)) /= nz) then - write(0,*) 'PSB_COINS: insert discarded items ' + + if ((nza+nz)>isza) then + call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info) + if(info /= izero) then + info=4010 + ch_err='psb_sp_reall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + endif + + call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,& + & imin,imax,jmin,jmax,info,gtl,ng) + + if(info /= izero) then + info=4010 + ch_err='psb_inner_ins' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + if (debug) then + if ((nza - a%infoa(psb_nnz_)) /= nz) then + write(0,*) 'PSB_COINS: insert discarded items ' + end if end if - end if - if ((nza - psb_sp_getifld(psb_nnz_,a,info)) /= nz) then - call psb_sp_setifld(nza,psb_del_bnd_,a,info) - endif - call psb_sp_setifld(nza,psb_nnz_,a,info) + if ((nza - psb_sp_getifld(psb_nnz_,a,info)) /= nz) then + call psb_sp_setifld(nza,psb_del_bnd_,a,info) + endif + call psb_sp_setifld(nza,psb_nnz_,a,info) + + case(psb_spmat_upd_) + + iupd = psb_sp_getifld(psb_upd_,a,info) + select case (iupd) + case (psb_upd_perm_) + ip1 = psb_sp_getifld(psb_upd_pnt_,a,info) + nzl = psb_sp_getifld(psb_del_bnd_,a,info) + nza = a%ia2(ip1+psb_nnz_) + + + call psb_inner_upd(nz,ia,ja,val,nza,a%aspk,& + & imin,imax,jmin,jmax,nzl,info,gtl,ng) + + if (info /= izero) then + info=4010 + ch_err='psb_inner_upd' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + if (debug) then + if ((nza - a%ia2(ip1+psb_nnz_)) /= nz) then + write(0,*) 'PSB_COINS: update discarded items ' + end if + end if + a%ia2(ip1+psb_nnz_) = nza + + if (debug) write(0,*) 'From COINS(UPD) : NZA:',nza + + case (psb_upd_dflt_, psb_upd_srch_) + + call psb_srch_upd(nz,ia,ja,val,nza,a,& + & imin,imax,jmin,jmax,nzl,info,gtl,ng) + + if (info > 0) then + if (rebuild_) then + if (debug) write(0,*)& + & 'COINS: Going through rebuild_ fingers crossed!' + irst = info + call psb_nullify_sp(tmp) + tmp%fida='COO' + call psb_csdp(a,tmp,info) + call psb_sp_setifld(psb_spmat_bld_,psb_state_,tmp,info) + if (debug) then + write(0,*) 'COINS Rebuild: size',tmp%infoa(psb_nnz_) ,irst + endif + call psb_sp_transfer(tmp,a,info) + call psb_spinfo(psb_nztotreq_,a,nza,info) + call psb_spinfo(psb_nzsizereq_,a,isza,info) + if(info /= izero) then + info=4010 + ch_err='psb_spinfo' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + + if (debug) write(0,*)& + & 'COINS: Reinserting',a%fida,nza,isza + if ((nza+nz)>isza) then + call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info) + if(info /= izero) then + info=4010 + ch_err='psb_sp_reall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + endif + if (irst <= nz) then + call psb_inner_ins((nz-irst+1),ia(irst:nz),ja(irst:nz),val(irst:nz),& + & nza,a%ia1,a%ia2,a%aspk,imin,imax,jmin,jmax,info,gtl,ng) + call psb_sp_setifld(nza,psb_del_bnd_,a,info) + call psb_sp_setifld(nza,psb_nnz_,a,info) + end if + + else + info = 2231 + call psb_errpush(info,name) + goto 9999 + end if + else if (info < 0) then + info = 2230 + call psb_errpush(info,name) + goto 9999 + + end if + + case default + + info = 2231 + call psb_errpush(info,name) + goto 9999 + end select + + + case default + info = 2232 + call psb_errpush(info,name) + goto 9999 + end select + + - case(psb_spmat_upd_) - iupd = psb_sp_getifld(psb_upd_,a,info) - select case (iupd) - case (psb_upd_perm_) - ip1 = psb_sp_getifld(psb_upd_pnt_,a,info) - nzl = psb_sp_getifld(psb_del_bnd_,a,info) - nza = a%ia2(ip1+psb_nnz_) - nza = a%ia2(ip1+psb_nnz_) - nzl = a%infoa(psb_del_bnd_) + else - call psb_inner_upd(nz,ia,ja,val,nza,a%aspk,gtl,ng,& - & imin,imax,jmin,jmax,nzl,info) + ng = -1 - if (info /= izero) then + select case(spstate) + case(psb_spmat_bld_) + if ((ufida /= 'COO').and.(ufida/='COI')) then + info = 134 + ch_err(1:3)=ufida(1:3) + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_spinfo(psb_nztotreq_,a,nza,info) + call psb_spinfo(psb_nzsizereq_,a,isza,info) + if(info /= izero) then info=4010 - ch_err='psb_inner_upd' + ch_err='psb_spinfo' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + + if ((nza+nz)>isza) then + call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info) + if(info /= izero) then + info=4010 + ch_err='psb_sp_reall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + endif + + call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,& + & imin,imax,jmin,jmax,info) + + if(info /= izero) then + info=4010 + ch_err='psb_inner_ins' call psb_errpush(info,name,a_err=ch_err) goto 9999 endif if (debug) then - if ((nza - a%ia2(ip1+psb_nnz_)) /= nz) then - write(0,*) 'PSB_COINS: update discarded items ' + if ((nza - a%infoa(psb_nnz_)) /= nz) then + write(0,*) 'PSB_COINS: insert discarded items ' end if end if - a%ia2(ip1+psb_nnz_) = nza + if ((nza - psb_sp_getifld(psb_nnz_,a,info)) /= nz) then + call psb_sp_setifld(nza,psb_del_bnd_,a,info) + endif + call psb_sp_setifld(nza,psb_nnz_,a,info) - if (debug) write(0,*) 'From COINS(UPD) : NZA:',nza + case(psb_spmat_upd_) - case (psb_upd_dflt_, psb_upd_srch_) + iupd = psb_sp_getifld(psb_upd_,a,info) + select case (iupd) + case (psb_upd_perm_) + ip1 = psb_sp_getifld(psb_upd_pnt_,a,info) + nzl = psb_sp_getifld(psb_del_bnd_,a,info) + nza = a%ia2(ip1+psb_nnz_) - call psb_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,& + + call psb_inner_upd(nz,ia,ja,val,nza,a%aspk,& + & imin,imax,jmin,jmax,nzl,info) + + if (info /= izero) then + info=4010 + ch_err='psb_inner_upd' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + if (debug) then + if ((nza - a%ia2(ip1+psb_nnz_)) /= nz) then + write(0,*) 'PSB_COINS: update discarded items ' + end if + end if + a%ia2(ip1+psb_nnz_) = nza + + if (debug) write(0,*) 'From COINS(UPD) : NZA:',nza + + case (psb_upd_dflt_, psb_upd_srch_) + + call psb_srch_upd(nz,ia,ja,val,nza,a,& & imin,imax,jmin,jmax,nzl,info) - if (info > 0) then - if (rebuild_) then - if (debug) write(0,*)& - & 'COINS: Going through rebuild_ fingers crossed!' - irst = info - call psb_nullify_sp(tmp) - tmp%fida='COO' - call psb_csdp(a,tmp,info) - call psb_sp_setifld(psb_spmat_bld_,psb_state_,tmp,info) - if (debug) then - write(0,*) 'COINS Rebuild: size',tmp%infoa(psb_nnz_) ,irst - endif - call psb_sp_transfer(tmp,a,info) - call psb_spinfo(psb_nztotreq_,a,nza,info) - call psb_spinfo(psb_nzsizereq_,a,isza,info) - if(info /= izero) then - info=4010 - ch_err='psb_spinfo' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - endif - - if (debug) write(0,*)& - & 'COINS: Reinserting',a%fida,nza,isza - if ((nza+nz)>isza) then - call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info) + + if (info > 0) then + if (rebuild_) then + if (debug) write(0,*)& + & 'COINS: Going through rebuild_ fingers crossed!' + irst = info + call psb_nullify_sp(tmp) + tmp%fida='COO' + call psb_csdp(a,tmp,info) + call psb_sp_setifld(psb_spmat_bld_,psb_state_,tmp,info) + if (debug) then + write(0,*) 'COINS Rebuild: size',tmp%infoa(psb_nnz_) ,irst + endif + call psb_sp_transfer(tmp,a,info) + call psb_spinfo(psb_nztotreq_,a,nza,info) + call psb_spinfo(psb_nzsizereq_,a,isza,info) if(info /= izero) then info=4010 - ch_err='psb_sp_reall' + ch_err='psb_spinfo' call psb_errpush(info,name,a_err=ch_err) goto 9999 endif - endif - if (irst <= nz) then - call psb_inner_ins((nz-irst+1),ia(irst:nz),ja(irst:nz),val(irst:nz),& - & nza,a%ia1,a%ia2,a%aspk,gtl,ng,imin,imax,jmin,jmax,info) - call psb_sp_setifld(nza,psb_del_bnd_,a,info) - call psb_sp_setifld(nza,psb_nnz_,a,info) + + if (debug) write(0,*)& + & 'COINS: Reinserting',a%fida,nza,isza + if ((nza+nz)>isza) then + call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info) + if(info /= izero) then + info=4010 + ch_err='psb_sp_reall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + endif + if (irst <= nz) then + call psb_inner_ins((nz-irst+1),ia(irst:nz),ja(irst:nz),val(irst:nz),& + & nza,a%ia1,a%ia2,a%aspk,imin,imax,jmin,jmax,info) + call psb_sp_setifld(nza,psb_del_bnd_,a,info) + call psb_sp_setifld(nza,psb_nnz_,a,info) + end if + + else + info = 2231 + call psb_errpush(info,name) + goto 9999 end if - - else - info = 2231 + else if (info < 0) then + info = 2230 call psb_errpush(info,name) goto 9999 + end if - else if (info < 0) then - info = 2230 + + case default + + info = 2231 call psb_errpush(info,name) goto 9999 + end select - end if - - case default - info = 2231 + case default + info = 2232 call psb_errpush(info,name) - goto 9999 + goto 9999 end select + endif - case default - info = 2232 - call psb_errpush(info,name) - goto 9999 - end select return call psb_erractionrestore(err_act) @@ -257,70 +417,112 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild) contains - subroutine psb_inner_upd(nz,ia,ja,val,nza,aspk,gtl,ng,& - & imin,imax,jmin,jmax,nzl,info) + subroutine psb_inner_upd(nz,ia,ja,val,nza,aspk,& + & imin,imax,jmin,jmax,nzl,info,gtl,ng) implicit none - integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng - integer, intent(in) :: ia(*),ja(*),gtl(*) + integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl + integer, intent(in) :: ia(*),ja(*) integer, intent(inout) :: nza real(kind(1.d0)), intent(in) :: val(*) real(kind(1.d0)), intent(inout) :: aspk(*) integer, intent(out) :: info + integer, intent(in), optional :: ng,gtl(*) integer :: i,ir,ic - if (nza >= nzl) then - do i=1, nz - nza = nza + 1 - aspk(nza) = val(i) - end do + if (present(gtl)) then + if (.not.present(ng)) then + info = -1 + return + endif + if (nza >= nzl) then + do i=1, nz + nza = nza + 1 + aspk(nza) = val(i) + end do + else + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then + nza = nza + 1 + aspk(nza) = val(i) + end if + end if + end do + end if else - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - ic = gtl(ic) + if (nza >= nzl) then + do i=1, nz + nza = nza + 1 + aspk(nza) = val(i) + end do + else + do i=1, nz + ir = ia(i) + ic = ja(i) if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then nza = nza + 1 aspk(nza) = val(i) end if - end if - end do + end do + end if end if - end subroutine psb_inner_upd - subroutine psb_inner_ins(nz,ia,ja,val,nza,ia1,ia2,aspk,gtl,ng,& - & imin,imax,jmin,jmax,info) + subroutine psb_inner_ins(nz,ia,ja,val,nza,ia1,ia2,aspk,& + & imin,imax,jmin,jmax,info,gtl,ng) implicit none - integer, intent(in) :: nz, imin,imax,jmin,jmax,ng - integer, intent(in) :: ia(*),ja(*),gtl(*) + integer, intent(in) :: nz, imin,imax,jmin,jmax + integer, intent(in) :: ia(*),ja(*) integer, intent(inout) :: nza,ia1(*),ia2(*) real(kind(1.d0)), intent(in) :: val(*) real(kind(1.d0)), intent(inout) :: aspk(*) integer, intent(out) :: info + integer, intent(in), optional :: ng,gtl(*) integer :: i,ir,ic info = 0 - do i=1, nz - ir = ia(i) - ic = ja(i) + if (present(gtl)) then + if (.not.present(ng)) then + info = -1 + return + endif + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then + nza = nza + 1 + ia1(nza) = ir + ia2(nza) = ic + aspk(nza) = val(i) + end if + end if + end do + else - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - ic = gtl(ic) + do i=1, nz + ir = ia(i) + ic = ja(i) if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then nza = nza + 1 ia1(nza) = ir ia2(nza) = ic aspk(nza) = val(i) end if - end if - end do + end do + end if end subroutine psb_inner_ins + + end subroutine psb_dcoins diff --git a/src/serial/psb_update_mod.f90 b/src/serial/psb_update_mod.f90 index 2201f609..828033de 100644 --- a/src/serial/psb_update_mod.f90 +++ b/src/serial/psb_update_mod.f90 @@ -51,79 +51,129 @@ module psb_update_mod contains - subroutine psb_d_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,& - & imin,imax,jmin,jmax,nzl,info) + subroutine psb_d_srch_upd(nz,ia,ja,val,nza,a,& + & imin,imax,jmin,jmax,nzl,info,gtl,ng) implicit none type(psb_dspmat_type), intent(inout) :: a - integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng - integer, intent(in) :: ia(*),ja(*),gtl(*) + integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl + integer, intent(in) :: ia(*),ja(*) integer, intent(inout) :: nza real(kind(1.d0)), intent(in) :: val(*) integer, intent(out) :: info + integer, intent(in), optional :: ng,gtl(*) info = 0 - select case(toupper(a%fida)) - case ('CSR') + if (present(gtl)) then + if (.not.present(ng)) then + info = -1 + return + endif + + select case(toupper(a%fida)) + case ('CSR') !!$ write(0,*) 'Calling csr_srch_upd' - call csr_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,& - & imin,imax,jmin,jmax,nzl,info) + call csr_srch_upd(nz,ia,ja,val,nza,a,& + & imin,imax,jmin,jmax,nzl,info,gtl,ng) !!$ write(0,*) 'From csr_srch_upd:',info - case ('COO') - call coo_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,& - & imin,imax,jmin,jmax,nzl,info) - - case default - - info = -9 - - end select + case ('COO') + call coo_srch_upd(nz,ia,ja,val,nza,a,& + & imin,imax,jmin,jmax,nzl,info,gtl,ng) + + case default + + info = -9 + + end select + else + select case(toupper(a%fida)) + case ('CSR') +!!$ write(0,*) 'Calling csr_srch_upd' + call csr_srch_upd(nz,ia,ja,val,nza,a,& + & imin,imax,jmin,jmax,nzl,info) +!!$ write(0,*) 'From csr_srch_upd:',info + case ('COO') + call coo_srch_upd(nz,ia,ja,val,nza,a,& + & imin,imax,jmin,jmax,nzl,info) + + case default + + info = -9 + + end select + end if end subroutine psb_d_srch_upd - subroutine psb_z_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,& - & imin,imax,jmin,jmax,nzl,info) + subroutine psb_z_srch_upd(nz,ia,ja,val,nza,a,& + & imin,imax,jmin,jmax,nzl,info,gtl,ng) implicit none type(psb_zspmat_type), intent(inout) :: a - integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng - integer, intent(in) :: ia(*),ja(*),gtl(*) + integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl + integer, intent(in) :: ia(*),ja(*) integer, intent(inout) :: nza complex(kind(1.d0)), intent(in) :: val(*) integer, intent(out) :: info - + integer, intent(in), optional :: ng,gtl(*) + info = 0 - select case(toupper(a%fida)) - case ('CSR') + if (present(gtl)) then + if (.not.present(ng)) then + info = -1 + return + endif + + select case(toupper(a%fida)) + case ('CSR') !!$ write(0,*) 'Calling csr_srch_upd' - call csr_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,& - & imin,imax,jmin,jmax,nzl,info) + call csr_srch_upd(nz,ia,ja,val,nza,a,& + & imin,imax,jmin,jmax,nzl,info,gtl,ng) !!$ write(0,*) 'From csr_srch_upd:',info - case ('COO') - call coo_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,& - & imin,imax,jmin,jmax,nzl,info) - - case default - - info = -9 - - end select + case ('COO') + call coo_srch_upd(nz,ia,ja,val,nza,a,& + & imin,imax,jmin,jmax,nzl,info,gtl,ng) + + case default + + info = -9 + + end select + else + select case(toupper(a%fida)) + case ('CSR') +!!$ write(0,*) 'Calling csr_srch_upd' + call csr_srch_upd(nz,ia,ja,val,nza,a,& + & imin,imax,jmin,jmax,nzl,info) +!!$ write(0,*) 'From csr_srch_upd:',info + case ('COO') + call coo_srch_upd(nz,ia,ja,val,nza,a,& + & imin,imax,jmin,jmax,nzl,info) + + case default + + info = -9 + + end select + end if end subroutine psb_z_srch_upd - subroutine d_csr_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,& - & imin,imax,jmin,jmax,nzl,info) + subroutine d_csr_srch_upd(nz,ia,ja,val,nza,a,& + & imin,imax,jmin,jmax,nzl,info,gtl,ng) implicit none type(psb_dspmat_type), intent(inout) :: a - integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng - integer, intent(in) :: ia(*),ja(*),gtl(*) + integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl + integer, intent(in) :: ia(*),ja(*) integer, intent(inout) :: nza real(kind(1.d0)), intent(in) :: val(*) integer, intent(out) :: info + integer, intent(in), optional :: ng,gtl(*) + integer :: i,ir,ic,check_flag, ilr, ilc, ip, & & i1,i2,nc,lb,ub,m,nnz,dupl @@ -131,19 +181,112 @@ contains dupl = psb_sp_getifld(psb_dupl_,a,info) - select case(dupl) - case(psb_dupl_ovwrt_,psb_dupl_err_) - ! Overwrite. - ! Cannot test for error, should have been caught earlier. - - ilr = -1 - ilc = -1 - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - ic = gtl(ic) + + if (present(gtl)) then + if (.not.present(ng)) then + info = -1 + return + endif + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + i1 = a%ia2(ir) + i2 = a%ia2(ir+1) + nc=i2-i1 + + + if (.true.) then + call issrch(ip,ic,nc,a%ia1(i1:i2-1)) + if (ip>0) then + a%aspk(i1+ip-1) = val(i) + else + write(0,*)'Was searching ',ic,' in: ',i1,i2,' : ',a%ia1(i1:i2-1) + info = i + return + end if + + else +!!$ + ip = -1 + lb = i1 + ub = i2-1 + do + if (lb > ub) exit + m = (lb+ub)/2 +!!$ write(0,*) 'Debug: ',lb,m,ub + if (ic == a%ia1(m)) then + ip = m + lb = ub + 1 + else if (ic < a%ia1(m)) then + ub = m-1 + else + lb = m + 1 + end if + enddo + + if (ip>0) then + a%aspk(ip) = val(i) + else + write(0,*)'Was searching ',ic,' in: ',i1,i2,' : ',a%ia1(i1:i2-1) + info = i + return + end if + + end if + end if + end do + + case(psb_dupl_add_) + ! Add + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + i1 = a%ia2(ir) + i2 = a%ia2(ir+1) + nc = i2-i1 + call issrch(ip,ic,nc,a%ia1(i1:i2-1)) + if (ip>0) then + a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) + else + info = i + return + end if + end if + end do + + case default + info = -3 + write(0,*) 'Duplicate handling: ',dupl + end select + + else + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) i1 = a%ia2(ir) i2 = a%ia2(ir+1) nc=i2-i1 @@ -187,19 +330,15 @@ contains end if end if - end if - end do - - case(psb_dupl_add_) - ! Add - ilr = -1 - ilc = -1 - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - ic = gtl(ic) + end do + + case(psb_dupl_add_) + ! Add + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) i1 = a%ia2(ir) i2 = a%ia2(ir+1) nc = i2-i1 @@ -210,25 +349,27 @@ contains info = i return end if - end if - end do + end do + + case default + info = -3 + write(0,*) 'Duplicate handling: ',dupl + end select + end if - case default - info = -3 - write(0,*) 'Duplicate handling: ',dupl - end select end subroutine d_csr_srch_upd - subroutine d_coo_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,& - & imin,imax,jmin,jmax,nzl,info) + subroutine d_coo_srch_upd(nz,ia,ja,val,nza,a,& + & imin,imax,jmin,jmax,nzl,info,gtl,ng) implicit none type(psb_dspmat_type), intent(inout) :: a - integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng - integer, intent(in) :: ia(*),ja(*),gtl(*) + integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl + integer, intent(in) :: ia(*),ja(*) integer, intent(inout) :: nza real(kind(1.d0)), intent(in) :: val(*) integer, intent(out) :: info + integer, intent(in), optional :: ng,gtl(*) integer :: i,ir,ic,check_flag, ilr, ilc, ip, & & i1,i2,nc,lb,ub,m,nnz,dupl,isrt @@ -246,16 +387,95 @@ contains nnz = psb_sp_getifld(psb_nnz_,a,info) - select case(dupl) - case(psb_dupl_ovwrt_,psb_dupl_err_) - ! Overwrite. - ! Cannot test for error, should have been caught earlier. - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - ic = gtl(ic) + if (present(gtl)) then + if (.not.present(ng)) then + info = -1 + return + endif + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if (ir /= ilr) then + call ibsrch(i1,ir,nnz,a%ia1) + i2 = i1 + do + if (i2+1 > nnz) exit + if (a%ia1(i2+1) /= a%ia1(i2)) exit + i2 = i2 + 1 + end do + do + if (i1-1 < 1) exit + if (a%ia1(i1-1) /= a%ia1(i1)) exit + i1 = i1 - 1 + end do + ilr = ir + end if + nc = i2-i1+1 + call issrch(ip,ic,nc,a%ia2(i1:i2)) + if (ip>0) then + a%aspk(i1+ip-1) = val(i) + else + info = i + return + end if + end if + end do + case(psb_dupl_add_) + ! Add + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if (ir /= ilr) then + call ibsrch(i1,ir,nnz,a%ia1) + i2 = i1 + do + if (i2+1 > nnz) exit + if (a%ia1(i2+1) /= a%ia1(i2)) exit + i2 = i2 + 1 + end do + do + if (i1-1 < 1) exit + if (a%ia1(i1-1) /= a%ia1(i1)) exit + i1 = i1 - 1 + end do + ilr = ir + end if + nc = i2-i1+1 + call issrch(ip,ic,nc,a%ia2(i1:i2)) + if (ip>0) then + a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) + else + info = i + return + end if + end if + end do + + case default + info = -3 + write(0,*) 'Duplicate handling: ',dupl + end select + + else + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + do i=1, nz + ir = ia(i) + ic = ja(i) if (ir /= ilr) then call ibsrch(i1,ir,nnz,a%ia1) i2 = i1 @@ -279,16 +499,13 @@ contains info = i return end if - end if - end do - case(psb_dupl_add_) - ! Add - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - ic = gtl(ic) + end do + + case(psb_dupl_add_) + ! Add + do i=1, nz + ir = ia(i) + ic = ja(i) if (ir /= ilr) then call ibsrch(i1,ir,nnz,a%ia1) i2 = i1 @@ -312,28 +529,30 @@ contains info = i return end if - end if - end do + end do - case default - info = -3 - write(0,*) 'Duplicate handling: ',dupl - end select + case default + info = -3 + write(0,*) 'Duplicate handling: ',dupl + end select - end subroutine d_coo_srch_upd + end if + end subroutine d_coo_srch_upd - subroutine z_csr_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,& - & imin,imax,jmin,jmax,nzl,info) + subroutine z_csr_srch_upd(nz,ia,ja,val,nza,a,& + & imin,imax,jmin,jmax,nzl,info,gtl,ng) implicit none type(psb_zspmat_type), intent(inout) :: a - integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng - integer, intent(in) :: ia(*),ja(*),gtl(*) + integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl + integer, intent(in) :: ia(*),ja(*) integer, intent(inout) :: nza complex(kind(1.d0)), intent(in) :: val(*) integer, intent(out) :: info + integer, intent(in), optional :: ng,gtl(*) + integer :: i,ir,ic,check_flag, ilr, ilc, ip, & & i1,i2,nc,lb,ub,m,nnz,dupl @@ -341,19 +560,112 @@ contains dupl = psb_sp_getifld(psb_dupl_,a,info) - select case(dupl) - case(psb_dupl_ovwrt_,psb_dupl_err_) - ! Overwrite. - ! Cannot test for error, should have been caught earlier. - - ilr = -1 - ilc = -1 - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - ic = gtl(ic) + + if (present(gtl)) then + if (.not.present(ng)) then + info = -1 + return + endif + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + i1 = a%ia2(ir) + i2 = a%ia2(ir+1) + nc=i2-i1 + + + if (.true.) then + call issrch(ip,ic,nc,a%ia1(i1:i2-1)) + if (ip>0) then + a%aspk(i1+ip-1) = val(i) + else + write(0,*)'Was searching ',ic,' in: ',i1,i2,' : ',a%ia1(i1:i2-1) + info = i + return + end if + + else +!!$ + ip = -1 + lb = i1 + ub = i2-1 + do + if (lb > ub) exit + m = (lb+ub)/2 +!!$ write(0,*) 'Debug: ',lb,m,ub + if (ic == a%ia1(m)) then + ip = m + lb = ub + 1 + else if (ic < a%ia1(m)) then + ub = m-1 + else + lb = m + 1 + end if + enddo + + if (ip>0) then + a%aspk(ip) = val(i) + else + write(0,*)'Was searching ',ic,' in: ',i1,i2,' : ',a%ia1(i1:i2-1) + info = i + return + end if + + end if + end if + end do + + case(psb_dupl_add_) + ! Add + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + i1 = a%ia2(ir) + i2 = a%ia2(ir+1) + nc = i2-i1 + call issrch(ip,ic,nc,a%ia1(i1:i2-1)) + if (ip>0) then + a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) + else + info = i + return + end if + end if + end do + + case default + info = -3 + write(0,*) 'Duplicate handling: ',dupl + end select + + else + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) i1 = a%ia2(ir) i2 = a%ia2(ir+1) nc=i2-i1 @@ -397,19 +709,15 @@ contains end if end if - end if - end do - - case(psb_dupl_add_) - ! Add - ilr = -1 - ilc = -1 - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - ic = gtl(ic) + end do + + case(psb_dupl_add_) + ! Add + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) i1 = a%ia2(ir) i2 = a%ia2(ir+1) nc = i2-i1 @@ -420,25 +728,27 @@ contains info = i return end if - end if - end do + end do + + case default + info = -3 + write(0,*) 'Duplicate handling: ',dupl + end select + end if - case default - info = -3 - write(0,*) 'Duplicate handling: ',dupl - end select end subroutine z_csr_srch_upd - subroutine z_coo_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,& - & imin,imax,jmin,jmax,nzl,info) + subroutine z_coo_srch_upd(nz,ia,ja,val,nza,a,& + & imin,imax,jmin,jmax,nzl,info,gtl,ng) implicit none type(psb_zspmat_type), intent(inout) :: a - integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng - integer, intent(in) :: ia(*),ja(*),gtl(*) + integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl + integer, intent(in) :: ia(*),ja(*) integer, intent(inout) :: nza complex(kind(1.d0)), intent(in) :: val(*) integer, intent(out) :: info + integer, intent(in), optional :: ng,gtl(*) integer :: i,ir,ic,check_flag, ilr, ilc, ip, & & i1,i2,nc,lb,ub,m,nnz,dupl,isrt @@ -456,16 +766,95 @@ contains nnz = psb_sp_getifld(psb_nnz_,a,info) - select case(dupl) - case(psb_dupl_ovwrt_,psb_dupl_err_) - ! Overwrite. - ! Cannot test for error, should have been caught earlier. - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - ic = gtl(ic) + if (present(gtl)) then + if (.not.present(ng)) then + info = -1 + return + endif + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if (ir /= ilr) then + call ibsrch(i1,ir,nnz,a%ia1) + i2 = i1 + do + if (i2+1 > nnz) exit + if (a%ia1(i2+1) /= a%ia1(i2)) exit + i2 = i2 + 1 + end do + do + if (i1-1 < 1) exit + if (a%ia1(i1-1) /= a%ia1(i1)) exit + i1 = i1 - 1 + end do + ilr = ir + end if + nc = i2-i1+1 + call issrch(ip,ic,nc,a%ia2(i1:i2)) + if (ip>0) then + a%aspk(i1+ip-1) = val(i) + else + info = i + return + end if + end if + end do + case(psb_dupl_add_) + ! Add + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if (ir /= ilr) then + call ibsrch(i1,ir,nnz,a%ia1) + i2 = i1 + do + if (i2+1 > nnz) exit + if (a%ia1(i2+1) /= a%ia1(i2)) exit + i2 = i2 + 1 + end do + do + if (i1-1 < 1) exit + if (a%ia1(i1-1) /= a%ia1(i1)) exit + i1 = i1 - 1 + end do + ilr = ir + end if + nc = i2-i1+1 + call issrch(ip,ic,nc,a%ia2(i1:i2)) + if (ip>0) then + a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) + else + info = i + return + end if + end if + end do + + case default + info = -3 + write(0,*) 'Duplicate handling: ',dupl + end select + + else + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + do i=1, nz + ir = ia(i) + ic = ja(i) if (ir /= ilr) then call ibsrch(i1,ir,nnz,a%ia1) i2 = i1 @@ -476,7 +865,6 @@ contains end do do if (i1-1 < 1) exit - if (a%ia1(i1-1) /= a%ia1(i1)) exit i1 = i1 - 1 end do @@ -487,19 +875,16 @@ contains if (ip>0) then a%aspk(i1+ip-1) = val(i) else - info = i + info = i return end if - end if - end do - case(psb_dupl_add_) - ! Add - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - ic = gtl(ic) + end do + + case(psb_dupl_add_) + ! Add + do i=1, nz + ir = ia(i) + ic = ja(i) if (ir /= ilr) then call ibsrch(i1,ir,nnz,a%ia1) i2 = i1 @@ -520,19 +905,19 @@ contains if (ip>0) then a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) else - info = i + info = i return end if - end if - end do + end do - case default - info = -3 - write(0,*) 'Duplicate handling: ',dupl - end select + case default + info = -3 + write(0,*) 'Duplicate handling: ',dupl + end select - end subroutine z_coo_srch_upd + end if + end subroutine z_coo_srch_upd end module psb_update_mod diff --git a/src/serial/psb_zcoins.f90 b/src/serial/psb_zcoins.f90 index 86943164..959f64dd 100644 --- a/src/serial/psb_zcoins.f90 +++ b/src/serial/psb_zcoins.f90 @@ -31,7 +31,7 @@ ! File: psbzcoins.f90 ! Subroutine: ! Parameters: -subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild) +subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild) use psb_spmat_type use psb_const_mod @@ -43,11 +43,12 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild) implicit none integer, intent(in) :: nz, imin,imax,jmin,jmax - integer, intent(in) :: ia(:),ja(:),gtl(:) + integer, intent(in) :: ia(:),ja(:) complex(kind(1.d0)), intent(in) :: val(:) type(psb_zspmat_type), intent(inout) :: a integer, intent(out) :: info - logical, optional, intent(in) :: rebuild + integer, intent(in), optional :: gtl(:) + logical, intent(in), optional :: rebuild character(len=5) :: ufida integer :: i,j,ir,ic,nr,nc, ng, nza, isza,spstate, nnz,& @@ -95,153 +96,312 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild) end if call touppers(a%fida,ufida) - ng = size(gtl) spstate = psb_sp_getifld(psb_state_,a,info) - select case(spstate) - case(psb_spmat_bld_) - if ((ufida /= 'COO').and.(ufida/='COI')) then - info = 134 - ch_err(1:3)=ufida(1:3) - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - call psb_spinfo(psb_nztotreq_,a,nza,info) - call psb_spinfo(psb_nzsizereq_,a,isza,info) - if(info /= izero) then - info=4010 - ch_err='psb_spinfo' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - endif + if (present(gtl)) then - if ((nza+nz)>isza) then - call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info) + ng = size(gtl) + select case(spstate) + case(psb_spmat_bld_) + if ((ufida /= 'COO').and.(ufida/='COI')) then + info = 134 + ch_err(1:3)=ufida(1:3) + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_spinfo(psb_nztotreq_,a,nza,info) + call psb_spinfo(psb_nzsizereq_,a,isza,info) if(info /= izero) then info=4010 - ch_err='psb_sp_reall' + ch_err='psb_spinfo' call psb_errpush(info,name,a_err=ch_err) goto 9999 endif - endif - call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,gtl,ng,& - & imin,imax,jmin,jmax,info) - if(info /= izero) then - info=4010 - ch_err='psb_inner_ins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - endif - if (debug) then - if ((nza - a%infoa(psb_nnz_)) /= nz) then - write(0,*) 'PSB_COINS: insert discarded items ' + + if ((nza+nz)>isza) then + call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info) + if(info /= izero) then + info=4010 + ch_err='psb_sp_reall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + endif + + call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,& + & imin,imax,jmin,jmax,info,gtl,ng) + + if(info /= izero) then + info=4010 + ch_err='psb_inner_ins' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + if (debug) then + if ((nza - a%infoa(psb_nnz_)) /= nz) then + write(0,*) 'PSB_COINS: insert discarded items ' + end if end if - end if - if ((nza - psb_sp_getifld(psb_nnz_,a,info)) /= nz) then - call psb_sp_setifld(nza,psb_del_bnd_,a,info) - endif - call psb_sp_setifld(nza,psb_nnz_,a,info) + if ((nza - psb_sp_getifld(psb_nnz_,a,info)) /= nz) then + call psb_sp_setifld(nza,psb_del_bnd_,a,info) + endif + call psb_sp_setifld(nza,psb_nnz_,a,info) + + case(psb_spmat_upd_) + + iupd = psb_sp_getifld(psb_upd_,a,info) + select case (iupd) + case (psb_upd_perm_) + ip1 = psb_sp_getifld(psb_upd_pnt_,a,info) + nzl = psb_sp_getifld(psb_del_bnd_,a,info) + nza = a%ia2(ip1+psb_nnz_) + + + call psb_inner_upd(nz,ia,ja,val,nza,a%aspk,& + & imin,imax,jmin,jmax,nzl,info,gtl,ng) + + if (info /= izero) then + info=4010 + ch_err='psb_inner_upd' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + if (debug) then + if ((nza - a%ia2(ip1+psb_nnz_)) /= nz) then + write(0,*) 'PSB_COINS: update discarded items ' + end if + end if + a%ia2(ip1+psb_nnz_) = nza + + if (debug) write(0,*) 'From COINS(UPD) : NZA:',nza + + case (psb_upd_dflt_, psb_upd_srch_) + + call psb_srch_upd(nz,ia,ja,val,nza,a,& + & imin,imax,jmin,jmax,nzl,info,gtl,ng) + + if (info > 0) then + if (rebuild_) then + if (debug) write(0,*)& + & 'COINS: Going through rebuild_ fingers crossed!' + irst = info + call psb_nullify_sp(tmp) + tmp%fida='COO' + call psb_csdp(a,tmp,info) + call psb_sp_setifld(psb_spmat_bld_,psb_state_,tmp,info) + if (debug) then + write(0,*) 'COINS Rebuild: size',tmp%infoa(psb_nnz_) ,irst + endif + call psb_sp_transfer(tmp,a,info) + call psb_spinfo(psb_nztotreq_,a,nza,info) + call psb_spinfo(psb_nzsizereq_,a,isza,info) + if(info /= izero) then + info=4010 + ch_err='psb_spinfo' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + + if (debug) write(0,*)& + & 'COINS: Reinserting',a%fida,nza,isza + if ((nza+nz)>isza) then + call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info) + if(info /= izero) then + info=4010 + ch_err='psb_sp_reall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + endif + if (irst <= nz) then + call psb_inner_ins((nz-irst+1),ia(irst:nz),ja(irst:nz),val(irst:nz),& + & nza,a%ia1,a%ia2,a%aspk,imin,imax,jmin,jmax,info,gtl,ng) + call psb_sp_setifld(nza,psb_del_bnd_,a,info) + call psb_sp_setifld(nza,psb_nnz_,a,info) + end if + + else + info = 2231 + call psb_errpush(info,name) + goto 9999 + end if + else if (info < 0) then + info = 2230 + call psb_errpush(info,name) + goto 9999 + + end if + + case default + + info = 2231 + call psb_errpush(info,name) + goto 9999 + end select + + + case default + info = 2232 + call psb_errpush(info,name) + goto 9999 + end select + - case(psb_spmat_upd_) - iupd = psb_sp_getifld(psb_upd_,a,info) - select case (iupd) - case (psb_upd_perm_) - ip1 = psb_sp_getifld(psb_upd_pnt_,a,info) - nzl = psb_sp_getifld(psb_del_bnd_,a,info) - nza = a%ia2(ip1+psb_nnz_) - nza = a%ia2(ip1+psb_nnz_) - nzl = a%infoa(psb_del_bnd_) - call psb_inner_upd(nz,ia,ja,val,nza,a%aspk,gtl,ng,& - & imin,imax,jmin,jmax,nzl,info) + else + + ng = -1 - if (info /= izero) then + select case(spstate) + case(psb_spmat_bld_) + if ((ufida /= 'COO').and.(ufida/='COI')) then + info = 134 + ch_err(1:3)=ufida(1:3) + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_spinfo(psb_nztotreq_,a,nza,info) + call psb_spinfo(psb_nzsizereq_,a,isza,info) + if(info /= izero) then + info=4010 + ch_err='psb_spinfo' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + + if ((nza+nz)>isza) then + call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info) + if(info /= izero) then + info=4010 + ch_err='psb_sp_reall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + endif + + call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,& + & imin,imax,jmin,jmax,info) + + if(info /= izero) then info=4010 - ch_err='psb_inner_upd' + ch_err='psb_inner_ins' call psb_errpush(info,name,a_err=ch_err) goto 9999 endif if (debug) then - if ((nza - a%ia2(ip1+psb_nnz_)) /= nz) then - write(0,*) 'PSB_COINS: update discarded items ' + if ((nza - a%infoa(psb_nnz_)) /= nz) then + write(0,*) 'PSB_COINS: insert discarded items ' end if end if - a%ia2(ip1+psb_nnz_) = nza + if ((nza - psb_sp_getifld(psb_nnz_,a,info)) /= nz) then + call psb_sp_setifld(nza,psb_del_bnd_,a,info) + endif + call psb_sp_setifld(nza,psb_nnz_,a,info) + + case(psb_spmat_upd_) - if (debug) write(0,*) 'From COINS(UPD) : NZA:',nza + iupd = psb_sp_getifld(psb_upd_,a,info) + select case (iupd) + case (psb_upd_perm_) + ip1 = psb_sp_getifld(psb_upd_pnt_,a,info) + nzl = psb_sp_getifld(psb_del_bnd_,a,info) + nza = a%ia2(ip1+psb_nnz_) - case (psb_upd_dflt_, psb_upd_srch_) - call psb_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,& + call psb_inner_upd(nz,ia,ja,val,nza,a%aspk,& & imin,imax,jmin,jmax,nzl,info) - if (info > 0) then - if (rebuild_) then - if (debug) write(0,*)& - & 'COINS: Going through rebuild_ fingers crossed!' - irst = info - call psb_nullify_sp(tmp) - tmp%fida='COO' - call psb_csdp(a,tmp,info) - call psb_sp_setifld(psb_spmat_bld_,psb_state_,tmp,info) - if (debug) then - write(0,*) 'COINS Rebuild: size',tmp%infoa(psb_nnz_) ,irst - endif - call psb_sp_transfer(tmp,a,info) - call psb_spinfo(psb_nztotreq_,a,nza,info) - call psb_spinfo(psb_nzsizereq_,a,isza,info) - if(info /= izero) then - info=4010 - ch_err='psb_spinfo' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - endif - - if (debug) write(0,*)& - & 'COINS: Reinserting',a%fida,nza,isza - if ((nza+nz)>isza) then - call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info) + + if (info /= izero) then + info=4010 + ch_err='psb_inner_upd' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + if (debug) then + if ((nza - a%ia2(ip1+psb_nnz_)) /= nz) then + write(0,*) 'PSB_COINS: update discarded items ' + end if + end if + a%ia2(ip1+psb_nnz_) = nza + + if (debug) write(0,*) 'From COINS(UPD) : NZA:',nza + + case (psb_upd_dflt_, psb_upd_srch_) + + call psb_srch_upd(nz,ia,ja,val,nza,a,& + & imin,imax,jmin,jmax,nzl,info) + + if (info > 0) then + if (rebuild_) then + if (debug) write(0,*)& + & 'COINS: Going through rebuild_ fingers crossed!' + irst = info + call psb_nullify_sp(tmp) + tmp%fida='COO' + call psb_csdp(a,tmp,info) + call psb_sp_setifld(psb_spmat_bld_,psb_state_,tmp,info) + if (debug) then + write(0,*) 'COINS Rebuild: size',tmp%infoa(psb_nnz_) ,irst + endif + call psb_sp_transfer(tmp,a,info) + call psb_spinfo(psb_nztotreq_,a,nza,info) + call psb_spinfo(psb_nzsizereq_,a,isza,info) if(info /= izero) then info=4010 - ch_err='psb_sp_reall' + ch_err='psb_spinfo' call psb_errpush(info,name,a_err=ch_err) goto 9999 endif - endif - if (irst <= nz) then - call psb_inner_ins((nz-irst+1),ia(irst:nz),ja(irst:nz),val(irst:nz),& - & nza,a%ia1,a%ia2,a%aspk,gtl,ng,imin,imax,jmin,jmax,info) - call psb_sp_setifld(nza,psb_del_bnd_,a,info) - call psb_sp_setifld(nza,psb_nnz_,a,info) + + if (debug) write(0,*)& + & 'COINS: Reinserting',a%fida,nza,isza + if ((nza+nz)>isza) then + call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info) + if(info /= izero) then + info=4010 + ch_err='psb_sp_reall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + endif + if (irst <= nz) then + call psb_inner_ins((nz-irst+1),ia(irst:nz),ja(irst:nz),val(irst:nz),& + & nza,a%ia1,a%ia2,a%aspk,imin,imax,jmin,jmax,info) + call psb_sp_setifld(nza,psb_del_bnd_,a,info) + call psb_sp_setifld(nza,psb_nnz_,a,info) + end if + + else + info = 2231 + call psb_errpush(info,name) + goto 9999 end if - - else - info = 2231 + else if (info < 0) then + info = 2230 call psb_errpush(info,name) goto 9999 + end if - else if (info < 0) then - info = 2230 + + case default + + info = 2231 call psb_errpush(info,name) goto 9999 + end select - end if - - case default - info = 2231 + case default + info = 2232 call psb_errpush(info,name) - goto 9999 + goto 9999 end select + endif - case default - info = 2232 - call psb_errpush(info,name) - goto 9999 - end select return call psb_erractionrestore(err_act) @@ -257,69 +417,108 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild) contains - subroutine psb_inner_upd(nz,ia,ja,val,nza,aspk,gtl,ng,& - & imin,imax,jmin,jmax,nzl,info) + subroutine psb_inner_upd(nz,ia,ja,val,nza,aspk,& + & imin,imax,jmin,jmax,nzl,info,gtl,ng) implicit none - integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng - integer, intent(in) :: ia(*),ja(*),gtl(*) + integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl + integer, intent(in) :: ia(*),ja(*) integer, intent(inout) :: nza complex(kind(1.d0)), intent(in) :: val(*) complex(kind(1.d0)), intent(inout) :: aspk(*) integer, intent(out) :: info + integer, intent(in), optional :: ng,gtl(*) integer :: i,ir,ic - if (nza >= nzl) then - do i=1, nz - nza = nza + 1 - aspk(nza) = val(i) - end do + if (present(gtl)) then + if (.not.present(ng)) then + info = -1 + return + endif + if (nza >= nzl) then + do i=1, nz + nza = nza + 1 + aspk(nza) = val(i) + end do + else + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then + nza = nza + 1 + aspk(nza) = val(i) + end if + end if + end do + end if else - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - ic = gtl(ic) + if (nza >= nzl) then + do i=1, nz + nza = nza + 1 + aspk(nza) = val(i) + end do + else + do i=1, nz + ir = ia(i) + ic = ja(i) if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then nza = nza + 1 aspk(nza) = val(i) end if - end if - end do + end do + end if end if - end subroutine psb_inner_upd - subroutine psb_inner_ins(nz,ia,ja,val,nza,ia1,ia2,aspk,gtl,ng,& - & imin,imax,jmin,jmax,info) + subroutine psb_inner_ins(nz,ia,ja,val,nza,ia1,ia2,aspk,& + & imin,imax,jmin,jmax,info,gtl,ng) implicit none - integer, intent(in) :: nz, imin,imax,jmin,jmax,ng - integer, intent(in) :: ia(*),ja(*),gtl(*) + integer, intent(in) :: nz, imin,imax,jmin,jmax + integer, intent(in) :: ia(*),ja(*) integer, intent(inout) :: nza,ia1(*),ia2(*) complex(kind(1.d0)), intent(in) :: val(*) complex(kind(1.d0)), intent(inout) :: aspk(*) integer, intent(out) :: info - + integer, intent(in), optional :: ng,gtl(*) integer :: i,ir,ic info = 0 - do i=1, nz - ir = ia(i) - ic = ja(i) + if (present(gtl)) then + if (.not.present(ng)) then + info = -1 + return + endif + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then + nza = nza + 1 + ia1(nza) = ir + ia2(nza) = ic + aspk(nza) = val(i) + end if + end if + end do + else - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - ic = gtl(ic) + do i=1, nz + ir = ia(i) + ic = ja(i) if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then nza = nza + 1 ia1(nza) = ir ia2(nza) = ic aspk(nza) = val(i) end if - end if - end do + end do + end if end subroutine psb_inner_ins diff --git a/src/tools/psb_dspins.f90 b/src/tools/psb_dspins.f90 index 682bec37..8b51a11c 100644 --- a/src/tools/psb_dspins.f90 +++ b/src/tools/psb_dspins.f90 @@ -145,7 +145,7 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild) ncol = desc_a%matrix_data(psb_n_col_) if (spstate == psb_spmat_bld_) then - call psb_coins(nz,ia,ja,val,a,desc_a%glob_to_loc,1,nrow,1,ncol,info) + call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,info,gtl=desc_a%glob_to_loc) if (info /= 0) then info=4010 ch_err='psb_coins' @@ -160,8 +160,8 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild) else if (psb_is_asb_dec(dectype)) then nrow = desc_a%matrix_data(psb_n_row_) ncol = desc_a%matrix_data(psb_n_col_) - call psb_coins(nz,ia,ja,val,a,desc_a%glob_to_loc,1,nrow,1,ncol,& - & info,rebuild=rebuild_) + call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,& + & info,gtl=desc_a%glob_to_loc,rebuild=rebuild_) if (info /= 0) then info=4010 ch_err='psb_coins' diff --git a/src/tools/psb_zspins.f90 b/src/tools/psb_zspins.f90 index 9e79edd9..e11fd3ae 100644 --- a/src/tools/psb_zspins.f90 +++ b/src/tools/psb_zspins.f90 @@ -145,7 +145,7 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild) ncol = desc_a%matrix_data(psb_n_col_) if (spstate == psb_spmat_bld_) then - call psb_coins(nz,ia,ja,val,a,desc_a%glob_to_loc,1,nrow,1,ncol,info) + call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,info,gtl=desc_a%glob_to_loc) if (info /= 0) then info=4010 ch_err='psb_coins' @@ -160,8 +160,8 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild) else if (psb_is_asb_dec(dectype)) then nrow = desc_a%matrix_data(psb_n_row_) ncol = desc_a%matrix_data(psb_n_col_) - call psb_coins(nz,ia,ja,val,a,desc_a%glob_to_loc,1,nrow,1,ncol,& - & info,rebuild=rebuild_) + call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,& + & info,gtl=desc_a%glob_to_loc,rebuild=rebuild_) if (info /= 0) then info=4010 ch_err='psb_coins'