From d1c2be389b20345cd7d6d8fb8039d66289d8dd3e Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Mon, 9 Apr 2007 17:19:14 +0000 Subject: [PATCH] Fixed neigh to use getrow, not inside code, to avoid pointless code replications. Defined imsru: sort unique. --- base/serial/aux/Makefile | 2 +- base/serial/aux/imsru.f90 | 61 ++++++++ base/serial/psb_dneigh.f90 | 273 +++++------------------------------- base/serial/psb_zneigh.f90 | 275 +++++-------------------------------- 4 files changed, 131 insertions(+), 480 deletions(-) create mode 100644 base/serial/aux/imsru.f90 diff --git a/base/serial/aux/Makefile b/base/serial/aux/Makefile index 4400b293..10938cf7 100644 --- a/base/serial/aux/Makefile +++ b/base/serial/aux/Makefile @@ -4,7 +4,7 @@ include ../../../Make.inc # FOBJS = isr.o isrx.o msort_up.o msort_dw.o\ - isaperm.o ibsrch.o issrch.o imsr.o imsrx.o + isaperm.o ibsrch.o issrch.o imsr.o imsrx.o imsru.o OBJS=$(FOBJS) diff --git a/base/serial/aux/imsru.f90 b/base/serial/aux/imsru.f90 new file mode 100644 index 00000000..dbd273fd --- /dev/null +++ b/base/serial/aux/imsru.f90 @@ -0,0 +1,61 @@ +!!$ +!!$ 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: imsru.f90 + ! Subroutine: + ! Parameters: +subroutine imsru(n,x,idir,nout) + use psb_serial_mod + integer :: n, idir,nout + integer :: x(n) + + + integer, allocatable :: iaux(:) + + integer :: iswap, iret, info, lp, k + integer :: lswap + + if (n<0) then +!!$ write(0,*) 'Error: IMSR: N<0' + return + endif + + if (n<=1) return + + call imsr(n,x,idir) + nout = 1 + do k=2,n + if (x(k) /= x(nout)) then + nout = nout + 1 + x(nout) = x(k) + endif + enddo + return +end subroutine imsru diff --git a/base/serial/psb_dneigh.f90 b/base/serial/psb_dneigh.f90 index ebd23c53..6a1afb59 100644 --- a/base/serial/psb_dneigh.f90 +++ b/base/serial/psb_dneigh.f90 @@ -37,6 +37,7 @@ subroutine psb_dneigh(a,idx,neigh,n,info,lev) use psb_realloc_mod use psb_const_mod use psb_spmat_type + use psb_serial_mod, psb_protect_name => psb_dneigh implicit none @@ -47,9 +48,11 @@ subroutine psb_dneigh(a,idx,neigh,n,info,lev) integer, optional :: lev ! level of neighbours to find integer, allocatable :: tmpn(:) - integer :: level, dim, i, j, k, r, c, brow,& - & elem_pt, ii, n1, col_idx, ne, err_act, nn, nidx + integer :: lev_, dim, i, j, k, r, c, brow,nl, ifl,ill,& + & elem_pt, ii, n1, col_idx, ne, err_act, nn, nidx,ntl character(len=20) :: name, ch_err + integer, allocatable :: ia(:), ja(:) + real(kind(1.d0)), allocatable :: val(:) name='psb_dneigh' info = 0 @@ -58,34 +61,39 @@ subroutine psb_dneigh(a,idx,neigh,n,info,lev) n = 0 info = 0 if(present(lev)) then - if(lev.le.2) then - level=lev - else - write(0,'("Too many levels!!!")') - return - endif + lev_=lev else - level=1 + lev_=1 end if - call psb_dneigh1l(a,idx,neigh,n) - - if(level.eq.2) then - n1=n - allocate(tmpn(max(10,2*n))) - if (size(neigh) 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 - - - i=0 - dim=0 - blkfnd: do - i=i+1 - if(ia1(i).eq.iidx) then - blk=i - dim=dim+ia3(i)-ia2(i) - ipx = ia1(i) ! the first row index of the block - rb = iidx-ipx ! the row offset within the block - row = ia3(i)+rb - dim = dim+ja(row+1)-ja(row) - exit blkfnd - else if(ia1(i).gt.iidx) then - blk=i-1 - dim=dim+ia3(i-1)-ia2(i-1) - ipx = ia1(i-1) ! the first row index of the block - rb = iidx-ipx ! the row offset within the block - row = ia3(i-1)+rb - dim = dim+ja(row+1)-ja(row) - exit blkfnd - end if - end do blkfnd - - if(dim.gt.size(neigh)) call psb_realloc(dim,neigh,info) - - 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 = iidx-ipx ! the row offset within the block - npg = ja(k_pt+1)-ja(k_pt) ! the number of rows in the block - - k=0 - do col = ia2(blk), ia3(blk)-1 - k=k+1 - neigh(k) = ka(ja(col)+rb) - end do - - ! extract second part of the row from the csr tail - row=ia3(blk)+rb - do j=ja(row), ja(row+1)-1 - k=k+1 - neigh(k) = ka(j) - end do - - n=k - - end subroutine jad_dneigh1l - - subroutine psb_dneigh1l(a,idx,neigh,n) - - use psb_realloc_mod - use psb_const_mod - use psb_spmat_type - use psb_string_mod - implicit none - - - type(psb_dspmat_type), intent(in) :: a ! the sparse matrix - integer, intent(in) :: idx ! the index whose neighbours we want to find - integer, intent(out) :: n ! the number of neighbours and the info - integer, allocatable :: neigh(:) ! the neighbours - - - select case(toupper(a%fida(1:3))) - case('CSR') - call csr_dneigh1l(a,idx,neigh,n) - case('COO') - call coo_dneigh1l(a,idx,neigh,n) - case('JAD') - call jad_dneigh1l(a,idx,neigh,n) - end select - - end subroutine psb_dneigh1l end subroutine psb_dneigh diff --git a/base/serial/psb_zneigh.f90 b/base/serial/psb_zneigh.f90 index 7a22723c..ef3ec432 100644 --- a/base/serial/psb_zneigh.f90 +++ b/base/serial/psb_zneigh.f90 @@ -37,19 +37,22 @@ subroutine psb_zneigh(a,idx,neigh,n,info,lev) use psb_realloc_mod use psb_const_mod use psb_spmat_type + use psb_serial_mod, psb_protect_name => psb_zneigh implicit none type(psb_zspmat_type), intent(in) :: a ! the sparse matrix integer, intent(in) :: idx ! the index whose neighbours we want to find integer, intent(out) :: n, info ! the number of neighbours and the info - integer, allocatable :: neigh(:) ! the neighbours + integer, allocatable :: neigh(:) ! the neighbours integer, optional :: lev ! level of neighbours to find integer, allocatable :: tmpn(:) - integer :: level, dim, i, j, k, r, c, brow,& - & elem_pt, ii, n1, col_idx, ne, err_act, nn, nidx + integer :: lev_, dim, i, j, k, r, c, brow,nl, ifl,ill,& + & elem_pt, ii, n1, col_idx, ne, err_act, nn, nidx,ntl character(len=20) :: name, ch_err + integer, allocatable :: ia(:), ja(:) + complex(kind(1.d0)), allocatable :: val(:) name='psb_zneigh' info = 0 @@ -58,34 +61,39 @@ subroutine psb_zneigh(a,idx,neigh,n,info,lev) n = 0 info = 0 if(present(lev)) then - if(lev.le.2) then - level=lev - else - write(0,'("Too many levels!!!")') - return - endif + lev_ = lev else - level=1 + lev_=1 end if - call psb_zneigh1l(a,idx,neigh,n) - - if(level.eq.2) then - n1=n - allocate(tmpn(max(10,2*n))) - if (size(neigh) 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 - - - i=0 - dim=0 - blkfnd: do - i=i+1 - if(ia1(i).eq.iidx) then - blk=i - dim=dim+ia3(i)-ia2(i) - ipx = ia1(i) ! the first row index of the block - rb = iidx-ipx ! the row offset within the block - row = ia3(i)+rb - dim = dim+ja(row+1)-ja(row) - exit blkfnd - else if(ia1(i).gt.iidx) then - blk=i-1 - dim=dim+ia3(i-1)-ia2(i-1) - ipx = ia1(i-1) ! the first row index of the block - rb = iidx-ipx ! the row offset within the block - row = ia3(i-1)+rb - dim = dim+ja(row+1)-ja(row) - exit blkfnd - end if - end do blkfnd - - if(dim.gt.size(neigh)) call psb_realloc(dim,neigh,info) - - 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 = iidx-ipx ! the row offset within the block - npg = ja(k_pt+1)-ja(k_pt) ! the number of rows in the block - - k=0 - do col = ia2(blk), ia3(blk)-1 - k=k+1 - neigh(k) = ka(ja(col)+rb) - end do - - ! extract second part of the row from the csr tail - row=ia3(blk)+rb - do j=ja(row), ja(row+1)-1 - k=k+1 - neigh(k) = ka(j) - end do - - n=k - - end subroutine jad_zneigh1l - - subroutine psb_zneigh1l(a,idx,neigh,n) - - use psb_realloc_mod - use psb_const_mod - use psb_spmat_type - use psb_string_mod - implicit none - - - type(psb_zspmat_type), intent(in) :: a ! the sparse matrix - integer, intent(in) :: idx ! the index whose neighbours we want to find - integer, intent(out) :: n ! the number of neighbours and the info - integer, allocatable :: neigh(:) ! the neighbours - - - select case(toupper(a%fida(1:3))) - case('CSR') - call csr_zneigh1l(a,idx,neigh,n) - case('COO') - call coo_zneigh1l(a,idx,neigh,n) - case('JAD') - call jad_zneigh1l(a,idx,neigh,n) - end select - - end subroutine psb_zneigh1l end subroutine psb_zneigh