mirror of
https://github.com/sfilippone/psblas3.git
synced 2026-10-20 21:02:26 +00:00
Changed implementation and internals of GETROW/GETBLK.
This commit is contained in:
Salvatore Filippone
@@ -12,7 +12,7 @@ FOBJS = psb_cest.o psb_dcoins.o psb_dcsdp.o psb_dcsmm.o psb_dcsmv.o \
|
||||
psb_zfixcoo.o psb_zipcoo2csr.o psb_zipcsr2coo.o psb_zipcoo2csc.o \
|
||||
psb_zcoins.o psb_zcsprt.o psb_zneigh.o psb_ztransp.o psb_ztransc.o\
|
||||
psb_zrwextd.o psb_zsymbmm.o psb_znumbmm.o psb_zspscal.o psb_zspclip.o\
|
||||
psb_getifield.o psb_setifield.o psb_update_mod.o
|
||||
psb_getifield.o psb_setifield.o psb_update_mod.o psb_getrow_mod.o
|
||||
|
||||
LIBDIR=..
|
||||
MODDIR=../modules
|
||||
@@ -26,6 +26,7 @@ lib: auxd cood csrd jadd f77d dpd lib1
|
||||
lib1: $(FOBJS)
|
||||
|
||||
psb_dcoins.o psb_zcoins.o: psb_update_mod.o
|
||||
psb_dspgetrow.o: psb_getrow_mod.o
|
||||
|
||||
auxd:
|
||||
(cd aux; make lib)
|
||||
|
||||
@@ -54,11 +54,12 @@ How do you add a new storage format? Here is your checklist.
|
||||
11. You have to provide the query facilities used in PSB_XSPINFO
|
||||
(included in base/modules/psb_spmat_type).
|
||||
12. You have to provide the functionality to extract a block of
|
||||
rows. This is used in the GETBLK/GETROW/CLIP chain (TO BE
|
||||
REVISED). GETROW is (all that is) used by the ILU factorization.
|
||||
rows in psb_getrow_mod.f90. This is used in the GETBLK/GETROW/CLIP
|
||||
chain; GETROW is (all that is) used by the ILU factorization.
|
||||
13. You have to provide the GETDIAG functionality.
|
||||
14. You have to provide the functionality to NEIGH.
|
||||
WATCH OUT: THIS COULD BE CHANGED TO USE GETROW!!!!
|
||||
15. RWEXTD: what do we do here? Should we switch to/from COO?????
|
||||
15. RWEXTD: what do we do here? Should we switch to/from COO?????
|
||||
The current interface is ugly!
|
||||
|
||||
Is this complete? I sure hope so....
|
||||
@@ -74,7 +74,7 @@ subroutine psb_dipcoo2csr(a,info,rwshr)
|
||||
call psb_fixcoo(a,info)
|
||||
nr = a%m
|
||||
nza = a%infoa(psb_nnz_)
|
||||
allocate(iaux(nr+1),stat=info)
|
||||
allocate(iaux(max(nr+1,1)),stat=info)
|
||||
if (info /= 0) then
|
||||
call psb_errpush(4010,name,a_err='Allocate')
|
||||
goto 9999
|
||||
@@ -137,6 +137,7 @@ subroutine psb_dipcoo2csr(a,info,rwshr)
|
||||
if (nr < itemp(nza)) then
|
||||
write(0,*) 'IPCOO2CSR: RWSHR=.false. : ',&
|
||||
&nr,itemp(nza),' Expect trouble!'
|
||||
info = 12
|
||||
end if
|
||||
|
||||
|
||||
@@ -168,6 +169,7 @@ subroutine psb_dipcoo2csr(a,info,rwshr)
|
||||
!
|
||||
if (j /= (nza+1)) then
|
||||
write(0,*) 'IPCOO2CSR : Problem from loop :',j,nza
|
||||
info = 13
|
||||
endif
|
||||
do
|
||||
if (i>nr) exit
|
||||
|
||||
@@ -104,7 +104,7 @@ contains
|
||||
minmn = min(m,n)
|
||||
do i = 1,n
|
||||
|
||||
call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info,bw=w)
|
||||
call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info)
|
||||
|
||||
do jj=1, nazr
|
||||
j=iacl(jj)
|
||||
@@ -115,7 +115,7 @@ contains
|
||||
return
|
||||
|
||||
endif
|
||||
call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info,bw=w)
|
||||
call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info)
|
||||
do k=1,nbzr
|
||||
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
|
||||
write(0,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn
|
||||
@@ -139,8 +139,6 @@ contains
|
||||
end do
|
||||
|
||||
|
||||
|
||||
|
||||
end subroutine inner_numbmm
|
||||
|
||||
end subroutine psb_dnumbmm
|
||||
|
||||
@@ -109,7 +109,8 @@ subroutine psb_dspclip(a,b,info,imin,imax,jmin,jmax,rscale,cscale)
|
||||
|
||||
sizeb = psb_sp_get_nnzeros(a)
|
||||
call psb_sp_all(mb,kb,b,sizeb,info)
|
||||
b%fida='COO'
|
||||
b%fida = 'COO'
|
||||
b%descra = a%descra
|
||||
nzb = 0
|
||||
do i=imin_, imax_, irbk
|
||||
nrt = min(irbk,imax_-i+1)
|
||||
|
||||
@@ -27,138 +27,95 @@
|
||||
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
||||
!!$ POSSIBILITY OF SUCH DAMAGE.
|
||||
!!$
|
||||
!!$
|
||||
!!$
|
||||
! File: psb_dspgetrow.f90
|
||||
! Subroutine: psb_dspgetrow
|
||||
! Gets one or more rows from a sparse matrix.
|
||||
! Parameters:
|
||||
|
||||
!*****************************************************************************
|
||||
!* *
|
||||
!* Takes a specified row from matrix A and copies into NZ,IA,JA,VAL in COO *
|
||||
!* format. *
|
||||
!* *
|
||||
!*****************************************************************************
|
||||
subroutine psb_dspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw,bw)
|
||||
subroutine psb_dspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw,append,nzin)
|
||||
! Output is always in COO format
|
||||
use psb_spmat_type
|
||||
use psb_const_mod
|
||||
use psb_getrow_mod
|
||||
use psb_string_mod
|
||||
use psb_serial_mod, psb_protect_name => psb_dspgetrow
|
||||
|
||||
type(psb_dspmat_type), intent(in) :: a
|
||||
integer, intent(in) :: irw
|
||||
integer, intent(out) :: nz
|
||||
integer, intent(inout) :: ia(:), ja(:)
|
||||
real(kind(1.d0)), intent(inout) :: val(:)
|
||||
integer, intent(in), target, optional :: iren(:)
|
||||
integer, intent(in), optional :: lrw
|
||||
integer, intent(out) :: info
|
||||
type(psb_dspmat_type), intent(inout), optional, target :: bw
|
||||
implicit none
|
||||
|
||||
integer :: lrw_, ierr(5), err_act
|
||||
type(psb_dspmat_type), target,save :: b
|
||||
type(psb_dspmat_type), pointer :: b_
|
||||
type(psb_dspmat_type), intent(in) :: a
|
||||
integer, intent(in) :: irw
|
||||
integer, intent(out) :: nz
|
||||
integer, allocatable, intent(inout) :: ia(:), ja(:)
|
||||
real(kind(1.d0)), allocatable, intent(inout) :: val(:)
|
||||
integer,intent(out) :: info
|
||||
logical, intent(in), optional :: append
|
||||
integer, intent(in), optional :: iren(:)
|
||||
integer, intent(in), optional :: lrw, nzin
|
||||
|
||||
integer, pointer :: iren_(:)
|
||||
character(len=20) :: name, ch_err
|
||||
logical :: append_
|
||||
integer :: i,j,k,ip,jp,nr,idx,iret,nzin_, nza, lrw_, irw_, err_act
|
||||
character(len=20) :: name, ch_err
|
||||
|
||||
|
||||
name = 'psb_sp_getrow'
|
||||
name='psb_spgetrow'
|
||||
info = 0
|
||||
!!$ call psb_erractionsave(err_act)
|
||||
!!$ call psb_set_erraction(0)
|
||||
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
irw_ = irw
|
||||
if (present(lrw)) then
|
||||
lrw_ = lrw
|
||||
else
|
||||
lrw_ = irw
|
||||
endif
|
||||
if (lrw_ < irw) then
|
||||
write(0,*) 'SPGETROW input error: fixing lrw',irw,lrw_
|
||||
write(0,*) 'SPGTBLK input error: fixing lrw',irw,lrw_
|
||||
lrw_ = irw
|
||||
end if
|
||||
if (present(bw)) then
|
||||
b_ => bw
|
||||
if (present(append)) then
|
||||
append_=append
|
||||
else
|
||||
b_ => b
|
||||
end if
|
||||
|
||||
call psb_nullify_sp(b_)
|
||||
if (.not.(allocated(b_%aspk).and.allocated(b_%ia1).and.&
|
||||
& allocated(b_%ia2))) then
|
||||
!!$ write(0,*) 'First allocation for B in SPGETROW'
|
||||
call psb_sp_all(lrw_-irw+1,lrw_-irw+1,b_,info)
|
||||
end if
|
||||
if (present(iren)) then
|
||||
call psb_sp_getblk(irw,a,b_,info,iren=iren,lrw=lrw_)
|
||||
else
|
||||
call psb_sp_getblk(irw,a,b_,info,lrw=lrw_)
|
||||
end if
|
||||
|
||||
|
||||
|
||||
if (info /= 0) then
|
||||
info=136
|
||||
ch_err=a%fida(1:3)
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (toupper(b_%fida) /= 'COO') then
|
||||
info=4010
|
||||
ch_err=a%fida(1:3)
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
append_=.false.
|
||||
endif
|
||||
|
||||
nz = b_%infoa(psb_nnz_)
|
||||
|
||||
if (size(ia)>= nz) then
|
||||
ia(1:nz) = b_%ia1(1:nz)
|
||||
if ((append_).and.(present(nzin))) then
|
||||
nzin_ = nzin
|
||||
else
|
||||
info = 135
|
||||
ierr(1) = 4
|
||||
ierr(2) = size(ia)
|
||||
call psb_errpush(info,name,i_err=ierr)
|
||||
goto 9999
|
||||
nzin_ = 0
|
||||
endif
|
||||
|
||||
if (size(ja)>= nz) then
|
||||
ja(1:nz) = b_%ia2(1:nz)
|
||||
if (toupper(a%fida) == 'CSR') then
|
||||
call csr_getrow(irw_,a,nz,ia,ja,val,nzin_,append_,lrw_,info,iren)
|
||||
|
||||
else if (toupper(a%fida) == 'COO') then
|
||||
call coo_getrow(irw_,a,nz,ia,ja,val,nzin_,append_,lrw_,info,iren)
|
||||
|
||||
else if (toupper(a%fida) == 'JAD') then
|
||||
call jad_getrow(irw_,a,nz,ia,ja,val,nzin_,append_,lrw_,info,iren)
|
||||
|
||||
else
|
||||
info = 135
|
||||
ierr(1) = 5
|
||||
ierr(2) = size(ja)
|
||||
call psb_errpush(info,name,i_err=ierr)
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (size(val)>= nz) then
|
||||
val(1:nz) = b_%aspk(1:nz)
|
||||
else
|
||||
info = 135
|
||||
ierr(1) = 6
|
||||
ierr(2) = size(val)
|
||||
call psb_errpush(info,name,i_err=ierr)
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
|
||||
!!$ call psb_sp_free(b,info)
|
||||
|
||||
info=136
|
||||
ch_err=a%fida(1:3)
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (info /= 0) goto 9999
|
||||
!!$ call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
|
||||
9999 continue
|
||||
!!$ call psb_erractionrestore(err_act)
|
||||
call psb_erractionsave(err_act)
|
||||
if (err_act.eq.psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
return
|
||||
|
||||
|
||||
|
||||
end subroutine psb_dspgetrow
|
||||
|
||||
|
||||
+21
-421
@@ -35,8 +35,7 @@
|
||||
!*****************************************************************************
|
||||
!* *
|
||||
!* Takes a specified row from matrix A and copies into matrix B (possibly *
|
||||
!* appending to B). Output is always COO. Input might be anything, once *
|
||||
!* we get to actually write the code..... *
|
||||
!* appending to B). Output is always COO. Input might be anything, *
|
||||
!* *
|
||||
!*****************************************************************************
|
||||
subroutine psb_dspgtblk(irw,a,b,info,append,iren,lrw)
|
||||
@@ -44,6 +43,8 @@ subroutine psb_dspgtblk(irw,a,b,info,append,iren,lrw)
|
||||
! the input format
|
||||
use psb_spmat_type
|
||||
use psb_const_mod
|
||||
use psb_string_mod
|
||||
use psb_serial_mod, psb_protect_name => psb_dspgtblk
|
||||
implicit none
|
||||
|
||||
type(psb_dspmat_type), intent(in) :: a
|
||||
@@ -55,7 +56,6 @@ subroutine psb_dspgtblk(irw,a,b,info,append,iren,lrw)
|
||||
integer, intent(in), optional :: lrw
|
||||
|
||||
logical :: append_
|
||||
integer, pointer :: iren_(:)
|
||||
integer :: i,j,k,ip,jp,nr,idx, nz,iret,nzb, nza, lrw_, irw_, err_act
|
||||
character(len=20) :: name, ch_err
|
||||
|
||||
@@ -74,15 +74,10 @@ subroutine psb_dspgtblk(irw,a,b,info,append,iren,lrw)
|
||||
lrw_ = irw
|
||||
end if
|
||||
if (present(append)) then
|
||||
append_=append
|
||||
append_ = append
|
||||
else
|
||||
append_=.false.
|
||||
append_ = .false.
|
||||
endif
|
||||
if (present(iren)) then
|
||||
iren_ => iren
|
||||
else
|
||||
iren_ => null()
|
||||
end if
|
||||
|
||||
|
||||
if (append_) then
|
||||
@@ -91,25 +86,24 @@ subroutine psb_dspgtblk(irw,a,b,info,append,iren,lrw)
|
||||
nzb = 0
|
||||
b%m = 0
|
||||
b%k = 0
|
||||
b%descra = a%descra
|
||||
endif
|
||||
b%fida = 'COO'
|
||||
|
||||
call psb_sp_getrow(irw,a,nz,b%ia1,b%ia2,b%aspk,info,iren=iren,&
|
||||
& lrw=lrw_,append=append_,nzin=nzb)
|
||||
if (.not.allocated(b%pl)) then
|
||||
allocate(b%pl(1),stat=info)
|
||||
b%pl = 0
|
||||
endif
|
||||
if (.not.allocated(b%pr)) then
|
||||
allocate(b%pr(1),stat=info)
|
||||
b%pr = 0
|
||||
endif
|
||||
b%infoa(psb_nnz_) = nzb+nz
|
||||
b%m = b%m+lrw_-irw+1
|
||||
b%k = max(b%k,a%k)
|
||||
|
||||
if (a%fida == 'CSR') then
|
||||
call csr_dspgtblk(irw_,a,b,append_,iren_,lrw_)
|
||||
|
||||
else if (a%fida == 'COO') then
|
||||
call coo_dspgtblk(irw_,a,b,append_,iren_,lrw_)
|
||||
|
||||
else if (a%fida == 'JAD') then
|
||||
call jad_dspgtblk(irw_,a,b,append_,iren_,lrw_)
|
||||
|
||||
else
|
||||
info=136
|
||||
ch_err=a%fida(1:3)
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
|
||||
!!$ call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
@@ -121,400 +115,6 @@ subroutine psb_dspgtblk(irw,a,b,info,append,iren,lrw)
|
||||
return
|
||||
end if
|
||||
return
|
||||
|
||||
contains
|
||||
|
||||
subroutine csr_dspgtblk(irw,a,b,append,iren,lrw)
|
||||
|
||||
use psb_spmat_type
|
||||
use psb_const_mod
|
||||
implicit none
|
||||
|
||||
type(psb_dspmat_type), intent(in) :: a
|
||||
integer :: irw
|
||||
type(psb_dspmat_type), intent(inout) :: b
|
||||
logical, intent(in) :: append
|
||||
integer, pointer :: iren(:)
|
||||
integer :: lrw
|
||||
|
||||
integer :: idx,i,j ,nr,nz,nzb, row_idx
|
||||
integer, pointer :: indices(:)
|
||||
|
||||
if (append) then
|
||||
nzb = b%infoa(psb_nnz_)
|
||||
else
|
||||
nzb = 0
|
||||
endif
|
||||
|
||||
if (a%pl(1) /= 0) then
|
||||
|
||||
nr = lrw - irw + 1
|
||||
allocate(indices(nr))
|
||||
nz = 0
|
||||
do i=1,nr
|
||||
indices(i)=a%pl(irw+i-1)
|
||||
nz=nz+a%ia2(indices(i)+1)-a%ia2(indices(i))
|
||||
end do
|
||||
|
||||
if (min(size(b%ia1),size(b%ia2),size(b%aspk)) < nzb+nz) then
|
||||
!!$ write(0,*) 'Realloc in gtblk 1',size(b%ia1),size(b%ia2),size(b%aspk),nzb,nz
|
||||
call psb_sp_reall(b,nzb+nz,iret)
|
||||
endif
|
||||
|
||||
k=0
|
||||
if(associated(iren)) then
|
||||
do i=1,nr
|
||||
row_idx=indices(i)
|
||||
do j=a%ia2(row_idx),a%ia2(row_idx+1)-1
|
||||
k = k + 1
|
||||
b%aspk(nzb+k) = a%aspk(j)
|
||||
b%ia1(nzb+k) = iren(row_idx)
|
||||
b%ia2(nzb+k) = iren(a%ia1(j))
|
||||
end do
|
||||
end do
|
||||
else
|
||||
do i=1,nr
|
||||
row_idx=indices(i)
|
||||
do j=a%ia2(row_idx),a%ia2(row_idx+1)-1
|
||||
k = k + 1
|
||||
b%aspk(nzb+k) = a%aspk(j)
|
||||
b%ia1(nzb+k) = row_idx
|
||||
b%ia2(nzb+k) = a%ia1(j)
|
||||
end do
|
||||
end do
|
||||
end if
|
||||
|
||||
b%infoa(psb_nnz_) = nzb+k
|
||||
b%m = b%m+nr
|
||||
b%k = max(b%k,a%k)
|
||||
|
||||
else
|
||||
idx = irw
|
||||
|
||||
if (idx<0) then
|
||||
write(0,*) ' spgtblk Error : idx no good ',idx
|
||||
return
|
||||
end if
|
||||
nr = lrw - irw + 1
|
||||
nz = a%ia2(idx+nr) - a%ia2(idx)
|
||||
|
||||
if (min(size(b%ia1),size(b%ia2),size(b%aspk)) < nzb+nz) then
|
||||
!!$ write(0,*) 'Realloc in gtblk 2',size(b%ia1),size(b%ia2),size(b%aspk),nzb,nz
|
||||
call psb_sp_reall(b,nzb+nz,iret)
|
||||
endif
|
||||
b%fida='COO'
|
||||
|
||||
if (associated(iren)) then
|
||||
k=0
|
||||
do i=irw,lrw
|
||||
do j=a%ia2(i),a%ia2(i+1)-1
|
||||
k = k + 1
|
||||
b%aspk(nzb+k) = a%aspk(j)
|
||||
b%ia1(nzb+k) = iren(i)
|
||||
b%ia2(nzb+k) = iren(a%ia1(j))
|
||||
end do
|
||||
enddo
|
||||
else
|
||||
k=0
|
||||
|
||||
do i=irw,lrw
|
||||
do j=a%ia2(i),a%ia2(i+1)-1
|
||||
k = k + 1
|
||||
b%ia1(nzb+k) = i
|
||||
b%ia2(nzb+k) = a%ia1(j)
|
||||
b%aspk(nzb+k) = a%aspk(j)
|
||||
!!$ write(0,*) 'csr_gtblk: in:',a%aspk(j),i,a%ia1(j)
|
||||
end do
|
||||
enddo
|
||||
end if
|
||||
b%infoa(psb_nnz_) = nzb+nz
|
||||
if (a%pr(1) /= 0) then
|
||||
write(0,*) 'Feeling lazy today, Right Permutation will have to wait'
|
||||
endif
|
||||
b%m = b%m+lrw-irw+1
|
||||
b%k = max(b%k,a%k)
|
||||
|
||||
endif
|
||||
|
||||
end subroutine csr_dspgtblk
|
||||
|
||||
subroutine coo_dspgtblk(irw,a,b,append,iren,lrw)
|
||||
|
||||
use psb_spmat_type
|
||||
use psb_const_mod
|
||||
implicit none
|
||||
|
||||
type(psb_dspmat_type), intent(in) :: a
|
||||
integer :: irw
|
||||
type(psb_dspmat_type), intent(inout) :: b
|
||||
logical, intent(in) :: append
|
||||
integer, pointer :: iren(:)
|
||||
integer :: lrw
|
||||
|
||||
nza = a%infoa(psb_nnz_)
|
||||
if (a%pl(1) /= 0) then
|
||||
write(0,*) 'Fatal error in SPGTBLK: do not feed a permuted mat so far!'
|
||||
idx = -1
|
||||
else
|
||||
idx = irw
|
||||
endif
|
||||
if (idx<0) then
|
||||
write(0,*) ' spgtblk Error : idx no good ',idx
|
||||
return
|
||||
end if
|
||||
|
||||
if (a%infoa(psb_srtd_) == psb_isrtdcoo_) then
|
||||
! In this case we can do a binary search.
|
||||
do
|
||||
call ibsrch(ip,irw,nza,a%ia1)
|
||||
if (ip /= -1) exit
|
||||
irw = irw + 1
|
||||
if (irw > lrw) then
|
||||
write(0,*) 'Warning : did not find any rows. Is this an error? ',irw,lrw,idx
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
if (ip /= -1) then
|
||||
! expand [ip,jp] to contain all row entries.
|
||||
do
|
||||
if (ip < 2) exit
|
||||
if (a%ia1(ip-1) == irw) then
|
||||
ip = ip -1
|
||||
else
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
end if
|
||||
|
||||
do
|
||||
call ibsrch(jp,lrw,nza,a%ia1)
|
||||
if (jp /= -1) exit
|
||||
lrw = lrw - 1
|
||||
if (irw > lrw) then
|
||||
write(0,*) 'Warning : did not find any rows. Is this an error?'
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
if (jp /= -1) then
|
||||
! expand [ip,jp] to contain all row entries.
|
||||
do
|
||||
if (jp == nza) exit
|
||||
if (a%ia1(jp+1) == lrw) then
|
||||
jp = jp + 1
|
||||
else
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
if ((ip /= -1) .and.(jp /= -1)) then
|
||||
! Now do the copy.
|
||||
nz = jp - ip +1
|
||||
if (size(b%ia1) < nzb+nz) then
|
||||
!!$ write(0,*) 'Realloc in gtblk 3',size(b%ia1),size(b%ia2),size(b%aspk),nzb,nz
|
||||
call psb_sp_reall(b,nzb+nz,iret)
|
||||
endif
|
||||
b%fida='COO'
|
||||
if (associated(iren)) then
|
||||
do i=ip,jp
|
||||
nzb = nzb + 1
|
||||
b%aspk(nzb) = a%aspk(i)
|
||||
b%ia1(nzb) = iren(a%ia1(i))
|
||||
b%ia2(nzb) = iren(a%ia2(i))
|
||||
enddo
|
||||
else
|
||||
do i=ip,jp
|
||||
nzb = nzb + 1
|
||||
b%aspk(nzb) = a%aspk(i)
|
||||
b%ia1(nzb) = a%ia1(i)
|
||||
b%ia2(nzb) = a%ia2(i)
|
||||
enddo
|
||||
end if
|
||||
end if
|
||||
|
||||
else
|
||||
|
||||
nz = (nza*(lrw-irw+1))/max(a%m,1)
|
||||
|
||||
if (size(b%ia1) < nzb+nz) then
|
||||
!!$ write(0,*) 'Realloc in gtblk 4',size(b%ia1),size(b%ia2),size(b%aspk),nzb,nz
|
||||
call psb_sp_reall(b,nzb+nz,iret)
|
||||
endif
|
||||
|
||||
if (associated(iren)) then
|
||||
k = 0
|
||||
do i=1,a%infoa(psb_nnz_)
|
||||
if ((a%ia1(i)>=irw).and.(a%ia1(i)<=lrw)) then
|
||||
k = k + 1
|
||||
if (k > nz) then
|
||||
nz = k
|
||||
!!$ write(0,*) 'Realloc in gtblk 5',size(b%ia1),size(b%ia2),size(b%aspk),nzb,nz
|
||||
call psb_sp_reall(b,nzb+nz,iret)
|
||||
end if
|
||||
b%aspk(nzb+k) = a%aspk(i)
|
||||
b%ia1(nzb+k) = iren(a%ia1(i))
|
||||
b%ia2(nzb+k) = iren(a%ia2(i))
|
||||
endif
|
||||
enddo
|
||||
else
|
||||
k = 0
|
||||
do i=1,a%infoa(psb_nnz_)
|
||||
if ((a%ia1(i)>=irw).and.(a%ia1(i)<=lrw)) then
|
||||
k = k + 1
|
||||
if (k > nz) then
|
||||
nz = k
|
||||
!!$ write(0,*) 'Realloc in gtblk 6',size(b%ia1),size(b%ia2),size(b%aspk),nzb,nz
|
||||
call psb_sp_reall(b,nzb+nz,iret)
|
||||
end if
|
||||
b%aspk(nzb+k) = a%aspk(i)
|
||||
b%ia1(nzb+k) = (a%ia1(i))
|
||||
b%ia2(nzb+k) = (a%ia2(i))
|
||||
endif
|
||||
enddo
|
||||
nzb=nzb+k
|
||||
end if
|
||||
end if
|
||||
|
||||
b%infoa(psb_nnz_) = nzb
|
||||
b%m = b%m+lrw-irw+1
|
||||
b%k = max(b%k,a%k)
|
||||
end subroutine coo_dspgtblk
|
||||
|
||||
|
||||
|
||||
|
||||
subroutine jad_dspgtblk(irw,a,b,append,iren,lrw)
|
||||
|
||||
type(psb_dspmat_type), intent(in), target :: a
|
||||
integer :: irw
|
||||
type(psb_dspmat_type), intent(inout) :: b
|
||||
logical, intent(in) :: append
|
||||
integer, pointer :: iren(:)
|
||||
integer :: lrw
|
||||
|
||||
integer, pointer :: ia1(:), ia2(:), ia3(:),&
|
||||
& ja(:), ka(:), indices(:), blks(:)
|
||||
integer :: png, pia, pja, ipx, blk, rb, row, k_pt, npg, col, ng
|
||||
|
||||
|
||||
png = a%ia2(1) ! points to the number of blocks
|
||||
pia = a%ia2(2) ! points to the beginning of ia(3,png)
|
||||
pja = a%ia2(3) ! points to the beginning of ja(:)
|
||||
|
||||
ng = a%ia2(png) ! the number of blocks
|
||||
ja => a%ia2(pja:) ! the array containing the pointers to ka and aspk
|
||||
ka => a%ia1(:) ! the array containing the column indices
|
||||
ia1 => a%ia2(pia:pja-1:3) ! the array containing the first row index of each block
|
||||
ia2 => a%ia2(pia+1:pja-1:3) ! the array containing a pointer to the pos. in ja to the first jad column
|
||||
ia3 => a%ia2(pia+2:pja-1:3) ! the array containing a pointer to the pos. in ja to the first csr column
|
||||
|
||||
if (append) then
|
||||
nzb = b%infoa(psb_nnz_)
|
||||
else
|
||||
nzb = 0
|
||||
endif
|
||||
|
||||
if (a%pl(1) /= 0) then
|
||||
|
||||
nr = lrw - irw + 1
|
||||
allocate(indices(nr),blks(nr))
|
||||
nz = 0
|
||||
|
||||
do i=1,nr
|
||||
indices(i)=a%pl(irw+i-1)
|
||||
j=0
|
||||
blkfnd: do
|
||||
j=j+1
|
||||
if(ia1(j).eq.indices(i)) then
|
||||
blks(i)=j
|
||||
nz=nz+ia3(j)-ia2(j)
|
||||
ipx = ia1(j) ! the first row index of the block
|
||||
rb = indices(i)-ipx ! the row offset within the block
|
||||
row = ia3(j)+rb
|
||||
nz = nz+ja(row+1)-ja(row)
|
||||
exit blkfnd
|
||||
else if(ia1(j).gt.indices(i)) then
|
||||
blks(i)=j-1
|
||||
nz=nz+ia3(j-1)-ia2(j-1)
|
||||
ipx = ia1(j-1) ! the first row index of the block
|
||||
rb = indices(i)-ipx ! the row offset within the block
|
||||
row = ia3(j-1)+rb
|
||||
nz = nz+ja(row+1)-ja(row)
|
||||
exit blkfnd
|
||||
end if
|
||||
end do blkfnd
|
||||
end do
|
||||
|
||||
if (size(b%ia1) < nzb+nz) then
|
||||
call psb_sp_reall(b,nzb+nz,iret)
|
||||
endif
|
||||
|
||||
k=0
|
||||
! cycle over rows
|
||||
do i=1,nr
|
||||
|
||||
! find which block the row belongs to
|
||||
blk = blks(i)
|
||||
|
||||
! extract first part of the row from the jad block
|
||||
ipx = ia1(blk) ! the first row index of the block
|
||||
k_pt= ia2(blk) ! the pointer to the beginning of a column in ja
|
||||
rb = indices(i)-ipx ! the row offset within the block
|
||||
npg = ja(k_pt+1)-ja(k_pt) ! the number of rows in the block
|
||||
|
||||
if(associated(iren))then
|
||||
do col = ia2(blk), ia3(blk)-1
|
||||
k=k+1
|
||||
b%aspk(nzb+k) = a%aspk(ja(col)+rb)
|
||||
b%ia1(nzb+k) = iren(irw+i-1)
|
||||
b%ia2(nzb+k) = iren(ka(ja(col)+rb))
|
||||
end do
|
||||
else
|
||||
do col = ia2(blk), ia3(blk)-1
|
||||
k=k+1
|
||||
b%aspk(nzb+k) = a%aspk(ja(col)+rb)
|
||||
b%ia1(nzb+k) = irw+i-1
|
||||
b%ia2(nzb+k) = ka(ja(col)+rb)
|
||||
end do
|
||||
end if
|
||||
! extract second part of the row from the csr tail
|
||||
row=ia3(blk)+rb
|
||||
if(associated(iren))then
|
||||
do j=ja(row), ja(row+1)-1
|
||||
k=k+1
|
||||
b%aspk(nzb+k) = a%aspk(j)
|
||||
b%ia1(nzb+k) = iren(irw+i-1)
|
||||
b%ia2(nzb+k) = iren(ka(j))
|
||||
end do
|
||||
else
|
||||
do j=ja(row), ja(row+1)-1
|
||||
k=k+1
|
||||
b%aspk(nzb+k) = a%aspk(j)
|
||||
b%ia1(nzb+k) = irw+i-1
|
||||
b%ia2(nzb+k) = ka(j)
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
|
||||
b%infoa(psb_nnz_) = nzb+k
|
||||
b%m = b%m+lrw-irw+1
|
||||
b%k = max(b%k,a%k)
|
||||
b%fida='COO'
|
||||
|
||||
else
|
||||
! There might be some problems
|
||||
info=134
|
||||
ch_err=a%fida(1:3)
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
end if
|
||||
|
||||
|
||||
|
||||
end subroutine jad_dspgtblk
|
||||
|
||||
|
||||
|
||||
|
||||
end subroutine psb_dspgtblk
|
||||
|
||||
@@ -143,7 +143,7 @@ contains
|
||||
main: do i=1,n
|
||||
istart=-1
|
||||
length=0
|
||||
call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info,bw=w)
|
||||
call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info)
|
||||
do jj=1, nazr
|
||||
|
||||
j=iacl(jj)
|
||||
@@ -153,7 +153,7 @@ contains
|
||||
info = 1
|
||||
return
|
||||
endif
|
||||
call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info,bw=w)
|
||||
call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info)
|
||||
do k=1,nbzr
|
||||
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
|
||||
write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn
|
||||
|
||||
@@ -0,0 +1,828 @@
|
||||
module psb_getrow_mod
|
||||
|
||||
interface csr_getrow
|
||||
module procedure csr_dspgtrow, csr_zspgtrow
|
||||
end interface
|
||||
interface coo_getrow
|
||||
module procedure coo_dspgtrow, coo_zspgtrow
|
||||
end interface
|
||||
interface jad_getrow
|
||||
module procedure jad_dspgtrow, jad_zspgtrow
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
subroutine csr_dspgtrow(irw,a,nz,ia,ja,val,nzin,append,lrw,info,iren)
|
||||
|
||||
use psb_spmat_type
|
||||
use psb_const_mod
|
||||
implicit none
|
||||
|
||||
type(psb_dspmat_type), intent(in) :: a
|
||||
integer :: irw
|
||||
integer, intent(out) :: nz
|
||||
integer, allocatable, intent(inout) :: ia(:), ja(:)
|
||||
real(kind(1.d0)), allocatable, intent(inout) :: val(:)
|
||||
integer :: nzin
|
||||
logical, intent(in) :: append
|
||||
integer :: lrw,info
|
||||
integer, optional :: iren(:)
|
||||
|
||||
integer :: idx,i,j, k, nr, row_idx, nzin_
|
||||
integer, allocatable :: indices(:)
|
||||
|
||||
if (append) then
|
||||
nzin_ = nzin
|
||||
else
|
||||
nzin_ = 0
|
||||
endif
|
||||
|
||||
if (a%pl(1) /= 0) then
|
||||
|
||||
nr = lrw - irw + 1
|
||||
allocate(indices(nr))
|
||||
nz = 0
|
||||
do i=1,nr
|
||||
indices(i)=a%pl(irw+i-1)
|
||||
nz=nz+a%ia2(indices(i)+1)-a%ia2(indices(i))
|
||||
end do
|
||||
|
||||
call psb_ensure_size(nzin_+nz,ia,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nz,ja,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nz,val,info)
|
||||
if (info /= 0) return
|
||||
|
||||
k=0
|
||||
if(present(iren)) then
|
||||
do i=1,nr
|
||||
row_idx=indices(i)
|
||||
do j=a%ia2(row_idx),a%ia2(row_idx+1)-1
|
||||
k = k + 1
|
||||
val(nzin_+k) = a%aspk(j)
|
||||
ia(nzin_+k) = iren(row_idx)
|
||||
ja(nzin_+k) = iren(a%ia1(j))
|
||||
end do
|
||||
end do
|
||||
else
|
||||
do i=1,nr
|
||||
row_idx=indices(i)
|
||||
do j=a%ia2(row_idx),a%ia2(row_idx+1)-1
|
||||
k = k + 1
|
||||
val(nzin_+k) = a%aspk(j)
|
||||
ia(nzin_+k) = row_idx
|
||||
ja(nzin_+k) = a%ia1(j)
|
||||
end do
|
||||
end do
|
||||
end if
|
||||
|
||||
else
|
||||
|
||||
idx = irw
|
||||
|
||||
if (idx<0) then
|
||||
write(0,*) ' spgtrow Error : idx no good ',idx
|
||||
return
|
||||
end if
|
||||
nr = lrw - irw + 1
|
||||
nz = a%ia2(idx+nr) - a%ia2(idx)
|
||||
|
||||
call psb_ensure_size(nzin_+nz,ia,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nz,ja,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nz,val,info)
|
||||
if (info /= 0) return
|
||||
|
||||
|
||||
if (present(iren)) then
|
||||
k=0
|
||||
do i=irw,lrw
|
||||
do j=a%ia2(i),a%ia2(i+1)-1
|
||||
k = k + 1
|
||||
val(nzin_+k) = a%aspk(j)
|
||||
ia(nzin_+k) = iren(i)
|
||||
ja(nzin_+k) = iren(a%ia1(j))
|
||||
end do
|
||||
enddo
|
||||
else
|
||||
k=0
|
||||
|
||||
do i=irw,lrw
|
||||
do j=a%ia2(i),a%ia2(i+1)-1
|
||||
k = k + 1
|
||||
ia(nzin_+k) = i
|
||||
ja(nzin_+k) = a%ia1(j)
|
||||
val(nzin_+k) = a%aspk(j)
|
||||
end do
|
||||
enddo
|
||||
end if
|
||||
!!$ if (nz /= k) then
|
||||
!!$ write(0,*) 'csr getrow Size mismatch ',nz,k
|
||||
!!$ endif
|
||||
if (a%pr(1) /= 0) then
|
||||
write(0,*) 'Feeling lazy today, Right Permutation will have to wait'
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
end subroutine csr_dspgtrow
|
||||
|
||||
subroutine coo_dspgtrow(irw,a,nz,ia,ja,val,nzin,append,lrw,info,iren)
|
||||
|
||||
use psb_spmat_type
|
||||
use psb_const_mod
|
||||
implicit none
|
||||
|
||||
type(psb_dspmat_type), intent(in) :: a
|
||||
integer :: irw
|
||||
integer, intent(out) :: nz
|
||||
integer, allocatable, intent(inout) :: ia(:), ja(:)
|
||||
real(kind(1.d0)), allocatable, intent(inout) :: val(:)
|
||||
integer :: nzin
|
||||
logical, intent(in) :: append
|
||||
integer :: lrw,info
|
||||
integer, optional :: iren(:)
|
||||
integer :: nzin_, nza, idx,ip,jp,i,j,k, nzt
|
||||
logical, parameter :: debug=.false.
|
||||
|
||||
nza = a%infoa(psb_nnz_)
|
||||
if (a%pl(1) /= 0) then
|
||||
write(0,*) 'Fatal error in SPGTROW: do not feed a permuted mat so far!'
|
||||
idx = -1
|
||||
else
|
||||
idx = irw
|
||||
endif
|
||||
if (idx<0) then
|
||||
write(0,*) ' spgtrow Error : idx no good ',idx
|
||||
return
|
||||
end if
|
||||
|
||||
if (append) then
|
||||
nzin_ = nzin
|
||||
else
|
||||
nzin_ = 0
|
||||
endif
|
||||
|
||||
if (a%infoa(psb_srtd_) == psb_isrtdcoo_) then
|
||||
! In this case we can do a binary search.
|
||||
if (debug) write(0,*) 'coo_getrow: srtdcoo '
|
||||
do
|
||||
call ibsrch(ip,irw,nza,a%ia1)
|
||||
if (ip /= -1) exit
|
||||
irw = irw + 1
|
||||
if (irw > lrw) then
|
||||
write(0,*) 'Warning : did not find any rows. Is this an error? ',irw,lrw,idx
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
if (ip /= -1) then
|
||||
! expand [ip,jp] to contain all row entries.
|
||||
do
|
||||
if (ip < 2) exit
|
||||
if (a%ia1(ip-1) == irw) then
|
||||
ip = ip -1
|
||||
else
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
end if
|
||||
|
||||
do
|
||||
call ibsrch(jp,lrw,nza,a%ia1)
|
||||
if (jp /= -1) exit
|
||||
lrw = lrw - 1
|
||||
if (irw > lrw) then
|
||||
write(0,*) 'Warning : did not find any rows. Is this an error?'
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
if (jp /= -1) then
|
||||
! expand [ip,jp] to contain all row entries.
|
||||
do
|
||||
if (jp == nza) exit
|
||||
if (a%ia1(jp+1) == lrw) then
|
||||
jp = jp + 1
|
||||
else
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
if (debug) write(0,*) 'coo_getrow: ip jp',ip,jp,nza
|
||||
if ((ip /= -1) .and.(jp /= -1)) then
|
||||
! Now do the copy.
|
||||
nz = jp - ip +1
|
||||
|
||||
call psb_ensure_size(nzin_+nz,ia,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nz,ja,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nz,val,info)
|
||||
if (info /= 0) return
|
||||
|
||||
if (present(iren)) then
|
||||
do i=ip,jp
|
||||
nzin_ = nzin_ + 1
|
||||
val(nzin_) = a%aspk(i)
|
||||
ia(nzin_) = iren(a%ia1(i))
|
||||
ja(nzin_) = iren(a%ia2(i))
|
||||
enddo
|
||||
else
|
||||
do i=ip,jp
|
||||
nzin_ = nzin_ + 1
|
||||
val(nzin_) = a%aspk(i)
|
||||
ia(nzin_) = a%ia1(i)
|
||||
ja(nzin_) = a%ia2(i)
|
||||
enddo
|
||||
end if
|
||||
else
|
||||
nz = 0
|
||||
end if
|
||||
|
||||
else
|
||||
if (debug) write(0,*) 'coo_getrow: unsorted '
|
||||
nzt = (nza*(lrw-irw+1))/max(a%m,1)
|
||||
|
||||
call psb_ensure_size(nzin_+nzt,ia,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nzt,ja,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nzt,val,info)
|
||||
if (info /= 0) return
|
||||
|
||||
if (present(iren)) then
|
||||
k = 0
|
||||
do i=1, a%infoa(psb_nnz_)
|
||||
if ((a%ia1(i)>=irw).and.(a%ia1(i)<=lrw)) then
|
||||
k = k + 1
|
||||
if (k > nzt) then
|
||||
nzt = k
|
||||
call psb_ensure_size(nzin_+nzt,ia,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nzt,ja,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nzt,val,info)
|
||||
if (info /= 0) return
|
||||
end if
|
||||
val(nzin_+k) = a%aspk(i)
|
||||
ia(nzin_+k) = iren(a%ia1(i))
|
||||
ja(nzin_+k) = iren(a%ia2(i))
|
||||
endif
|
||||
enddo
|
||||
else
|
||||
k = 0
|
||||
do i=1,a%infoa(psb_nnz_)
|
||||
if ((a%ia1(i)>=irw).and.(a%ia1(i)<=lrw)) then
|
||||
k = k + 1
|
||||
if (k > nzt) then
|
||||
nzt = k
|
||||
call psb_ensure_size(nzin_+nzt,ia,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nzt,ja,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nzt,val,info)
|
||||
if (info /= 0) return
|
||||
|
||||
end if
|
||||
val(nzin_+k) = a%aspk(i)
|
||||
ia(nzin_+k) = (a%ia1(i))
|
||||
ja(nzin_+k) = (a%ia2(i))
|
||||
endif
|
||||
enddo
|
||||
nzin_=nzin_+k
|
||||
end if
|
||||
nz = k
|
||||
end if
|
||||
|
||||
end subroutine coo_dspgtrow
|
||||
|
||||
|
||||
subroutine jad_dspgtrow(irw,a,nz,ia,ja,val,nzin,append,lrw,info,iren)
|
||||
|
||||
use psb_spmat_type
|
||||
use psb_const_mod
|
||||
implicit none
|
||||
|
||||
type(psb_dspmat_type), intent(in), target :: a
|
||||
integer :: irw
|
||||
integer, intent(out) :: nz
|
||||
integer, allocatable, intent(inout) :: ia(:), ja(:)
|
||||
real(kind(1.d0)), allocatable, intent(inout) :: val(:)
|
||||
integer :: nzin
|
||||
logical, intent(in) :: append
|
||||
integer, optional :: iren(:)
|
||||
integer :: lrw,info
|
||||
|
||||
integer, pointer :: ia1(:), ia2(:), ia3(:),&
|
||||
& ja_(:), ka_(:), indices(:), blks(:)
|
||||
integer :: png, pia, pja, ipx, blk, rb, row, k_pt, npg, col, ng, nzin_,&
|
||||
& i,j,k,nr
|
||||
|
||||
|
||||
png = a%ia2(1) ! points to the number of blocks
|
||||
pia = a%ia2(2) ! points to the beginning of ia(3,png)
|
||||
pja = a%ia2(3) ! points to the beginning of ja(:)
|
||||
|
||||
ng = a%ia2(png) ! the number of blocks
|
||||
ja_ => a%ia2(pja:) ! the array containing the pointers to ka and aspk
|
||||
ka_ => a%ia1(:) ! the array containing the column indices
|
||||
ia1 => a%ia2(pia:pja-1:3) ! the array containing the first row index of each block
|
||||
ia2 => a%ia2(pia+1:pja-1:3) ! the array containing a pointer to the pos. in ja to the first jad column
|
||||
ia3 => a%ia2(pia+2:pja-1:3) ! the array containing a pointer to the pos. in ja to the first csr column
|
||||
|
||||
if (append) then
|
||||
nzin_ = nzin
|
||||
else
|
||||
nzin_ = 0
|
||||
endif
|
||||
|
||||
if (a%pl(1) /= 0) then
|
||||
|
||||
nr = lrw - irw + 1
|
||||
allocate(indices(nr),blks(nr))
|
||||
nz = 0
|
||||
|
||||
do i=1,nr
|
||||
indices(i)=a%pl(irw+i-1)
|
||||
j=0
|
||||
blkfnd: do
|
||||
j=j+1
|
||||
if(ia1(j).eq.indices(i)) then
|
||||
blks(i)=j
|
||||
nz=nz+ia3(j)-ia2(j)
|
||||
ipx = ia1(j) ! the first row index of the block
|
||||
rb = indices(i)-ipx ! the row offset within the block
|
||||
row = ia3(j)+rb
|
||||
nz = nz+ja_(row+1)-ja_(row)
|
||||
exit blkfnd
|
||||
else if(ia1(j).gt.indices(i)) then
|
||||
blks(i)=j-1
|
||||
nz=nz+ia3(j-1)-ia2(j-1)
|
||||
ipx = ia1(j-1) ! the first row index of the block
|
||||
rb = indices(i)-ipx ! the row offset within the block
|
||||
row = ia3(j-1)+rb
|
||||
nz = nz+ja_(row+1)-ja_(row)
|
||||
exit blkfnd
|
||||
end if
|
||||
end do blkfnd
|
||||
end do
|
||||
|
||||
|
||||
call psb_ensure_size(nzin_+nz,ia,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nz,ja,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nz,val,info)
|
||||
if (info /= 0) return
|
||||
|
||||
k=0
|
||||
! cycle over rows
|
||||
do i=1,nr
|
||||
|
||||
! find which block the row belongs to
|
||||
blk = blks(i)
|
||||
|
||||
! extract first part of the row from the jad block
|
||||
ipx = ia1(blk) ! the first row index of the block
|
||||
k_pt= ia2(blk) ! the pointer to the beginning of a column in ja
|
||||
rb = indices(i)-ipx ! the row offset within the block
|
||||
npg = ja_(k_pt+1)-ja_(k_pt) ! the number of rows in the block
|
||||
|
||||
if(present(iren))then
|
||||
do col = ia2(blk), ia3(blk)-1
|
||||
k=k+1
|
||||
val(nzin_+k) = a%aspk(ja_(col)+rb)
|
||||
ia(nzin_+k) = iren(irw+i-1)
|
||||
ja(nzin_+k) = iren(ka_(ja_(col)+rb))
|
||||
end do
|
||||
else
|
||||
do col = ia2(blk), ia3(blk)-1
|
||||
k=k+1
|
||||
val(nzin_+k) = a%aspk(ja_(col)+rb)
|
||||
ia(nzin_+k) = irw+i-1
|
||||
ja(nzin_+k) = ka_(ja_(col)+rb)
|
||||
end do
|
||||
end if
|
||||
! extract second part of the row from the csr tail
|
||||
row=ia3(blk)+rb
|
||||
if(present(iren))then
|
||||
do j=ja_(row), ja_(row+1)-1
|
||||
k=k+1
|
||||
val(nzin_+k) = a%aspk(j)
|
||||
ia(nzin_+k) = iren(irw+i-1)
|
||||
ja(nzin_+k) = iren(ka_(j))
|
||||
end do
|
||||
else
|
||||
do j=ja_(row), ja_(row+1)-1
|
||||
k=k+1
|
||||
val(nzin_+k) = a%aspk(j)
|
||||
ia(nzin_+k) = irw+i-1
|
||||
ja(nzin_+k) = ka_(j)
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
|
||||
else
|
||||
! There might be some problems
|
||||
info=134
|
||||
end if
|
||||
|
||||
end subroutine jad_dspgtrow
|
||||
|
||||
|
||||
subroutine csr_zspgtrow(irw,a,nz,ia,ja,val,nzin,append,lrw,info,iren)
|
||||
|
||||
use psb_spmat_type
|
||||
use psb_const_mod
|
||||
implicit none
|
||||
|
||||
type(psb_zspmat_type), intent(in) :: a
|
||||
integer :: irw
|
||||
integer, intent(out) :: nz
|
||||
integer, allocatable, intent(inout) :: ia(:), ja(:)
|
||||
complex(kind(1.d0)), allocatable, intent(inout) :: val(:)
|
||||
integer :: nzin
|
||||
logical, intent(in) :: append
|
||||
integer :: lrw,info
|
||||
integer, optional :: iren(:)
|
||||
|
||||
integer :: idx,i,j, k, nr, row_idx, nzin_
|
||||
integer, allocatable :: indices(:)
|
||||
|
||||
if (append) then
|
||||
nzin_ = nzin
|
||||
else
|
||||
nzin_ = 0
|
||||
endif
|
||||
|
||||
if (a%pl(1) /= 0) then
|
||||
|
||||
nr = lrw - irw + 1
|
||||
allocate(indices(nr))
|
||||
nz = 0
|
||||
do i=1,nr
|
||||
indices(i)=a%pl(irw+i-1)
|
||||
nz=nz+a%ia2(indices(i)+1)-a%ia2(indices(i))
|
||||
end do
|
||||
|
||||
call psb_ensure_size(nzin_+nz,ia,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nz,ja,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nz,val,info)
|
||||
if (info /= 0) return
|
||||
|
||||
k=0
|
||||
if(present(iren)) then
|
||||
do i=1,nr
|
||||
row_idx=indices(i)
|
||||
do j=a%ia2(row_idx),a%ia2(row_idx+1)-1
|
||||
k = k + 1
|
||||
val(nzin_+k) = a%aspk(j)
|
||||
ia(nzin_+k) = iren(row_idx)
|
||||
ja(nzin_+k) = iren(a%ia1(j))
|
||||
end do
|
||||
end do
|
||||
else
|
||||
do i=1,nr
|
||||
row_idx=indices(i)
|
||||
do j=a%ia2(row_idx),a%ia2(row_idx+1)-1
|
||||
k = k + 1
|
||||
val(nzin_+k) = a%aspk(j)
|
||||
ia(nzin_+k) = row_idx
|
||||
ja(nzin_+k) = a%ia1(j)
|
||||
end do
|
||||
end do
|
||||
end if
|
||||
|
||||
else
|
||||
idx = irw
|
||||
|
||||
if (idx<0) then
|
||||
write(0,*) ' spgtrow Error : idx no good ',idx
|
||||
return
|
||||
end if
|
||||
nr = lrw - irw + 1
|
||||
nz = a%ia2(idx+nr) - a%ia2(idx)
|
||||
|
||||
call psb_ensure_size(nzin_+nz,ia,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nz,ja,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nz,val,info)
|
||||
if (info /= 0) return
|
||||
|
||||
|
||||
if (present(iren)) then
|
||||
k=0
|
||||
do i=irw,lrw
|
||||
do j=a%ia2(i),a%ia2(i+1)-1
|
||||
k = k + 1
|
||||
val(nzin_+k) = a%aspk(j)
|
||||
ia(nzin_+k) = iren(i)
|
||||
ja(nzin_+k) = iren(a%ia1(j))
|
||||
end do
|
||||
enddo
|
||||
else
|
||||
k=0
|
||||
|
||||
do i=irw,lrw
|
||||
do j=a%ia2(i),a%ia2(i+1)-1
|
||||
k = k + 1
|
||||
ia(nzin_+k) = i
|
||||
ja(nzin_+k) = a%ia1(j)
|
||||
val(nzin_+k) = a%aspk(j)
|
||||
end do
|
||||
enddo
|
||||
end if
|
||||
if (a%pr(1) /= 0) then
|
||||
write(0,*) 'Feeling lazy today, Right Permutation will have to wait'
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
end subroutine csr_zspgtrow
|
||||
|
||||
|
||||
subroutine coo_zspgtrow(irw,a,nz,ia,ja,val,nzin,append,lrw,info,iren)
|
||||
|
||||
use psb_spmat_type
|
||||
use psb_const_mod
|
||||
implicit none
|
||||
|
||||
type(psb_zspmat_type), intent(in) :: a
|
||||
integer :: irw
|
||||
integer, intent(out) :: nz
|
||||
integer, allocatable, intent(inout) :: ia(:), ja(:)
|
||||
complex(kind(1.d0)), allocatable, intent(inout) :: val(:)
|
||||
integer :: nzin
|
||||
logical, intent(in) :: append
|
||||
integer :: lrw,info
|
||||
integer, optional :: iren(:)
|
||||
integer :: nzin_, nza, idx,ip,jp,i,j,k, nzt
|
||||
logical, parameter :: debug=.true.
|
||||
|
||||
nza = a%infoa(psb_nnz_)
|
||||
if (a%pl(1) /= 0) then
|
||||
write(0,*) 'Fatal error in SPGTROW: do not feed a permuted mat so far!'
|
||||
idx = -1
|
||||
else
|
||||
idx = irw
|
||||
endif
|
||||
if (idx<0) then
|
||||
write(0,*) ' spgtrow Error : idx no good ',idx
|
||||
return
|
||||
end if
|
||||
|
||||
if (append) then
|
||||
nzin_ = nzin
|
||||
else
|
||||
nzin_ = 0
|
||||
endif
|
||||
|
||||
if (a%infoa(psb_srtd_) == psb_isrtdcoo_) then
|
||||
! In this case we can do a binary search.
|
||||
if (debug) write(0,*) 'coo_getrow: srtdcoo '
|
||||
do
|
||||
call ibsrch(ip,irw,nza,a%ia1)
|
||||
if (ip /= -1) exit
|
||||
irw = irw + 1
|
||||
if (irw > lrw) then
|
||||
write(0,*) 'Warning : did not find any rows. Is this an error? ',irw,lrw,idx
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
if (ip /= -1) then
|
||||
! expand [ip,jp] to contain all row entries.
|
||||
do
|
||||
if (ip < 2) exit
|
||||
if (a%ia1(ip-1) == irw) then
|
||||
ip = ip -1
|
||||
else
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
end if
|
||||
|
||||
do
|
||||
call ibsrch(jp,lrw,nza,a%ia1)
|
||||
if (jp /= -1) exit
|
||||
lrw = lrw - 1
|
||||
if (irw > lrw) then
|
||||
write(0,*) 'Warning : did not find any rows. Is this an error?'
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
if (jp /= -1) then
|
||||
! expand [ip,jp] to contain all row entries.
|
||||
do
|
||||
if (jp == nza) exit
|
||||
if (a%ia1(jp+1) == lrw) then
|
||||
jp = jp + 1
|
||||
else
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
if (debug) write(0,*) 'coo_getrow: ip jp',ip,jp,nza
|
||||
if ((ip /= -1) .and.(jp /= -1)) then
|
||||
! Now do the copy.
|
||||
nz = jp - ip +1
|
||||
|
||||
call psb_ensure_size(nzin_+nz,ia,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nz,ja,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nz,val,info)
|
||||
if (info /= 0) return
|
||||
|
||||
if (present(iren)) then
|
||||
do i=ip,jp
|
||||
nzin_ = nzin_ + 1
|
||||
val(nzin_) = a%aspk(i)
|
||||
ia(nzin_) = iren(a%ia1(i))
|
||||
ja(nzin_) = iren(a%ia2(i))
|
||||
enddo
|
||||
else
|
||||
do i=ip,jp
|
||||
nzin_ = nzin_ + 1
|
||||
val(nzin_) = a%aspk(i)
|
||||
ia(nzin_) = a%ia1(i)
|
||||
ja(nzin_) = a%ia2(i)
|
||||
enddo
|
||||
end if
|
||||
else
|
||||
nz = 0
|
||||
end if
|
||||
|
||||
else
|
||||
if (debug) write(0,*) 'coo_getrow: unsorted '
|
||||
nzt = (nza*(lrw-irw+1))/max(a%m,1)
|
||||
|
||||
call psb_ensure_size(nzin_+nzt,ia,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nzt,ja,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nzt,val,info)
|
||||
if (info /= 0) return
|
||||
|
||||
if (present(iren)) then
|
||||
k = 0
|
||||
do i=1, a%infoa(psb_nnz_)
|
||||
if ((a%ia1(i)>=irw).and.(a%ia1(i)<=lrw)) then
|
||||
k = k + 1
|
||||
if (k > nzt) then
|
||||
nzt = k
|
||||
call psb_ensure_size(nzin_+nzt,ia,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nzt,ja,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nzt,val,info)
|
||||
if (info /= 0) return
|
||||
end if
|
||||
val(nzin_+k) = a%aspk(i)
|
||||
ia(nzin_+k) = iren(a%ia1(i))
|
||||
ja(nzin_+k) = iren(a%ia2(i))
|
||||
endif
|
||||
enddo
|
||||
else
|
||||
k = 0
|
||||
do i=1,a%infoa(psb_nnz_)
|
||||
if ((a%ia1(i)>=irw).and.(a%ia1(i)<=lrw)) then
|
||||
k = k + 1
|
||||
if (k > nzt) then
|
||||
nzt = k
|
||||
call psb_ensure_size(nzin_+nzt,ia,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nzt,ja,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nzt,val,info)
|
||||
if (info /= 0) return
|
||||
|
||||
end if
|
||||
val(nzin_+k) = a%aspk(i)
|
||||
ia(nzin_+k) = (a%ia1(i))
|
||||
ja(nzin_+k) = (a%ia2(i))
|
||||
endif
|
||||
enddo
|
||||
nzin_=nzin_+k
|
||||
end if
|
||||
nz = k
|
||||
end if
|
||||
|
||||
end subroutine coo_zspgtrow
|
||||
|
||||
|
||||
subroutine jad_zspgtrow(irw,a,nz,ia,ja,val,nzin,append,lrw,info,iren)
|
||||
|
||||
use psb_spmat_type
|
||||
use psb_const_mod
|
||||
implicit none
|
||||
|
||||
type(psb_zspmat_type), intent(in), target :: a
|
||||
integer :: irw
|
||||
integer, intent(out) :: nz
|
||||
integer, allocatable, intent(inout) :: ia(:), ja(:)
|
||||
complex(kind(1.d0)), allocatable, intent(inout) :: val(:)
|
||||
integer :: nzin
|
||||
logical, intent(in) :: append
|
||||
integer, optional :: iren(:)
|
||||
integer :: lrw,info
|
||||
|
||||
integer, pointer :: ia1(:), ia2(:), ia3(:),&
|
||||
& ja_(:), ka_(:), indices(:), blks(:)
|
||||
integer :: png, pia, pja, ipx, blk, rb, row, k_pt, npg, col, ng, nzin_,&
|
||||
& i,j,k,nr
|
||||
|
||||
|
||||
png = a%ia2(1) ! points to the number of blocks
|
||||
pia = a%ia2(2) ! points to the beginning of ia(3,png)
|
||||
pja = a%ia2(3) ! points to the beginning of ja(:)
|
||||
|
||||
ng = a%ia2(png) ! the number of blocks
|
||||
ja_ => a%ia2(pja:) ! the array containing the pointers to ka and aspk
|
||||
ka_ => a%ia1(:) ! the array containing the column indices
|
||||
ia1 => a%ia2(pia:pja-1:3) ! the array containing the first row index of each block
|
||||
ia2 => a%ia2(pia+1:pja-1:3) ! the array containing a pointer to the pos. in ja to the first jad column
|
||||
ia3 => a%ia2(pia+2:pja-1:3) ! the array containing a pointer to the pos. in ja to the first csr column
|
||||
|
||||
if (append) then
|
||||
nzin_ = nzin
|
||||
else
|
||||
nzin_ = 0
|
||||
endif
|
||||
|
||||
if (a%pl(1) /= 0) then
|
||||
|
||||
nr = lrw - irw + 1
|
||||
allocate(indices(nr),blks(nr))
|
||||
nz = 0
|
||||
|
||||
do i=1,nr
|
||||
indices(i)=a%pl(irw+i-1)
|
||||
j=0
|
||||
blkfnd: do
|
||||
j=j+1
|
||||
if(ia1(j).eq.indices(i)) then
|
||||
blks(i)=j
|
||||
nz=nz+ia3(j)-ia2(j)
|
||||
ipx = ia1(j) ! the first row index of the block
|
||||
rb = indices(i)-ipx ! the row offset within the block
|
||||
row = ia3(j)+rb
|
||||
nz = nz+ja_(row+1)-ja_(row)
|
||||
exit blkfnd
|
||||
else if(ia1(j).gt.indices(i)) then
|
||||
blks(i)=j-1
|
||||
nz=nz+ia3(j-1)-ia2(j-1)
|
||||
ipx = ia1(j-1) ! the first row index of the block
|
||||
rb = indices(i)-ipx ! the row offset within the block
|
||||
row = ia3(j-1)+rb
|
||||
nz = nz+ja_(row+1)-ja_(row)
|
||||
exit blkfnd
|
||||
end if
|
||||
end do blkfnd
|
||||
end do
|
||||
|
||||
|
||||
call psb_ensure_size(nzin_+nz,ia,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nz,ja,info)
|
||||
if (info==0) call psb_ensure_size(nzin_+nz,val,info)
|
||||
if (info /= 0) return
|
||||
|
||||
k=0
|
||||
! cycle over rows
|
||||
do i=1,nr
|
||||
|
||||
! find which block the row belongs to
|
||||
blk = blks(i)
|
||||
|
||||
! extract first part of the row from the jad block
|
||||
ipx = ia1(blk) ! the first row index of the block
|
||||
k_pt= ia2(blk) ! the pointer to the beginning of a column in ja
|
||||
rb = indices(i)-ipx ! the row offset within the block
|
||||
npg = ja_(k_pt+1)-ja_(k_pt) ! the number of rows in the block
|
||||
|
||||
if(present(iren))then
|
||||
do col = ia2(blk), ia3(blk)-1
|
||||
k=k+1
|
||||
val(nzin_+k) = a%aspk(ja_(col)+rb)
|
||||
ia(nzin_+k) = iren(irw+i-1)
|
||||
ja(nzin_+k) = iren(ka_(ja_(col)+rb))
|
||||
end do
|
||||
else
|
||||
do col = ia2(blk), ia3(blk)-1
|
||||
k=k+1
|
||||
val(nzin_+k) = a%aspk(ja_(col)+rb)
|
||||
ia(nzin_+k) = irw+i-1
|
||||
ja(nzin_+k) = ka_(ja_(col)+rb)
|
||||
end do
|
||||
end if
|
||||
! extract second part of the row from the csr tail
|
||||
row=ia3(blk)+rb
|
||||
if(present(iren))then
|
||||
do j=ja_(row), ja_(row+1)-1
|
||||
k=k+1
|
||||
val(nzin_+k) = a%aspk(j)
|
||||
ia(nzin_+k) = iren(irw+i-1)
|
||||
ja(nzin_+k) = iren(ka_(j))
|
||||
end do
|
||||
else
|
||||
do j=ja_(row), ja_(row+1)-1
|
||||
k=k+1
|
||||
val(nzin_+k) = a%aspk(j)
|
||||
ia(nzin_+k) = irw+i-1
|
||||
ja(nzin_+k) = ka_(j)
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
|
||||
else
|
||||
! There might be some problems
|
||||
info=134
|
||||
end if
|
||||
|
||||
end subroutine jad_zspgtrow
|
||||
|
||||
|
||||
end module psb_getrow_mod
|
||||
@@ -109,7 +109,8 @@ subroutine psb_zspclip(a,b,info,imin,imax,jmin,jmax,rscale,cscale)
|
||||
|
||||
sizeb = psb_sp_get_nnzeros(a)
|
||||
call psb_sp_all(mb,kb,b,sizeb,info)
|
||||
b%fida='COO'
|
||||
b%fida = 'COO'
|
||||
b%descra = a%descra
|
||||
nzb = 0
|
||||
do i=imin_, imax_, irbk
|
||||
nrt = min(irbk,imax_-i+1)
|
||||
|
||||
@@ -27,120 +27,95 @@
|
||||
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
||||
!!$ POSSIBILITY OF SUCH DAMAGE.
|
||||
!!$
|
||||
!!$
|
||||
!!$
|
||||
! File: psb_zspgetrow.f90
|
||||
! Subroutine: psb_zspgetrow
|
||||
! Gets one or more rows from a sparse matrix.
|
||||
! Parameters:
|
||||
|
||||
!*****************************************************************************
|
||||
!* *
|
||||
!* Takes a specified row from matrix A and copies into NZ,IA,JA,VAL in COO *
|
||||
!* format. *
|
||||
!* *
|
||||
!*****************************************************************************
|
||||
subroutine psb_zspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw)
|
||||
subroutine psb_zspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw,append,nzin)
|
||||
! Output is always in COO format
|
||||
use psb_spmat_type
|
||||
use psb_const_mod
|
||||
use psb_getrow_mod
|
||||
use psb_string_mod
|
||||
use psb_serial_mod, psb_protect_name => psb_zspgetrow
|
||||
type(psb_zspmat_type), intent(in) :: a
|
||||
integer, intent(in) :: irw
|
||||
integer, intent(out) :: nz
|
||||
integer, intent(inout) :: ia(:), ja(:)
|
||||
complex(kind(1.d0)), intent(inout) :: val(:)
|
||||
integer, intent(in), target, optional :: iren(:)
|
||||
integer, intent(in), optional :: lrw
|
||||
integer, intent(out) :: info
|
||||
implicit none
|
||||
|
||||
integer :: lrw_, ierr(5), err_act
|
||||
type(psb_zspmat_type) :: b
|
||||
integer, pointer :: iren_(:)
|
||||
character(len=20) :: name, ch_err
|
||||
type(psb_zspmat_type), intent(in) :: a
|
||||
integer, intent(in) :: irw
|
||||
integer, intent(out) :: nz
|
||||
integer, allocatable, intent(inout) :: ia(:), ja(:)
|
||||
complex(kind(1.d0)), allocatable, intent(inout) :: val(:)
|
||||
integer,intent(out) :: info
|
||||
logical, intent(in), optional :: append
|
||||
integer, intent(in), optional :: iren(:)
|
||||
integer, intent(in), optional :: lrw, nzin
|
||||
|
||||
logical :: append_
|
||||
integer :: i,j,k,ip,jp,nr,idx,iret,nzin_, nza, lrw_, irw_, err_act
|
||||
character(len=20) :: name, ch_err
|
||||
|
||||
name='psb_sp_getrow'
|
||||
name='psb_spgetrow'
|
||||
info = 0
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
call psb_set_erraction(0)
|
||||
|
||||
call psb_nullify_sp(b)
|
||||
|
||||
irw_ = irw
|
||||
if (present(lrw)) then
|
||||
lrw_ = lrw
|
||||
else
|
||||
lrw_ = irw
|
||||
endif
|
||||
if (lrw_ < irw) then
|
||||
write(0,*) 'SPGETROW input error: fixing lrw',irw,lrw_
|
||||
write(0,*) 'SPGTBLK input error: fixing lrw',irw,lrw_
|
||||
lrw_ = irw
|
||||
end if
|
||||
call psb_sp_all(lrw_-irw+1,lrw_-irw+1,b,info)
|
||||
if (present(append)) then
|
||||
append_=append
|
||||
else
|
||||
append_=.false.
|
||||
endif
|
||||
|
||||
if (present(iren)) then
|
||||
call psb_sp_getblk(irw,a,b,info,iren=iren,lrw=lrw_)
|
||||
else
|
||||
call psb_sp_getblk(irw,a,b,info,lrw=lrw_)
|
||||
|
||||
if ((append_).and.(present(nzin))) then
|
||||
nzin_ = nzin
|
||||
else
|
||||
nzin_ = 0
|
||||
endif
|
||||
|
||||
if (toupper(a%fida) == 'CSR') then
|
||||
call csr_getrow(irw_,a,nz,ia,ja,val,nzin_,append_,lrw_,info,iren)
|
||||
|
||||
else if (toupper(a%fida) == 'COO') then
|
||||
call coo_getrow(irw_,a,nz,ia,ja,val,nzin_,append_,lrw_,info,iren)
|
||||
|
||||
else if (toupper(a%fida) == 'JAD') then
|
||||
call jad_getrow(irw_,a,nz,ia,ja,val,nzin_,append_,lrw_,info,iren)
|
||||
|
||||
else
|
||||
info=136
|
||||
ch_err=a%fida(1:3)
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
if (info /= 0) then
|
||||
info=136
|
||||
ch_err=a%fida(1:3)
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (toupper(b%fida) /= 'COO') then
|
||||
info=4010
|
||||
ch_err=a%fida(1:3)
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
nz = b%infoa(psb_nnz_)
|
||||
|
||||
if (size(ia)>= nz) then
|
||||
ia(1:nz) = b%ia1(1:nz)
|
||||
else
|
||||
info = 135
|
||||
ierr(1) = 4
|
||||
ierr(2) = size(ia)
|
||||
call psb_errpush(info,name,i_err=ierr)
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (size(ja)>= nz) then
|
||||
ja(1:nz) = b%ia2(1:nz)
|
||||
else
|
||||
info = 135
|
||||
ierr(1) = 5
|
||||
ierr(2) = size(ja)
|
||||
call psb_errpush(info,name,i_err=ierr)
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (size(val)>= nz) then
|
||||
val(1:nz) = b%aspk(1:nz)
|
||||
else
|
||||
info = 135
|
||||
ierr(1) = 6
|
||||
ierr(2) = size(val)
|
||||
call psb_errpush(info,name,i_err=ierr)
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
|
||||
call psb_sp_free(b,info)
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
|
||||
if (info /= 0) goto 9999
|
||||
!!$ call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
|
||||
9999 continue
|
||||
call psb_erractionrestore(err_act)
|
||||
!!$ call psb_erractionrestore(err_act)
|
||||
call psb_erractionsave(err_act)
|
||||
if (err_act.eq.psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
return
|
||||
|
||||
|
||||
|
||||
end subroutine psb_zspgetrow
|
||||
|
||||
|
||||
+28
-423
@@ -32,12 +32,10 @@
|
||||
! Subroutine: psb_zspgtblk
|
||||
! Gets one or more rows from a sparse matrix.
|
||||
! Parameters:
|
||||
|
||||
!*****************************************************************************
|
||||
!* *
|
||||
!* Takes a specified row from matrix A and copies into matrix B (possibly *
|
||||
!* appending to B). Output is always COO. Input might be anything, once *
|
||||
!* we get to actually write the code..... *
|
||||
!* appending to B). Output is always COO. Input might be anything, *
|
||||
!* *
|
||||
!*****************************************************************************
|
||||
subroutine psb_zspgtblk(irw,a,b,info,append,iren,lrw)
|
||||
@@ -45,6 +43,8 @@ subroutine psb_zspgtblk(irw,a,b,info,append,iren,lrw)
|
||||
! the input format
|
||||
use psb_spmat_type
|
||||
use psb_const_mod
|
||||
use psb_string_mod
|
||||
use psb_serial_mod, psb_protect_name => psb_zspgtblk
|
||||
implicit none
|
||||
|
||||
type(psb_zspmat_type), intent(in) :: a
|
||||
@@ -56,13 +56,12 @@ subroutine psb_zspgtblk(irw,a,b,info,append,iren,lrw)
|
||||
integer, intent(in), optional :: lrw
|
||||
|
||||
logical :: append_
|
||||
integer, pointer :: iren_(:)
|
||||
integer :: i,j,k,ip,jp,nr,idx, nz,iret,nzb, nza, lrw_, irw_, err_act
|
||||
character(len=20) :: name, ch_err
|
||||
|
||||
name='psb_spgtblk'
|
||||
info = 0
|
||||
call psb_erractionsave(err_act)
|
||||
!!$ call psb_erractionsave(err_act)
|
||||
|
||||
irw_ = irw
|
||||
if (present(lrw)) then
|
||||
@@ -75,15 +74,10 @@ subroutine psb_zspgtblk(irw,a,b,info,append,iren,lrw)
|
||||
lrw_ = irw
|
||||
end if
|
||||
if (present(append)) then
|
||||
append_=append
|
||||
append_ = append
|
||||
else
|
||||
append_=.false.
|
||||
append_ = .false.
|
||||
endif
|
||||
if (present(iren)) then
|
||||
iren_=>iren
|
||||
else
|
||||
iren_ => null()
|
||||
end if
|
||||
|
||||
|
||||
if (append_) then
|
||||
@@ -92,425 +86,36 @@ subroutine psb_zspgtblk(irw,a,b,info,append,iren,lrw)
|
||||
nzb = 0
|
||||
b%m = 0
|
||||
b%k = 0
|
||||
b%descra = a%descra
|
||||
endif
|
||||
b%fida = 'COO'
|
||||
|
||||
call psb_sp_getrow(irw,a,nz,b%ia1,b%ia2,b%aspk,info,iren=iren,&
|
||||
& lrw=lrw_,append=append_,nzin=nzb)
|
||||
if (.not.allocated(b%pl)) then
|
||||
allocate(b%pl(1),stat=info)
|
||||
b%pl = 0
|
||||
endif
|
||||
if (.not.allocated(b%pr)) then
|
||||
allocate(b%pr(1),stat=info)
|
||||
b%pr = 0
|
||||
endif
|
||||
b%infoa(psb_nnz_) = nzb+nz
|
||||
b%m = b%m+lrw_-irw+1
|
||||
b%k = max(b%k,a%k)
|
||||
|
||||
if (a%fida == 'CSR') then
|
||||
call csr_zspgtblk(irw_,a,b,append_,iren_,lrw_)
|
||||
|
||||
else if (a%fida == 'COO') then
|
||||
call coo_zspgtblk(irw_,a,b,append_,iren_,lrw_)
|
||||
|
||||
else if (a%fida == 'JAD') then
|
||||
call jad_zspgtblk(irw_,a,b,append_,iren_,lrw_)
|
||||
|
||||
else
|
||||
info=136
|
||||
ch_err=a%fida(1:3)
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
!!$ call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
|
||||
9999 continue
|
||||
call psb_erractionrestore(err_act)
|
||||
!!$ call psb_erractionrestore(err_act)
|
||||
call psb_erractionsave(err_act)
|
||||
if (err_act.eq.psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
return
|
||||
|
||||
contains
|
||||
|
||||
subroutine csr_zspgtblk(irw,a,b,append,iren,lrw)
|
||||
|
||||
use psb_spmat_type
|
||||
use psb_const_mod
|
||||
implicit none
|
||||
|
||||
type(psb_zspmat_type), intent(in) :: a
|
||||
integer :: irw
|
||||
type(psb_zspmat_type), intent(inout) :: b
|
||||
logical, intent(in) :: append
|
||||
integer, pointer :: iren(:)
|
||||
integer :: lrw
|
||||
|
||||
integer :: idx,i,j ,nr,nz,nzb, row_idx
|
||||
integer, pointer :: indices(:)
|
||||
|
||||
if (append) then
|
||||
nzb = b%infoa(psb_nnz_)
|
||||
else
|
||||
nzb = 0
|
||||
endif
|
||||
|
||||
if (a%pl(1) /= 0) then
|
||||
|
||||
nr = lrw - irw + 1
|
||||
allocate(indices(nr))
|
||||
nz = 0
|
||||
do i=1,nr
|
||||
indices(i)=a%pl(irw+i-1)
|
||||
nz=nz+a%ia2(indices(i)+1)-a%ia2(indices(i))
|
||||
end do
|
||||
|
||||
if (min(size(b%ia1),size(b%ia2),size(b%aspk)) < nzb+nz) then
|
||||
call psb_sp_reall(b,nzb+nz,iret)
|
||||
endif
|
||||
|
||||
k=0
|
||||
if(associated(iren)) then
|
||||
do i=1,nr
|
||||
row_idx=indices(i)
|
||||
do j=a%ia2(row_idx),a%ia2(row_idx+1)-1
|
||||
k = k + 1
|
||||
b%aspk(nzb+k) = a%aspk(j)
|
||||
b%ia1(nzb+k) = iren(row_idx)
|
||||
b%ia2(nzb+k) = iren(a%ia1(j))
|
||||
end do
|
||||
end do
|
||||
else
|
||||
do i=1,nr
|
||||
row_idx=indices(i)
|
||||
do j=a%ia2(row_idx),a%ia2(row_idx+1)-1
|
||||
k = k + 1
|
||||
b%aspk(nzb+k) = a%aspk(j)
|
||||
b%ia1(nzb+k) = row_idx
|
||||
b%ia2(nzb+k) = a%ia1(j)
|
||||
end do
|
||||
end do
|
||||
end if
|
||||
|
||||
b%infoa(psb_nnz_) = nzb+k
|
||||
b%m = b%m+nr
|
||||
b%k = max(b%k,a%k)
|
||||
|
||||
else
|
||||
idx = irw
|
||||
|
||||
if (idx<0) then
|
||||
write(0,*) ' spgtblk Error : idx no good ',idx
|
||||
return
|
||||
end if
|
||||
nr = lrw - irw + 1
|
||||
nz = a%ia2(idx+nr) - a%ia2(idx)
|
||||
|
||||
if (min(size(b%ia1),size(b%ia2),size(b%aspk)) < nzb+nz) then
|
||||
call psb_sp_reall(b,nzb+nz,iret)
|
||||
endif
|
||||
b%fida='COO'
|
||||
|
||||
if (associated(iren)) then
|
||||
k=0
|
||||
do i=irw,lrw
|
||||
do j=a%ia2(i),a%ia2(i+1)-1
|
||||
k = k + 1
|
||||
b%aspk(nzb+k) = a%aspk(j)
|
||||
b%ia1(nzb+k) = iren(i)
|
||||
b%ia2(nzb+k) = iren(a%ia1(j))
|
||||
end do
|
||||
enddo
|
||||
else
|
||||
k=0
|
||||
|
||||
do i=irw,lrw
|
||||
do j=a%ia2(i),a%ia2(i+1)-1
|
||||
k = k + 1
|
||||
b%ia1(nzb+k) = i
|
||||
b%ia2(nzb+k) = a%ia1(j)
|
||||
b%aspk(nzb+k) = a%aspk(j)
|
||||
!!$ write(0,*) 'csr_gtblk: in:',a%aspk(j),i,a%ia1(j)
|
||||
end do
|
||||
enddo
|
||||
end if
|
||||
b%infoa(psb_nnz_) = nzb+nz
|
||||
if (a%pr(1) /= 0) then
|
||||
write(0,*) 'Feeling lazy today, Right Permutation will have to wait'
|
||||
endif
|
||||
b%m = b%m+lrw-irw+1
|
||||
b%k = max(b%k,a%k)
|
||||
|
||||
endif
|
||||
|
||||
end subroutine csr_zspgtblk
|
||||
|
||||
subroutine coo_zspgtblk(irw,a,b,append,iren,lrw)
|
||||
|
||||
use psb_spmat_type
|
||||
use psb_const_mod
|
||||
implicit none
|
||||
|
||||
type(psb_zspmat_type), intent(in) :: a
|
||||
integer :: irw
|
||||
type(psb_zspmat_type), intent(inout) :: b
|
||||
logical, intent(in) :: append
|
||||
integer, pointer :: iren(:)
|
||||
integer :: lrw
|
||||
|
||||
nza = a%infoa(psb_nnz_)
|
||||
if (a%pl(1) /= 0) then
|
||||
write(0,*) 'Fatal error in SPGTBLK: do not feed a permuted mat so far!'
|
||||
idx = -1
|
||||
else
|
||||
idx = irw
|
||||
endif
|
||||
if (idx<0) then
|
||||
write(0,*) ' spgtblk Error : idx no good ',idx
|
||||
return
|
||||
end if
|
||||
|
||||
if (a%infoa(psb_srtd_) == psb_isrtdcoo_) then
|
||||
! In this case we can do a binary search.
|
||||
do
|
||||
call ibsrch(ip,irw,nza,a%ia1)
|
||||
if (ip /= -1) exit
|
||||
irw = irw + 1
|
||||
if (irw > lrw) then
|
||||
write(0,*) 'Warning : did not find any rows. Is this an error? ',irw,lrw,idx
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
if (ip /= -1) then
|
||||
! expand [ip,jp] to contain all row entries.
|
||||
do
|
||||
if (ip < 2) exit
|
||||
if (a%ia1(ip-1) == irw) then
|
||||
ip = ip -1
|
||||
else
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
end if
|
||||
|
||||
do
|
||||
call ibsrch(jp,lrw,nza,a%ia1)
|
||||
if (jp /= -1) exit
|
||||
lrw = lrw - 1
|
||||
if (irw > lrw) then
|
||||
write(0,*) 'Warning : did not find any rows. Is this an error?'
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
if (jp /= -1) then
|
||||
! expand [ip,jp] to contain all row entries.
|
||||
do
|
||||
if (jp == nza) exit
|
||||
if (a%ia1(jp+1) == lrw) then
|
||||
jp = jp + 1
|
||||
else
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
if ((ip /= -1) .and.(jp /= -1)) then
|
||||
! Now do the copy.
|
||||
nz = jp - ip +1
|
||||
if (size(b%ia1) < nzb+nz) then
|
||||
call psb_sp_reall(b,nzb+nz,iret)
|
||||
endif
|
||||
b%fida='COO'
|
||||
if (associated(iren)) then
|
||||
do i=ip,jp
|
||||
nzb = nzb + 1
|
||||
b%aspk(nzb) = a%aspk(i)
|
||||
b%ia1(nzb) = iren(a%ia1(i))
|
||||
b%ia2(nzb) = iren(a%ia2(i))
|
||||
enddo
|
||||
else
|
||||
do i=ip,jp
|
||||
nzb = nzb + 1
|
||||
b%aspk(nzb) = a%aspk(i)
|
||||
b%ia1(nzb) = a%ia1(i)
|
||||
b%ia2(nzb) = a%ia2(i)
|
||||
enddo
|
||||
end if
|
||||
end if
|
||||
|
||||
else
|
||||
|
||||
nz = (nza*(lrw-irw+1))/max(a%m,1)
|
||||
|
||||
if (size(b%ia1) < nzb+nz) then
|
||||
call psb_sp_reall(b,nzb+nz,iret)
|
||||
endif
|
||||
|
||||
if (associated(iren)) then
|
||||
k = 0
|
||||
do i=1,a%infoa(psb_nnz_)
|
||||
if ((a%ia1(i)>=irw).and.(a%ia1(i)<=lrw)) then
|
||||
k = k + 1
|
||||
if (k > nz) then
|
||||
nz = k
|
||||
call psb_sp_reall(b,nzb+nz,iret)
|
||||
end if
|
||||
b%aspk(nzb+k) = a%aspk(i)
|
||||
b%ia1(nzb+k) = iren(a%ia1(i))
|
||||
b%ia2(nzb+k) = iren(a%ia2(i))
|
||||
endif
|
||||
enddo
|
||||
else
|
||||
k = 0
|
||||
do i=1,a%infoa(psb_nnz_)
|
||||
if ((a%ia1(i)>=irw).and.(a%ia1(i)<=lrw)) then
|
||||
k = k + 1
|
||||
if (k > nz) then
|
||||
nz = k
|
||||
call psb_sp_reall(b,nzb+nz,iret)
|
||||
end if
|
||||
b%aspk(nzb+k) = a%aspk(i)
|
||||
b%ia1(nzb+k) = (a%ia1(i))
|
||||
b%ia2(nzb+k) = (a%ia2(i))
|
||||
endif
|
||||
enddo
|
||||
nzb=nzb+k
|
||||
end if
|
||||
end if
|
||||
|
||||
b%infoa(psb_nnz_) = nzb
|
||||
b%m = b%m+lrw-irw+1
|
||||
b%k = max(b%k,a%k)
|
||||
end subroutine coo_zspgtblk
|
||||
|
||||
|
||||
|
||||
|
||||
subroutine jad_zspgtblk(irw,a,b,append,iren,lrw)
|
||||
|
||||
type(psb_zspmat_type), intent(in),target :: a
|
||||
integer :: irw
|
||||
type(psb_zspmat_type), intent(inout) :: b
|
||||
logical, intent(in) :: append
|
||||
integer, pointer :: iren(:)
|
||||
integer :: lrw
|
||||
|
||||
integer, pointer :: ia1(:), ia2(:), ia3(:),&
|
||||
& ja(:), ka(:), indices(:), blks(:)
|
||||
integer :: png, pia, pja, ipx, blk, rb, row, k_pt, npg, col, ng
|
||||
|
||||
|
||||
png = a%ia2(1) ! points to the number of blocks
|
||||
pia = a%ia2(2) ! points to the beginning of ia(3,png)
|
||||
pja = a%ia2(3) ! points to the beginning of ja(:)
|
||||
|
||||
ng = a%ia2(png) ! the number of blocks
|
||||
ja => a%ia2(pja:) ! the array containing the pointers to ka and aspk
|
||||
ka => a%ia1(:) ! the array containing the column indices
|
||||
ia1 => a%ia2(pia:pja-1:3) ! the array containing the first row index of each block
|
||||
ia2 => a%ia2(pia+1:pja-1:3) ! the array containing a pointer to the pos. in ja to the first jad column
|
||||
ia3 => a%ia2(pia+2:pja-1:3) ! the array containing a pointer to the pos. in ja to the first csr column
|
||||
|
||||
|
||||
if (append) then
|
||||
nzb = b%infoa(psb_nnz_)
|
||||
else
|
||||
nzb = 0
|
||||
endif
|
||||
|
||||
if (a%pl(1) /= 0) then
|
||||
|
||||
nr = lrw - irw + 1
|
||||
allocate(indices(nr),blks(nr))
|
||||
nz = 0
|
||||
|
||||
do i=1,nr
|
||||
indices(i)=a%pl(irw+i-1)
|
||||
j=0
|
||||
blkfnd: do
|
||||
j=j+1
|
||||
if(ia1(j).eq.indices(i)) then
|
||||
blks(i)=j
|
||||
nz=nz+ia3(j)-ia2(j)
|
||||
ipx = ia1(j) ! the first row index of the block
|
||||
rb = indices(i)-ipx ! the row offset within the block
|
||||
row = ia3(j)+rb
|
||||
nz = nz+ja(row+1)-ja(row)
|
||||
exit blkfnd
|
||||
else if(ia1(j).gt.indices(i)) then
|
||||
blks(i)=j-1
|
||||
nz=nz+ia3(j-1)-ia2(j-1)
|
||||
ipx = ia1(j-1) ! the first row index of the block
|
||||
rb = indices(i)-ipx ! the row offset within the block
|
||||
row = ia3(j-1)+rb
|
||||
nz = nz+ja(row+1)-ja(row)
|
||||
exit blkfnd
|
||||
end if
|
||||
end do blkfnd
|
||||
end do
|
||||
|
||||
if (size(b%ia1) < nzb+nz) then
|
||||
call psb_sp_reall(b,nzb+nz,iret)
|
||||
endif
|
||||
|
||||
k=0
|
||||
! cycle over rows
|
||||
do i=1,nr
|
||||
|
||||
! find which block the row belongs to
|
||||
blk = blks(i)
|
||||
|
||||
! extract first part of the row from the jad block
|
||||
ipx = ia1(blk) ! the first row index of the block
|
||||
k_pt= ia2(blk) ! the pointer to the beginning of a column in ja
|
||||
rb = indices(i)-ipx ! the row offset within the block
|
||||
npg = ja(k_pt+1)-ja(k_pt) ! the number of rows in the block
|
||||
|
||||
if(associated(iren))then
|
||||
do col = ia2(blk), ia3(blk)-1
|
||||
k=k+1
|
||||
b%aspk(nzb+k) = a%aspk(ja(col)+rb)
|
||||
b%ia1(nzb+k) = iren(irw+i-1)
|
||||
b%ia2(nzb+k) = iren(ka(ja(col)+rb))
|
||||
end do
|
||||
else
|
||||
do col = ia2(blk), ia3(blk)-1
|
||||
k=k+1
|
||||
b%aspk(nzb+k) = a%aspk(ja(col)+rb)
|
||||
b%ia1(nzb+k) = irw+i-1
|
||||
b%ia2(nzb+k) = ka(ja(col)+rb)
|
||||
end do
|
||||
end if
|
||||
! extract second part of the row from the csr tail
|
||||
row=ia3(blk)+rb
|
||||
if(associated(iren))then
|
||||
do j=ja(row), ja(row+1)-1
|
||||
k=k+1
|
||||
b%aspk(nzb+k) = a%aspk(j)
|
||||
b%ia1(nzb+k) = iren(irw+i-1)
|
||||
b%ia2(nzb+k) = iren(ka(j))
|
||||
end do
|
||||
else
|
||||
do j=ja(row), ja(row+1)-1
|
||||
k=k+1
|
||||
b%aspk(nzb+k) = a%aspk(j)
|
||||
b%ia1(nzb+k) = irw+i-1
|
||||
b%ia2(nzb+k) = ka(j)
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
|
||||
b%infoa(psb_nnz_) = nzb+k
|
||||
b%m = b%m+lrw-irw+1
|
||||
b%k = max(b%k,a%k)
|
||||
b%fida='COO'
|
||||
|
||||
else
|
||||
! There might be some problems
|
||||
info=134
|
||||
ch_err=a%fida(1:3)
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
end if
|
||||
|
||||
|
||||
|
||||
end subroutine jad_zspgtblk
|
||||
|
||||
|
||||
|
||||
|
||||
end subroutine psb_zspgtblk
|
||||
|
||||
|
||||
Reference in New Issue
Block a user