diff --git a/base/modules/psb_serial_mod.f90 b/base/modules/psb_serial_mod.f90 index c392fd7b..635f3449 100644 --- a/base/modules/psb_serial_mod.f90 +++ b/base/modules/psb_serial_mod.f90 @@ -425,28 +425,35 @@ module psb_serial_mod end interface interface psb_sp_getrow - 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 - 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 + + 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 end subroutine psb_dspgetrow - 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 - 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 + + 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 end subroutine psb_zspgetrow end interface diff --git a/base/serial/Makefile b/base/serial/Makefile index f7b02fdd..c4c0e6fd 100644 --- a/base/serial/Makefile +++ b/base/serial/Makefile @@ -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) diff --git a/base/serial/README.serial b/base/serial/README.serial index d0658adb..8ab424e5 100644 --- a/base/serial/README.serial +++ b/base/serial/README.serial @@ -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.... \ No newline at end of file diff --git a/base/serial/psb_dipcoo2csr.f90 b/base/serial/psb_dipcoo2csr.f90 index 2bd533c6..c765dd9c 100644 --- a/base/serial/psb_dipcoo2csr.f90 +++ b/base/serial/psb_dipcoo2csr.f90 @@ -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 diff --git a/base/serial/psb_dnumbmm.f90 b/base/serial/psb_dnumbmm.f90 index 61595928..a71abee0 100644 --- a/base/serial/psb_dnumbmm.f90 +++ b/base/serial/psb_dnumbmm.f90 @@ -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 diff --git a/base/serial/psb_dspclip.f90 b/base/serial/psb_dspclip.f90 index 0ef468dd..65dea802 100644 --- a/base/serial/psb_dspclip.f90 +++ b/base/serial/psb_dspclip.f90 @@ -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) diff --git a/base/serial/psb_dspgetrow.f90 b/base/serial/psb_dspgetrow.f90 index d1d7588f..b5fd42dc 100644 --- a/base/serial/psb_dspgetrow.f90 +++ b/base/serial/psb_dspgetrow.f90 @@ -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 - - integer :: lrw_, ierr(5), err_act - type(psb_dspmat_type), target,save :: b - type(psb_dspmat_type), pointer :: b_ - - integer, pointer :: iren_(:) - character(len=20) :: name, ch_err - - - name = 'psb_sp_getrow' + implicit none + + 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 + + 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_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) - else - info = 135 - ierr(1) = 5 - ierr(2) = size(ja) - call psb_errpush(info,name,i_err=ierr) - goto 9999 - 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) - 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 diff --git a/base/serial/psb_dspgtblk.f90 b/base/serial/psb_dspgtblk.f90 index 0777fd6c..31862114 100644 --- a/base/serial/psb_dspgtblk.f90 +++ b/base/serial/psb_dspgtblk.f90 @@ -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 diff --git a/base/serial/psb_dsymbmm.f90 b/base/serial/psb_dsymbmm.f90 index bcaeab5c..869a777f 100644 --- a/base/serial/psb_dsymbmm.f90 +++ b/base/serial/psb_dsymbmm.f90 @@ -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 diff --git a/base/serial/psb_getrow_mod.f90 b/base/serial/psb_getrow_mod.f90 new file mode 100644 index 00000000..9ee0e19e --- /dev/null +++ b/base/serial/psb_getrow_mod.f90 @@ -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 diff --git a/base/serial/psb_zspclip.f90 b/base/serial/psb_zspclip.f90 index 0f6e9c16..449cc63b 100644 --- a/base/serial/psb_zspclip.f90 +++ b/base/serial/psb_zspclip.f90 @@ -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) diff --git a/base/serial/psb_zspgetrow.f90 b/base/serial/psb_zspgetrow.f90 index d481e102..1a06d54c 100644 --- a/base/serial/psb_zspgetrow.f90 +++ b/base/serial/psb_zspgetrow.f90 @@ -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 - - integer :: lrw_, ierr(5), err_act - type(psb_zspmat_type) :: b - integer, pointer :: iren_(:) - character(len=20) :: name, ch_err - - - name='psb_sp_getrow' + implicit none + + 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_spgetrow' info = 0 - call psb_erractionsave(err_act) - call psb_set_erraction(0) - call psb_nullify_sp(b) + 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 - call psb_sp_all(lrw_-irw+1,lrw_-irw+1,b,info) - - 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 - endif - - nz = b%infoa(psb_nnz_) - - if (size(ia)>= nz) then - ia(1:nz) = b%ia1(1:nz) + if (present(append)) then + append_=append else - info = 135 - ierr(1) = 4 - ierr(2) = size(ia) - call psb_errpush(info,name,i_err=ierr) - goto 9999 + append_=.false. 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) + if ((append_).and.(present(nzin))) then + nzin_ = nzin else - info = 135 - ierr(1) = 6 - ierr(2) = size(val) - call psb_errpush(info,name,i_err=ierr) - goto 9999 + 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) - call psb_sp_free(b,info) - - call psb_erractionrestore(err_act) + 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) 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 diff --git a/base/serial/psb_zspgtblk.f90 b/base/serial/psb_zspgtblk.f90 index 32deedf0..a3b25faa 100644 --- a/base/serial/psb_zspgtblk.f90 +++ b/base/serial/psb_zspgtblk.f90 @@ -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 diff --git a/base/tools/psb_dcdovr.F90 b/base/tools/psb_dcdovr.F90 index 35741f82..5ef862d4 100644 --- a/base/tools/psb_dcdovr.F90 +++ b/base/tools/psb_dcdovr.F90 @@ -55,11 +55,11 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) use psb_realloc_mod use psi_mod #ifdef MPI_MOD - use mpi + use mpi #endif Implicit None #ifdef MPI_H - include 'mpif.h' + include 'mpif.h' #endif ! .. Array Arguments .. @@ -85,7 +85,7 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) & ictxt, lovr, lworks,lworkr, n_row,n_col, int_err(5),& & index_dim,elem_dim, l_tmp_ovr_idx,l_tmp_halo, nztot,nhalo Integer :: counter,counter_h, counter_o, counter_e,idx,gidx,proc,n_elem_recv,& - & n_elem_send,tot_recv,tot_elem,cntov_o,& + & n_elem_send,tot_recv,tot_elem,cntov_o,nz,& & counter_t,n_elem,i_ovr,jj,proc_id,isz, mglob, glx, & & idxr, idxs, lx, iszr, iszs, nxch, nsnd, nrcv,lidx,irsv, extype_ @@ -94,7 +94,6 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) Integer,allocatable :: halo(:),works(:),workr(:),t_halo_in(:),& & t_halo_out(:),temp(:),maskr(:) Integer,allocatable :: brvindx(:),rvsz(:), bsdindx(:),sdsz(:) - Logical,Parameter :: debug=.false. character(len=20) :: name, ch_err @@ -242,10 +241,10 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) gidx = desc_ov%loc_to_glob(idx) - call psb_check_size((cntov_o+3),orig_ovr,info,pad=-1) + call psb_ensure_size((cntov_o+3),orig_ovr,info,pad=-1) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if @@ -342,10 +341,10 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) gidx = desc_ov%loc_to_glob(idx) - call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1) + call psb_ensure_size((counter_o+3),tmp_ovr_idx,info,pad=-1) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if @@ -354,10 +353,10 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) tmp_ovr_idx(counter_o+2)=gidx tmp_ovr_idx(counter_o+3)=-1 counter_o=counter_o+3 - call psb_check_size((counter_h+3),tmp_halo,info,pad=-1) + call psb_ensure_size((counter_h+3),tmp_halo,info,pad=-1) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if @@ -383,10 +382,10 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) & write(0,*) me,i_ovr,'Out of local rows ',& & idx,psb_cd_get_local_rows(Desc_a) - call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1) + call psb_ensure_size((counter_o+3),tmp_ovr_idx,info,pad=-1) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if @@ -401,11 +400,10 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) ! If (i_ovr <= (novr)) Then n_elem = psb_sp_get_nnz_row(idx,a) - - call psb_check_size((idxs+tot_elem+n_elem),works,info) + call psb_ensure_size((idxs+tot_elem+n_elem),works,info) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if @@ -429,9 +427,11 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) goto 9999 end if !!$ write(0,*) me,'Iteration: ',j,i_ovr + Do jj=1,n_elem works(idxs+tot_elem+jj)=desc_ov%loc_to_glob(blk%ia2(jj)) End Do + tot_elem=tot_elem+n_elem End If @@ -518,10 +518,10 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) if (debug) write(0,*) 'ISZR :',iszr if (psb_is_large_desc(desc_ov)) then - call psb_check_size(iszr,maskr,info) + call psb_ensure_size(iszr,maskr,info) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if call psi_idx_cnv(iszr,workr,maskr,desc_ov,info) @@ -553,10 +553,10 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) ! proc_id = temp(i) - call psb_check_size((counter_t+3),t_halo_in,info,pad=-1) + call psb_ensure_size((counter_t+3),t_halo_in,info,pad=-1) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if @@ -585,20 +585,20 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) ! n_col=n_col+1 proc_id=-desc_ov%glob_to_loc(idx)-np-1 - call psb_check_size(n_col,desc_ov%loc_to_glob,info,pad=-1) + call psb_ensure_size(n_col,desc_ov%loc_to_glob,info,pad=-1) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if desc_ov%glob_to_loc(idx)=n_col desc_ov%loc_to_glob(n_col)=idx - call psb_check_size((counter_t+3),t_halo_in,info,pad=-1) + call psb_ensure_size((counter_t+3),t_halo_in,info,pad=-1) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if @@ -665,9 +665,9 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) ! desc_ov%matrix_data(psb_n_row_) = desc_a%matrix_data(psb_n_row_) call psb_transfer(orig_ovr,desc_ov%ovrlap_index,info) - call psb_check_size((counter_h+counter_t+1),tmp_halo,info,pad=-1) + call psb_ensure_size((counter_h+counter_t+1),tmp_halo,info,pad=-1) if (info /= 0) then - call psb_errpush(4010,name,a_err='psb_check_size') + call psb_errpush(4010,name,a_err='psb_ensure_size') goto 9999 end if tmp_halo(counter_h:counter_h+counter_t-1) = t_halo_in(1:counter_t) @@ -679,7 +679,7 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) call psb_errpush(4010,name,a_err='deallocate') goto 9999 end if - + case(psb_ovt_asov_) ! ! Build an overlapped descriptor for Additive Schwarz @@ -691,9 +691,9 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) ! 4. n_row(ov) current. ! 5. n_col(ov) current. ! - call psb_check_size((cntov_o+counter_o+1),orig_ovr,info,pad=-1) + call psb_ensure_size((cntov_o+counter_o+1),orig_ovr,info,pad=-1) if (info /= 0) then - call psb_errpush(4010,name,a_err='psb_check_size') + call psb_errpush(4010,name,a_err='psb_ensure_size') goto 9999 end if orig_ovr(cntov_o:cntov_o+counter_o-1) = tmp_ovr_idx(1:counter_o) diff --git a/base/tools/psb_dsphalo.F90 b/base/tools/psb_dsphalo.F90 index 611972af..cb7dbeab 100644 --- a/base/tools/psb_dsphalo.F90 +++ b/base/tools/psb_dsphalo.F90 @@ -225,6 +225,14 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,cliprow,outfmt,data) tmp%fida='COO' call psb_sp_setifld(psb_spmat_asb_,psb_state_,tmp,info) + + if (debugprt) then + open(20+me) + do i=1, psb_cd_get_local_cols(desc_a) + write(20+me,*) i,desc_a%loc_to_glob(i) + end do + close(20+me) + end if t2 = psb_wtime() l1 = 0 diff --git a/base/tools/psb_zcdovr.F90 b/base/tools/psb_zcdovr.F90 index d692878e..9930529a 100644 --- a/base/tools/psb_zcdovr.F90 +++ b/base/tools/psb_zcdovr.F90 @@ -240,10 +240,10 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) gidx = desc_ov%loc_to_glob(idx) - call psb_check_size((cntov_o+3),orig_ovr,info,pad=-1) + call psb_ensure_size((cntov_o+3),orig_ovr,info,pad=-1) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if @@ -340,10 +340,10 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) gidx = desc_ov%loc_to_glob(idx) - call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1) + call psb_ensure_size((counter_o+3),tmp_ovr_idx,info,pad=-1) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if @@ -352,10 +352,10 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) tmp_ovr_idx(counter_o+2)=gidx tmp_ovr_idx(counter_o+3)=-1 counter_o=counter_o+3 - call psb_check_size((counter_h+3),tmp_halo,info,pad=-1) + call psb_ensure_size((counter_h+3),tmp_halo,info,pad=-1) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if @@ -381,10 +381,10 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) & write(0,*) me,i_ovr,'Out of local rows ',& & idx,psb_cd_get_local_rows(Desc_a) - call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1) + call psb_ensure_size((counter_o+3),tmp_ovr_idx,info,pad=-1) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if @@ -400,10 +400,10 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) If (i_ovr <= (novr)) Then n_elem = psb_sp_get_nnz_row(idx,a) - call psb_check_size((idxs+tot_elem+n_elem),works,info) + call psb_ensure_size((idxs+tot_elem+n_elem),works,info) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if @@ -516,10 +516,10 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) if (debug) write(0,*) 'ISZR :',iszr if (psb_is_large_desc(desc_ov)) then - call psb_check_size(iszr,maskr,info) + call psb_ensure_size(iszr,maskr,info) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if call psi_idx_cnv(iszr,workr,maskr,desc_ov,info) @@ -551,10 +551,10 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) ! proc_id = temp(i) - call psb_check_size((counter_t+3),t_halo_in,info,pad=-1) + call psb_ensure_size((counter_t+3),t_halo_in,info,pad=-1) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if @@ -583,20 +583,20 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) ! n_col=n_col+1 proc_id=-desc_ov%glob_to_loc(idx)-np-1 - call psb_check_size(n_col,desc_ov%loc_to_glob,info,pad=-1) + call psb_ensure_size(n_col,desc_ov%loc_to_glob,info,pad=-1) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if desc_ov%glob_to_loc(idx)=n_col desc_ov%loc_to_glob(n_col)=idx - call psb_check_size((counter_t+3),t_halo_in,info,pad=-1) + call psb_ensure_size((counter_t+3),t_halo_in,info,pad=-1) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_check_size') + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if @@ -663,9 +663,9 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) ! desc_ov%matrix_data(psb_n_row_) = desc_a%matrix_data(psb_n_row_) call psb_transfer(orig_ovr,desc_ov%ovrlap_index,info) - call psb_check_size((counter_h+counter_t+1),tmp_halo,info,pad=-1) + call psb_ensure_size((counter_h+counter_t+1),tmp_halo,info,pad=-1) if (info /= 0) then - call psb_errpush(4010,name,a_err='psb_check_size') + call psb_errpush(4010,name,a_err='psb_ensure_size') goto 9999 end if tmp_halo(counter_h:counter_h+counter_t-1) = t_halo_in(1:counter_t) @@ -689,9 +689,9 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) ! 4. n_row(ov) current. ! 5. n_col(ov) current. ! - call psb_check_size((cntov_o+counter_o+1),orig_ovr,info,pad=-1) + call psb_ensure_size((cntov_o+counter_o+1),orig_ovr,info,pad=-1) if (info /= 0) then - call psb_errpush(4010,name,a_err='psb_check_size') + call psb_errpush(4010,name,a_err='psb_ensure_size') goto 9999 end if orig_ovr(cntov_o:cntov_o+counter_o-1) = tmp_ovr_idx(1:counter_o)