From df830a0b844b5051e6b4e8eda34eb955b1b4d958 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Wed, 14 Jun 2006 15:21:36 +0000 Subject: [PATCH] Added extrow routine. --- src/modules/psb_serial_mod.f90 | 25 +++++ src/modules/psb_spmat_type.f90 | 4 + src/serial/Makefile | 2 +- src/serial/psb_dspextrow.f90 | 162 +++++++++++++++++++++++++++++++++ src/serial/psb_dspgtrow.f90 | 147 +++++++++++++++--------------- src/serial/psb_zspextrow.f90 | 162 +++++++++++++++++++++++++++++++++ test/Fileread/mat_dist.f90 | 153 ++++++++++++++++--------------- 7 files changed, 505 insertions(+), 150 deletions(-) create mode 100644 src/serial/psb_dspextrow.f90 create mode 100644 src/serial/psb_zspextrow.f90 diff --git a/src/modules/psb_serial_mod.f90 b/src/modules/psb_serial_mod.f90 index 421adcb9..7be62e42 100644 --- a/src/modules/psb_serial_mod.f90 +++ b/src/modules/psb_serial_mod.f90 @@ -288,6 +288,31 @@ module psb_serial_mod end subroutine psb_zspgtrow end interface + interface psb_sp_extrow + subroutine psb_dspextrow(irw,a,nz,ia,ja,val,info,iren,lrw) + 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 + end subroutine psb_dspextrow + subroutine psb_zspextrow(irw,a,nz,ia,ja,val,info,iren,lrw) + 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 + end subroutine psb_zspextrow + end interface + interface psb_neigh subroutine psb_dneigh(a,idx,neigh,n,info,lev) use psb_spmat_type diff --git a/src/modules/psb_spmat_type.f90 b/src/modules/psb_spmat_type.f90 index 4b0f7725..b3d61f0a 100644 --- a/src/modules/psb_spmat_type.f90 +++ b/src/modules/psb_spmat_type.f90 @@ -125,6 +125,7 @@ contains type(psb_dspmat_type), intent(inout) :: mat nullify(mat%aspk,mat%ia1,mat%ia2,mat%pl,mat%pr) + mat%infoa(:) = 0 mat%m=0 mat%k=0 mat%fida='' @@ -614,6 +615,7 @@ contains end function psb_dspsizeof + subroutine psb_dsp_free(a,info) implicit none !....Parameters... @@ -649,6 +651,7 @@ contains type(psb_zspmat_type), intent(inout) :: mat nullify(mat%aspk,mat%ia1,mat%ia2,mat%pl,mat%pr) + mat%infoa(:) = 0 mat%m=0 mat%k=0 mat%fida='' @@ -1167,5 +1170,6 @@ contains End Subroutine psb_zsp_free + end module psb_spmat_type diff --git a/src/serial/Makefile b/src/serial/Makefile index d8151009..888fd428 100644 --- a/src/serial/Makefile +++ b/src/serial/Makefile @@ -6,7 +6,7 @@ FOBJS = psb_cest.o psb_dcoins.o psb_dcsdp.o psb_dcsmm.o psb_dcsmv.o \ psb_dfixcoo.o psb_dipcoo2csr.o psb_dipcsr2coo.o psb_dneigh.o \ psb_dnumbmm.o psb_drwextd.o psb_dspgtdiag.o psb_dspgtrow.o \ psb_dspinfo.o psb_dspscal.o psb_dsymbmm.o psb_dtransp.o \ - psb_dipcoo2csc.o lsame.o\ + psb_dipcoo2csc.o psb_dspextrow.o lsame.o psb_zspextrow.o\ psb_zcsmm.o psb_zcsmv.o psb_zspgtdiag.o psb_zspgtrow.o\ psb_zcsnmi.o psb_zcsrws.o psb_zcssm.o psb_zcssv.o psb_zcsdp.o\ psb_zfixcoo.o psb_zipcoo2csr.o psb_zipcsr2coo.o psb_zipcoo2csc.o \ diff --git a/src/serial/psb_dspextrow.f90 b/src/serial/psb_dspextrow.f90 new file mode 100644 index 00000000..95fbbdfa --- /dev/null +++ b/src/serial/psb_dspextrow.f90 @@ -0,0 +1,162 @@ +!!$ +!!$ Parallel Sparse BLAS v2.0 +!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the PSBLAS group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! File: psb_dspextrow.f90 +! Subroutine: psb_dspextrow +! 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_dspextrow(irw,a,nz,ia,ja,val,info,iren,lrw) + use psb_spmat_type + use psb_string_mod + 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 + interface psb_spgtrow + subroutine psb_dspgtrow(irw,a,b,info,append,iren,lrw) + ! Output is always in COO format into B, irrespective of + ! the input format + use psb_spmat_type + use psb_const_mod + implicit none + + type(psb_dspmat_type), intent(in) :: a + integer, intent(in) :: irw + type(psb_dspmat_type), intent(inout) :: b + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), target, optional :: iren(:) + integer, intent(in), optional :: lrw + end subroutine psb_dspgtrow + end interface + + integer :: lrw_, ierr(2), err_act + type(psb_dspmat_type) :: b + integer, pointer :: iren_(:) + character(len=20) :: name, ch_err + + + name='psb_sp_extrow' + info = 0 + call psb_erractionsave(err_act) + call psb_set_erraction(0) + + call psb_nullify_sp(b) + + if (present(lrw)) then + lrw_ = lrw + else + lrw_ = irw + endif + if (lrw_ < irw) then + write(0,*) 'SPEXTROW 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_spgtrow(irw,a,b,info,iren=iren,lrw=lrw_) + else + call psb_spgtrow(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) + 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) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.act_abort) then + call psb_error() + return + end if + return + + +end subroutine psb_dspextrow + diff --git a/src/serial/psb_dspgtrow.f90 b/src/serial/psb_dspgtrow.f90 index aafd1238..74a8114d 100644 --- a/src/serial/psb_dspgtrow.f90 +++ b/src/serial/psb_dspgtrow.f90 @@ -141,96 +141,95 @@ contains integer, pointer :: indices(:) if (append) then - nzb = b%infoa(psb_nnz_) + nzb = b%infoa(psb_nnz_) else - nzb = 0 + nzb = 0 endif if (a%pl(1) /= 0) then - nr = lrw - irw + 1 - allocate(indices(nr)) - do i=1,nr - indices(i)=a%pl(irw+i-1) - nz=nz+a%ia2(indices(i)+1)-a%ia2(indices(i)) - end do + nr = lrw - irw + 1 + allocate(indices(nr)) + 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 + 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 + 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 - 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 + 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 if + end do + end if - b%infoa(psb_nnz_) = nzb+k - b%m = b%m+nr - b%k = max(b%k,a%k) + 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,*) ' spgtrow 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' + idx = irw - 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 + 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) - do i=irw,lrw + 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' - do j=a%ia2(i),a%ia2(i+1)-1 - k = k + 1 - b%aspk(nzb+k) = a%aspk(j) - b%ia1(nzb+k) = i - b%ia2(nzb+k) = a%ia1(j) + 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_gtrow: 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) + 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 diff --git a/src/serial/psb_zspextrow.f90 b/src/serial/psb_zspextrow.f90 new file mode 100644 index 00000000..7a5adbd1 --- /dev/null +++ b/src/serial/psb_zspextrow.f90 @@ -0,0 +1,162 @@ +!!$ +!!$ Parallel Sparse BLAS v2.0 +!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the PSBLAS group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! File: psb_zspextrow.f90 +! Subroutine: psb_zspextrow +! 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_zspextrow(irw,a,nz,ia,ja,val,info,iren,lrw) + use psb_spmat_type + use psb_string_mod + 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 + interface psb_spgtrow + subroutine psb_zspgtrow(irw,a,b,info,append,iren,lrw) + ! Output is always in COO format into B, irrespective of + ! the input format + use psb_spmat_type + use psb_const_mod + implicit none + + type(psb_zspmat_type), intent(in) :: a + integer, intent(in) :: irw + type(psb_zspmat_type), intent(inout) :: b + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), target, optional :: iren(:) + integer, intent(in), optional :: lrw + end subroutine psb_zspgtrow + end interface + + integer :: lrw_, ierr(2), err_act + type(psb_zspmat_type) :: b + integer, pointer :: iren_(:) + character(len=20) :: name, ch_err + + + name='psb_sp_extrow' + info = 0 + call psb_erractionsave(err_act) + call psb_set_erraction(0) + + call psb_nullify_sp(b) + + if (present(lrw)) then + lrw_ = lrw + else + lrw_ = irw + endif + if (lrw_ < irw) then + write(0,*) 'SPEXTROW 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_spgtrow(irw,a,b,info,iren=iren,lrw=lrw_) + else + call psb_spgtrow(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) + 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) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.act_abort) then + call psb_error() + return + end if + return + + +end subroutine psb_zspextrow + diff --git a/test/Fileread/mat_dist.f90 b/test/Fileread/mat_dist.f90 index 56c1c964..f66e6053 100644 --- a/test/Fileread/mat_dist.f90 +++ b/test/Fileread/mat_dist.f90 @@ -122,7 +122,7 @@ contains integer :: np, iam integer :: ircode, length_row, i_count, j_count,& & k_count, blockdim, root, liwork, nrow, ncol, nnzero, nrhs,& - & i,j,k, ll, isize, iproc, nnr, err, err_act, int_err(5) + & i,j,k, ll, nz, isize, iproc, nnr, err, err_act, int_err(5) integer, pointer :: iwork(:) character :: afmt*5, atyp*5 integer, allocatable :: irow(:),icol(:) @@ -248,15 +248,19 @@ contains if (iam == root) then - ll=0 - do j = i_count, j_count - do k=a_glob%ia2(j),a_glob%ia2(j+1)-1 - ll = ll+1 - irow(ll) = j - icol(ll) = a_glob%ia1(k) - val(ll) = a_glob%aspk(k) - end do - enddo + ll = 0 + do i= i_count, j_count-1 + call psb_sp_extrow(i,a_glob,nz,& + & irow(ll+1:),icol(ll+1:),val(ll+1:), info) + if (info /= 0) then + if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then + write(0,*) 'Allocation failure? This should not happen!' + end if + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + ll = ll + nz + end do if (iproc == iam) then call psb_spins(ll,irow,icol,val,a,desc_a,info) @@ -332,15 +336,19 @@ contains k_count = iwork(j_count) if (iam == root) then - ll=0 - do j = i_count, i_count - do k=a_glob%ia2(j),a_glob%ia2(j+1)-1 - ll = ll+1 - irow(ll) = j - icol(ll) = a_glob%ia1(k) - val(ll) = a_glob%aspk(k) - end do - enddo + ll = 0 + do i= i_count, i_count + call psb_sp_extrow(i,a_glob,nz,& + & irow(ll+1:),icol(ll+1:),val(ll+1:), info) + if (info /= 0) then + if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then + write(0,*) 'Allocation failure? This should not happen!' + end if + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + ll = ll + nz + end do if (k_count == iam) then @@ -550,7 +558,7 @@ contains integer :: np, iam integer :: ircode, length_row, i_count, j_count,& & k_count, blockdim, root, liwork, nrow, ncol, nnzero, nrhs,& - & i,j,k, ll, isize, iproc, nnr, err, err_act, int_err(5) + & i,j,k, ll, nz, isize, iproc, nnr, err, err_act, int_err(5) integer, pointer :: iwork(:) character :: afmt*5, atyp*5 integer, allocatable :: irow(:),icol(:) @@ -661,26 +669,20 @@ contains nnr = j_count - i_count if (iam == root) then - ll = a_glob%ia2(j_count)-a_glob%ia2(i_count) - if (ll > size(val)) then - deallocate(val,irow,icol) - allocate(val(ll),irow(ll),icol(ll),stat=info) - if(info/=0) then - info=4010 - ch_err='Allocate' + + ll = 0 + do i= i_count, j_count-1 + call psb_sp_extrow(i,a_glob,nz,& + & irow(ll+1:),icol(ll+1:),val(ll+1:), info) + if (info /= 0) then + if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then + write(0,*) 'Allocation failure? This should not happen!' + end if call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - - endif - k = a_glob%ia2(i_count) - do i= i_count, j_count-1 - do j = a_glob%ia2(i),a_glob%ia2(i+1)-1 - irow(j-k+1) = i - icol(j-k+1) = a_glob%ia1(j) - val(j-k+1) = a_glob%aspk(j) - end do - enddo + ll = ll + nz + end do if (iproc == iam) then call psb_spins(ll,irow,icol,val,a,desc_a,info) @@ -897,7 +899,7 @@ contains integer :: np, iam integer :: ircode, length_row, i_count, j_count,& & k_count, blockdim, root, liwork, nrow, ncol, nnzero, nrhs,& - & i,j,k, ll, isize, iproc, nnr, err, err_act, int_err(5) + & i,j,k, ll, nz, isize, iproc, nnr, err, err_act, int_err(5) integer, pointer :: iwork(:) character :: afmt*5, atyp*5 integer, allocatable :: irow(:),icol(:) @@ -1021,15 +1023,19 @@ contains if (iam == root) then - ll=0 - do j = i_count, j_count - do k=a_glob%ia2(j),a_glob%ia2(j+1)-1 - ll = ll+1 - irow(ll) = j - icol(ll) = a_glob%ia1(k) - val(ll) = a_glob%aspk(k) - end do - enddo + ll = 0 + do i= i_count, j_count-1 + call psb_sp_extrow(i,a_glob,nz,& + & irow(ll+1:),icol(ll+1:),val(ll+1:), info) + if (info /= 0) then + if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then + write(0,*) 'Allocation failure? This should not happen!' + end if + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + ll = ll + nz + end do if (iproc == iam) then call psb_spins(ll,irow,icol,val,a,desc_a,info) @@ -1105,15 +1111,19 @@ contains k_count = iwork(j_count) if (iam == root) then - ll=0 - do j = i_count, i_count - do k=a_glob%ia2(j),a_glob%ia2(j+1)-1 - ll = ll+1 - irow(ll) = j - icol(ll) = a_glob%ia1(k) - val(ll) = a_glob%aspk(k) - end do - enddo + ll = 0 + do i= i_count, i_count + call psb_sp_extrow(i,a_glob,nz,& + & irow(ll+1:),icol(ll+1:),val(ll+1:), info) + if (info /= 0) then + if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then + write(0,*) 'Allocation failure? This should not happen!' + end if + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + ll = ll + nz + end do if (k_count == iam) then @@ -1322,7 +1332,7 @@ contains integer :: np, iam integer :: ircode, length_row, i_count, j_count,& & k_count, blockdim, root, liwork, nrow, ncol, nnzero, nrhs,& - & i,j,k, ll, isize, iproc, nnr, err, err_act, int_err(5) + & i,j,k, ll, nz, isize, iproc, nnr, err, err_act, int_err(5) integer, pointer :: iwork(:) character :: afmt*5, atyp*5 integer, allocatable :: irow(:),icol(:) @@ -1363,7 +1373,7 @@ contains call psb_errpush(info,name) goto 9999 endif - + nnzero = size(a_glob%aspk) nrhs = 1 end if @@ -1433,26 +1443,19 @@ contains nnr = j_count - i_count if (iam == root) then - ll = a_glob%ia2(j_count)-a_glob%ia2(i_count) - if (ll > size(val)) then - deallocate(val,irow,icol) - allocate(val(ll),irow(ll),icol(ll),stat=info) - if(info/=0) then - info=4010 - ch_err='Allocate' + ll = 0 + do i= i_count, j_count-1 + call psb_sp_extrow(i,a_glob,nz,& + & irow(ll+1:),icol(ll+1:),val(ll+1:), info) + if (info /= 0) then + if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then + write(0,*) 'Allocation failure? This should not happen!' + end if call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - - endif - k = a_glob%ia2(i_count) - do i= i_count, j_count-1 - do j = a_glob%ia2(i),a_glob%ia2(i+1)-1 - irow(j-k+1) = i - icol(j-k+1) = a_glob%ia1(j) - val(j-k+1) = a_glob%aspk(j) - end do - enddo + ll = ll + nz + end do if (iproc == iam) then call psb_spins(ll,irow,icol,val,a,desc_a,info)