Added JAD storage format specific code

psblas3-type-indexed
Alfredo Buttari 19 years ago
parent 4271226c8d
commit fcd0fdebc1

@ -16,7 +16,21 @@ subroutine psb_dneigh(a,idx,neigh,n,info,lev)
integer, pointer :: neigh(:) ! the neighbours integer, pointer :: neigh(:) ! the neighbours
integer, optional :: lev ! level of neighbours to find integer, optional :: lev ! level of neighbours to find
interface psb_spgtrow
subroutine psb_dspgtrow(irw,a,b,info,append,iren,lrw)
use psb_spmat_type
type(psb_dspmat_type), intent(in) :: a
integer, intent(in) :: irw
type(psb_dspmat_type), intent(inout) :: b
logical, intent(in), optional :: append
integer, intent(in), target, optional :: iren(:)
integer, intent(in), optional :: lrw
integer, intent(out) :: info
end subroutine psb_dspgtrow
end interface
type(psb_dspmat_type) :: atmp
integer, pointer :: tmpn(:)
integer :: level, dim, i, j, k, r, c, brow,& integer :: level, dim, i, j, k, r, c, brow,&
& elem_pt, ii, n1, col_idx, ne, err_act & elem_pt, ii, n1, col_idx, ne, err_act
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -25,12 +39,6 @@ subroutine psb_dneigh(a,idx,neigh,n,info,lev)
info = 0 info = 0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if ((a%fida /= 'CSR')) then
info=135
ch_err=a%fida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
n = 0 n = 0
info = 0 info = 0
if(present(lev)) then if(present(lev)) then
@ -44,7 +52,51 @@ subroutine psb_dneigh(a,idx,neigh,n,info,lev)
level=1 level=1
end if end if
if(a%fida.eq.'CSR') then if ((a%fida /= 'CSR')) then
call psb_spall(atmp,1,info)
if(info.ne.0) then
call psb_errpush(4010,name)
goto 9999
end if
call psb_spgtrow(idx,a,atmp,info)
if(info.ne.0) then
call psb_errpush(4010,name)
goto 9999
end if
dim=atmp%infoa(psb_nnz_)
allocate(tmpn(dim))
tmpn(1:dim)=atmp%ia2(1:dim)
if(level.eq.2) then
do i=1,dim
if ((1<=tmpn(i)).and.(tmpn(i)<=a%m).and.(tmpn(i).ne.idx)) then
call psb_spgtrow(tmpn(i),a,atmp,info,append=.true.)
if(info.ne.0) then
call psb_errpush(4010,name)
goto 9999
end if
end if
end do
end if
dim=atmp%infoa(psb_nnz_)
if(dim.gt.size(neigh)) &
& call psb_realloc(dim,neigh,info)
if(info.ne.0) then
call psb_errpush(4010,name)
goto 9999
end if
call psb_spfree(atmp,info)
if(info.ne.0) then
call psb_errpush(4010,name)
goto 9999
end if
call psb_nullify_sp(atmp)
deallocate(tmpn)
else if(a%fida.eq.'CSR') then
dim=0 dim=0
if(level.eq.1) dim=(a%ia2(idx+1)-a%ia2(idx)) if(level.eq.1) dim=(a%ia2(idx+1)-a%ia2(idx))

@ -21,7 +21,21 @@ subroutine psb_dspgtdiag(a,d,info)
real(kind(1.d0)), intent(inout) :: d(:) real(kind(1.d0)), intent(inout) :: d(:)
integer, intent(out) :: info integer, intent(out) :: info
integer :: i,j,k,nr, nz, err_act interface psb_spgtrow
subroutine psb_dspgtrow(irw,a,b,info,append,iren,lrw)
use psb_spmat_type
type(psb_dspmat_type), intent(in) :: a
integer, intent(in) :: irw
type(psb_dspmat_type), intent(inout) :: b
logical, intent(in), optional :: append
integer, intent(in), target, optional :: iren(:)
integer, intent(in), optional :: lrw
integer, intent(out) :: info
end subroutine psb_dspgtrow
end interface
type(psb_dspmat_type) :: tmpa
integer :: i,j,k,nr, nz, err_act, ii, rng, irb, nrb
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_dspgtdiag' name='psb_dspgtdiag'
@ -52,12 +66,29 @@ subroutine psb_dspgtdiag(a,d,info)
enddo enddo
else if (a%fida == 'JAD') then else if (a%fida == 'JAD') then
info=135
ch_err=a%fida(1:3) rng=min(a%m,a%k)
do i=1, rng, nrb
irb=min(i+nrb-1,rng)
call psb_spgtrow(i,a,tmpa,info,lrw=irb)
if(info.ne.0) then
info=4010
ch_err='psb_spgtrow'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
do ii=1,tmpa%infoa(psb_nnz_)
j=tmpa%ia1(ii)
if ((j==tmpa%ia2(ii)).and.(j <= rng) .and.(j>0)) then
d(j) = tmpa%aspk(ii)
endif
enddo
end do
end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

@ -232,14 +232,13 @@ contains
end if end if
if (a%infoa(psb_srtd_) == psb_isrtdcoo_) then if (a%infoa(psb_srtd_) == psb_isrtdcoo_) then
!!$ write(0,*) 'Gtrow_: srtd coo',irw
! In this case we can do a binary search. ! In this case we can do a binary search.
do do
call ibsrch(ip,irw,nza,a%ia1) call ibsrch(ip,irw,nza,a%ia1)
if (ip /= -1) exit if (ip /= -1) exit
irw = irw + 1 irw = irw + 1
if (irw > lrw) then if (irw > lrw) then
write(0,*) 'Warning : did not find any rows. Is this an error?' write(0,*) 'Warning : did not find any rows. Is this an error? ',irw,lrw,idx
exit exit
end if end if
end do end do

Loading…
Cancel
Save