fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

*** empty log message ***

Browse Source
psblas3-type-indexed
Salvatore Filippone 19 years ago
parent 00da2086a5
commit cbc36fdd48
6 changed files with 1244 additions and 456 deletions
Show Stats Download Patch File Download Diff File

10
src/modules/psb_serial_mod.f90
Unescape Escape View File

@ -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

240
src/serial/psb_dcoins.f90
Unescape Escape View File

@ -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,11 +96,166 @@ 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)
if (present(gtl)) then
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_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,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
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
else
ng = -1
select case(spstate)
case(psb_spmat_bld_)
if ((ufida /= 'COO').and.(ufida/='COI')) then
info = 134
@ -125,8 +281,10 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild)
goto 9999
endif
endif
call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,gtl,ng,&
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'
@ -152,10 +310,8 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild)
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,&
call psb_inner_upd(nz,ia,ja,val,nza,a%aspk,&
& imin,imax,jmin,jmax,nzl,info)
if (info /= izero) then
@ -175,8 +331,9 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild)
case (psb_upd_dflt_, psb_upd_srch_)
call psb_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,&
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,*)&
@ -212,7 +369,7 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild)
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)
& 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
@ -242,6 +399,9 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild)
call psb_errpush(info,name)
goto 9999
end select
endif
return
call psb_erractionrestore(err_act)
@ -257,18 +417,24 @@ 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 (present(gtl)) then
if (.not.present(ng)) then
info = -1
return
endif
if (nza >= nzl) then
do i=1, nz
nza = nza + 1
@ -288,26 +454,47 @@ contains
end if
end do
end if
else
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 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
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)
@ -319,8 +506,23 @@ contains
end if
end if
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
ia1(nza) = ir
ia2(nza) = ic
aspk(nza) = val(i)
end if
end do
end if
end subroutine psb_inner_ins
end subroutine psb_dcoins

447
src/serial/psb_update_mod.f90
Unescape Escape View File

@ -51,27 +51,50 @@ 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
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,&
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,&
& 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,gtl,ng,&
call coo_srch_upd(nz,ia,ja,val,nza,a,&
& imin,imax,jmin,jmax,nzl,info)
case default
@ -79,30 +102,54 @@ contains
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
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,&
& 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,&
& 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,gtl,ng,&
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,gtl,ng,&
call coo_srch_upd(nz,ia,ja,val,nza,a,&
& imin,imax,jmin,jmax,nzl,info)
case default
@ -110,20 +157,23 @@ contains
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,6 +181,13 @@ contains
dupl = psb_sp_getifld(psb_dupl_,a,info)
if (present(gtl)) then
if (.not.present(ng)) then
info = -1
return
endif
select case(dupl)
case(psb_dupl_ovwrt_,psb_dupl_err_)
! Overwrite.
@ -217,18 +274,102 @@ contains
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
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 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
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 do
case default
info = -3
write(0,*) 'Duplicate handling: ',dupl
end select
end if
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,6 +387,12 @@ contains
nnz = psb_sp_getifld(psb_nnz_,a,info)
if (present(gtl)) then
if (.not.present(ng)) then
info = -1
return
endif
select case(dupl)
case(psb_dupl_ovwrt_,psb_dupl_err_)
! Overwrite.
@ -320,20 +467,92 @@ contains
write(0,*) 'Duplicate handling: ',dupl
end select
end subroutine d_coo_srch_upd
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
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 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
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 do
subroutine z_csr_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
case default
info = -3
write(0,*) 'Duplicate handling: ',dupl
end select
end if
end subroutine d_coo_srch_upd
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,6 +560,13 @@ contains
dupl = psb_sp_getifld(psb_dupl_,a,info)
if (present(gtl)) then
if (.not.present(ng)) then
info = -1
return
endif
select case(dupl)
case(psb_dupl_ovwrt_,psb_dupl_err_)
! Overwrite.
@ -427,18 +653,102 @@ contains
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
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 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
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 do
case default
info = -3
write(0,*) 'Duplicate handling: ',dupl
end select
end if
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,6 +766,12 @@ contains
nnz = psb_sp_getifld(psb_nnz_,a,info)
if (present(gtl)) then
if (.not.present(ng)) then
info = -1
return
endif
select case(dupl)
case(psb_dupl_ovwrt_,psb_dupl_err_)
! Overwrite.
@ -476,7 +792,6 @@ contains
end do
do
if (i1-1 < 1) exit
if (a%ia1(i1-1) /= a%ia1(i1)) exit
i1 = i1 - 1
end do
@ -531,8 +846,78 @@ contains
write(0,*) 'Duplicate handling: ',dupl
end select
end subroutine z_coo_srch_upd
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
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 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
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 do
case default
info = -3
write(0,*) 'Duplicate handling: ',dupl
end select
end if
end subroutine z_coo_srch_upd
end module psb_update_mod

241
src/serial/psb_zcoins.f90
Unescape Escape View File

@ -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,11 +96,166 @@ 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)
if (present(gtl)) then
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_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,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
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
else
ng = -1
select case(spstate)
case(psb_spmat_bld_)
if ((ufida /= 'COO').and.(ufida/='COI')) then
info = 134
@ -125,8 +281,10 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild)
goto 9999
endif
endif
call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,gtl,ng,&
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'
@ -152,10 +310,8 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild)
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,&
call psb_inner_upd(nz,ia,ja,val,nza,a%aspk,&
& imin,imax,jmin,jmax,nzl,info)
if (info /= izero) then
@ -175,8 +331,9 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild)
case (psb_upd_dflt_, psb_upd_srch_)
call psb_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,&
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,*)&
@ -212,7 +369,7 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild)
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)
& 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
@ -242,6 +399,9 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info,rebuild)
call psb_errpush(info,name)
goto 9999
end select
endif
return
call psb_erractionrestore(err_act)
@ -257,18 +417,24 @@ 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 (present(gtl)) then
if (.not.present(ng)) then
info = -1
return
endif
if (nza >= nzl) then
do i=1, nz
nza = nza + 1
@ -288,27 +454,47 @@ contains
end if
end do
end if
else
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 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
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)
@ -320,6 +506,19 @@ contains
end if
end if
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
ia1(nza) = ir
ia2(nza) = ic
aspk(nza) = val(i)
end if
end do
end if
end subroutine psb_inner_ins

6
src/tools/psb_dspins.f90
Unescape Escape View File

@ -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'

6
src/tools/psb_zspins.f90
Unescape Escape View File

@ -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'

Write Preview
Loading…
Cancel
Save