fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Fixed makefile and use statements for psb_msort, psb_qsort etc.

Browse Source
psblas3-type-indexed
Salvatore Filippone 18 years ago
parent 49fce2d579
commit faa6867bb8
24 changed files with 1290 additions and 475 deletions
Show Stats Download Patch File Download Diff File

3
base/internals/psi_crea_bnd_elem.f90
Unescape Escape View File

@ -32,6 +32,7 @@ subroutine psi_crea_bnd_elem(bndel,desc_a,info)
use psb_realloc_mod use psb_realloc_mod
use psb_descriptor_type use psb_descriptor_type
use psb_error_mod use psb_error_mod
use psb_serial_mod
implicit none implicit none
integer, allocatable :: bndel(:) integer, allocatable :: bndel(:)
@ -67,7 +68,7 @@ subroutine psi_crea_bnd_elem(bndel,desc_a,info)
enddo enddo
if (i>0) then if (i>0) then
call isr(i,work) call psb_msort(work(1:i))
j=1 j=1
irv = work(1) irv = work(1)
do k=2, i do k=2, i

3
base/internals/psi_crea_ovr_elem.f90
Unescape Escape View File

@ -33,6 +33,7 @@ subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem,info)
use psb_realloc_mod use psb_realloc_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
use psb_serial_mod
implicit none implicit none
! ...parameter arrays.... ! ...parameter arrays....
@ -154,7 +155,7 @@ subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem,info)
i=i+2*desc_overlap(i)+2 i=i+2*desc_overlap(i)+2
enddo enddo
if (nel > 0) then if (nel > 0) then
call imsr(nel,telem(:,1)) call psb_msort(telem(1:nel,1))
iel = telem(1,1) iel = telem(1,1)
telem(1,2) = 2 telem(1,2) = 2
ix = 1 ix = 1

3
base/internals/psi_ldsc_pre_halo.f90
Unescape Escape View File

@ -113,7 +113,8 @@ subroutine psi_ldsc_pre_halo(desc,ext_hv,info)
idx = desc%hashv(i) idx = desc%hashv(i)
nh = desc%hashv(i+1) - desc%hashv(i) nh = desc%hashv(i+1) - desc%hashv(i)
if (nh > 1) then if (nh > 1) then
call imsrx(nh,desc%glb_lc(idx,1),desc%glb_lc(idx,2),1) call psb_msort(desc%glb_lc(idx:idx+nh-1,1),&
& ix=desc%glb_lc(idx:idx+nh-1,2),flag=psb_sort_keep_idx_)
end if end if
end do end do

5
base/internals/srtlist.f
Unescape Escape View File

@ -78,6 +78,7 @@ C
*********************************************************************** ***********************************************************************
SUBROUTINE SRTLIST(DEP_LIST,DL_LDA,LDL,NP,dg,dgp,upd, SUBROUTINE SRTLIST(DEP_LIST,DL_LDA,LDL,NP,dg,dgp,upd,
+ edges,idx,ich,INFO) + edges,idx,ich,INFO)
use psb_serial_mod
IMPLICIT NONE IMPLICIT NONE
INTEGER NP, DL_LDA, INFO INTEGER NP, DL_LDA, INFO
INTEGER DEP_LIST(DL_LDA,*), LDL(*),DG(*), DGP(*), IDX(*), INTEGER DEP_LIST(DL_LDA,*), LDL(*),DG(*), DGP(*), IDX(*),
@ -131,7 +132,7 @@ c$$$ write(0,*) 'SRTLIST Input :',i,ip
DGP(I) = -(DG(EDGES(1,I))+DG(EDGES(2,I))) DGP(I) = -(DG(EDGES(1,I))+DG(EDGES(2,I)))
ENDDO ENDDO
CALL ISRX(NEDGES-IST+1,DGP(IST),IDX(IST)) call psb_msort(dgp(ist:nedges),ix=idx(ist:nedges))
I1 = IST I1 = IST
NCH = 0 NCH = 0
DO I = IST, NEDGES DO I = IST, NEDGES
@ -156,7 +157,7 @@ c$$$ write(0,*) 'SRTLIST Input :',i,ip
info = 30 info = 30
return return
ENDIF ENDIF
CALL ISR(NCH,ICH) call psb_msort(ich(1:nch))
DO I=1, NCH DO I=1, NCH
ISWAP(1) = EDGES(1,IST) ISWAP(1) = EDGES(1,IST)
ISWAP(2) = EDGES(2,IST) ISWAP(2) = EDGES(2,IST)

4
base/modules/Makefile
Unescape Escape View File

@ -6,9 +6,9 @@ MODULES = psb_realloc_mod.o psb_string_mod.o psb_spmat_type.o \
psb_error_mod.o \ psb_error_mod.o \
psb_const_mod.o \ psb_const_mod.o \
psb_comm_mod.o psb_psblas_mod.o psi_mod.o \ psb_comm_mod.o psb_psblas_mod.o psi_mod.o \
psb_check_mod.o blacs_env.o psb_gps_mod.o psb_check_mod.o psb_gps_mod.o
# psb_methd_mod.o psb_prec_type.o psb_prec_mod.o # psb_methd_mod.o psb_prec_type.o psb_prec_mod.o blacs_env.o
LIBMOD=psb_base_mod$(.mod) LIBMOD=psb_base_mod$(.mod)
MPFOBJS=psb_penv_mod.o MPFOBJS=psb_penv_mod.o

61
base/modules/blacs_env.F90
Unescape Escape View File

@ -1,61 +0,0 @@
subroutine psb_set_coher(ictxt,isvch)
integer :: ictxt, isvch
! Ensure global coherence for convergence checks.
#ifdef NETLIB_BLACS
Call blacs_get(ictxt,16,isvch)
Call blacs_set(ictxt,16,1)
#endif
#ifdef ESSL_BLACS
! Do nothing: ESSL does coherence by default,
! and does not handle req=16
#endif
end subroutine psb_set_coher
subroutine psb_restore_coher(ictxt,isvch)
integer :: ictxt, isvch
! Ensure global coherence for convergence checks.
#ifdef NETLIB_BLACS
Call blacs_set(ictxt,16,isvch)
#endif
#ifdef ESSL_BLACS
! Do nothing: ESSL does coherence by default,
! and does not handle req=16
#endif
end subroutine psb_restore_coher
subroutine psb_get_mpicomm(ictxt,comm)
integer :: ictxt, comm
#if !defined(SERIAL_MPI)
call blacs_get(ictxt,10,comm)
#else
comm = ictxt
#endif
end subroutine psb_get_mpicomm
subroutine psb_get_rank(rank,ictxt,id)
integer :: rank,ictxt, id
integer :: blacs_pnum
#if defined(SERIAL_MPI)
rank = 0
#else
rank = blacs_pnum(ictxt,id,0)
#endif
end subroutine psb_get_rank
#if defined(ESSL_BLACS) || defined(SERIAL_MPI)
!
! Need these, as they are not in the ESSL implementation
! of the BLACS.
!
integer function krecvid(contxt,proc_to_comm,myrow)
integer contxt,proc_to_comm,myrow
krecvid=32766
return
end function krecvid
integer function ksendid(contxt,proc_to_comm,myrow)
integer contxt,proc_to_comm,myrow
ksendid=32766
return
end function ksendid
#endif

88
base/modules/psb_penv_mod.F90
Unescape Escape View File

@ -71,7 +71,7 @@ module psb_penv_mod
module procedure psb_ibcasts, psb_ibcastv, psb_ibcastm,& module procedure psb_ibcasts, psb_ibcastv, psb_ibcastm,&
& psb_dbcasts, psb_dbcastv, psb_dbcastm,& & psb_dbcasts, psb_dbcastv, psb_dbcastm,&
& psb_zbcasts, psb_zbcastv, psb_zbcastm,& & psb_zbcasts, psb_zbcastv, psb_zbcastm,&
& psb_hbcasts, psb_lbcasts & psb_hbcasts, psb_lbcasts, psb_lbcastv
end interface end interface
@ -544,6 +544,29 @@ contains
end subroutine psb_lbcasts end subroutine psb_lbcasts
subroutine psb_lbcastv(ictxt,dat,root)
use mpi
integer, intent(in) :: ictxt
logical, intent(inout) :: dat(:)
integer, intent(in), optional :: root
integer :: iam, np, root_,icomm,info
#if !defined(SERIAL_MPI)
if (present(root)) then
root_ = root
else
root_ = 0
endif
call psb_info(ictxt,iam,np)
call psb_get_mpicomm(ictxt,icomm)
call mpi_bcast(dat,size(dat),MPI_LOGICAL,root_,icomm,info)
#endif
end subroutine psb_lbcastv
subroutine psb_imaxs(ictxt,dat,root) subroutine psb_imaxs(ictxt,dat,root)
use mpi use mpi
@ -4787,4 +4810,67 @@ contains
#endif #endif
subroutine psb_set_coher(ictxt,isvch)
integer :: ictxt, isvch
! Ensure global coherence for convergence checks.
#ifdef NETLIB_BLACS
Call blacs_get(ictxt,16,isvch)
Call blacs_set(ictxt,16,1)
#endif
#ifdef ESSL_BLACS
! Do nothing: ESSL does coherence by default,
! and does not handle req=16
#endif
end subroutine psb_set_coher
subroutine psb_restore_coher(ictxt,isvch)
integer :: ictxt, isvch
! Ensure global coherence for convergence checks.
#ifdef NETLIB_BLACS
Call blacs_set(ictxt,16,isvch)
#endif
#ifdef ESSL_BLACS
! Do nothing: ESSL does coherence by default,
! and does not handle req=16
#endif
end subroutine psb_restore_coher
subroutine psb_get_mpicomm(ictxt,comm)
integer :: ictxt, comm
#if !defined(SERIAL_MPI)
call blacs_get(ictxt,10,comm)
#else
comm = ictxt
#endif
end subroutine psb_get_mpicomm
subroutine psb_get_rank(rank,ictxt,id)
integer :: rank,ictxt, id
integer :: blacs_pnum
#if defined(SERIAL_MPI)
rank = 0
#else
rank = blacs_pnum(ictxt,id,0)
#endif
end subroutine psb_get_rank
#if defined(ESSL_BLACS) || defined(SERIAL_MPI)
!
! Need these, as they are not in the ESSL implementation
! of the BLACS.
!
integer function krecvid(contxt,proc_to_comm,myrow)
integer contxt,proc_to_comm,myrow
krecvid=32766
return
end function krecvid
integer function ksendid(contxt,proc_to_comm,myrow)
integer contxt,proc_to_comm,myrow
ksendid=32766
return
end function ksendid
#endif
end module psb_penv_mod end module psb_penv_mod

128
base/modules/psb_serial_mod.f90
Unescape Escape View File

@ -419,6 +419,134 @@ module psb_serial_mod
end subroutine psb_zspgetrow end subroutine psb_zspgetrow
end interface end interface
interface psb_msort
module procedure imsort
end interface
interface psb_qsort
module procedure iqsort
end interface
integer, parameter :: psb_sort_up_=1, psb_sort_down_=-1
integer, parameter :: psb_sort_ovw_idx_=0, psb_sort_keep_idx_=1
contains
subroutine imsort(x,ix,dir,flag)
use psb_error_mod
implicit none
integer, intent(inout) :: x(:)
integer, optional, intent(in) :: dir, flag
integer, optional, intent(inout) :: ix(:)
integer :: dir_, flag_, n, err_act
character(len=20) :: name
name='psb_msort'
if (present(dir)) then
dir_ = dir
else
dir_= psb_sort_up_
end if
select case(dir_)
case( psb_sort_up_, psb_sort_down_)
! OK keep going
case default
call psb_errpush(30,name,i_err=(/3,dir_,0,0,0/))
goto 9999
end select
n = size(x)
if (present(ix)) then
if (size(ix) < n) then
call psb_errpush(35,name,i_err=(/2,size(ix),0,0,0/))
goto 9999
end if
if (present(flag)) then
flag_ = flag
else
flag_ = psb_sort_ovw_idx_
end if
select case(flag_)
case( psb_sort_ovw_idx_, psb_sort_keep_idx_)
! OK keep going
case default
call psb_errpush(30,name,i_err=(/4,flag_,0,0,0/))
goto 9999
end select
call imsrx(n,x,ix,dir_,flag_)
else
call imsr(n,x,dir_)
end if
9999 continue
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
end subroutine imsort
subroutine iqsort(x,ix,dir,flag)
use psb_error_mod
implicit none
integer, intent(inout) :: x(:)
integer, optional, intent(in) :: dir, flag
integer, optional, intent(inout) :: ix(:)
integer :: dir_, flag_, n, err_act
character(len=20) :: name
name='psb_qsort'
if (present(dir)) then
dir_ = dir
else
dir_= psb_sort_up_
end if
select case(dir_)
case( psb_sort_up_, psb_sort_down_)
! OK keep going
case default
call psb_errpush(30,name,i_err=(/3,dir_,0,0,0/))
goto 9999
end select
n = size(x)
if (present(ix)) then
if (size(ix) < n) then
call psb_errpush(35,name,i_err=(/2,size(ix),0,0,0/))
goto 9999
end if
if (present(flag)) then
flag_ = flag
else
flag_ = psb_sort_ovw_idx_
end if
select case(flag_)
case( psb_sort_ovw_idx_, psb_sort_keep_idx_)
! OK keep going
case default
call psb_errpush(30,name,i_err=(/4,flag_,0,0,0/))
goto 9999
end select
call isrx(n,x,ix,dir_,flag_)
else
call isr(n,x,dir_)
end if
9999 continue
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
end subroutine iqsort

2
base/serial/aux/Makefile
Unescape Escape View File

@ -4,7 +4,7 @@ include ../../../Make.inc
# #
FOBJS = isr.o isrx.o \ FOBJS = isr.o isrx.o \
mrgsrt.o isaperm.o ibsrch.o issrch.o imsr.o imsrx.o mrgsrt.o mrgsrtd.o isaperm.o ibsrch.o issrch.o imsr.o imsrx.o
OBJS=$(FOBJS) OBJS=$(FOBJS)

13
base/serial/aux/imsr.f90
Unescape Escape View File

@ -31,10 +31,12 @@
! File: imsr.f90 ! File: imsr.f90
! Subroutine: ! Subroutine:
! Parameters: ! Parameters:
subroutine imsr(n,x) subroutine imsr(n,x,idir)
integer :: n use psb_serial_mod
integer :: n, idir
integer :: x(n) integer :: x(n)
integer, allocatable :: iaux(:) integer, allocatable :: iaux(:)
integer :: iswap, iret, info, lp, k integer :: iswap, iret, info, lp, k
@ -53,8 +55,11 @@ subroutine imsr(n,x)
return return
endif endif
if (idir==psb_sort_up_) then
call mrgsrt(n,x,iaux,iret) call mrgsrt(n,x,iaux,iret)
else
call mrgsrtd(n,x,iaux,iret)
end if
if (iret == 0) then if (iret == 0) then
lp = iaux(0) lp = iaux(0)

28
base/serial/aux/imsrx.f90
Unescape Escape View File

@ -31,8 +31,9 @@
! File: imsrx.f90 ! File: imsrx.f90
! Subroutine: ! Subroutine:
! Parameters: ! Parameters:
subroutine imsrx(n,x,indx,flag) subroutine imsrx(n,x,indx,idir,flag)
integer :: n, flag use psb_serial_mod
integer :: n,idir,flag
integer :: x(n) integer :: x(n)
integer :: indx(n) integer :: indx(n)
@ -47,12 +48,14 @@ subroutine imsrx(n,x,indx,flag)
endif endif
if (n==0) return if (n==0) return
if (n==1) then
if (flag == 0) then if (flag == psb_sort_ovw_idx_) then
indx(1) = 1 do k=1,n
endif indx(k) = k
return enddo
endif end if
if (n==1) return
allocate(iaux(0:n+1),stat=info) allocate(iaux(0:n+1),stat=info)
if (info/=0) then if (info/=0) then
@ -60,12 +63,11 @@ subroutine imsrx(n,x,indx,flag)
return return
endif endif
if (flag == 0) then if (idir == psb_sort_up_) then
do k=1,n call mrgsrt(n,x,iaux,iret)
indx(k) = k else
enddo call mrgsrtd(n,x,iaux,iret)
end if end if
call mrgsrt(n,x,iaux,iret)
if (iret /= 1) then if (iret /= 1) then
lp = iaux(0) lp = iaux(0)

177
base/serial/aux/isr.f
Unescape Escape View File

@ -1,177 +0,0 @@
C
C Parallel Sparse BLAS v2.0
C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
C Alfredo Buttari University of Rome Tor Vergata
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C 1. Redistributions of source code must retain the above copyright
C notice, this list of conditions and the following disclaimer.
C 2. Redistributions in binary form must reproduce the above copyright
C notice, this list of conditions, and the following disclaimer in the
C documentation and/or other materials provided with the distribution.
C 3. The name of the PSBLAS group or the names of its contributors may
C not be used to endorse or promote products derived from this
C software without specific written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
C POSSIBILITY OF SUCH DAMAGE.
C
C
SUBROUTINE ISR(N,X)
C
C Quicksort.
C Adapted from a number of sources, including Don Knuth's TAOCP.
C
C .. Scalar Arguments ..
INTEGER N
C ..
C .. Array Arguments ..
INTEGER X(N)
C ..
C .. Local Scalars ..
INTEGER I, J, XX, ILX, IUX, ISTP, PIV, LPIV
INTEGER IT1, N1, N2
INTEGER MAXSTACK,NPARMS,ITHRS
PARAMETER (MAXSTACK=64,NPARMS=3,ITHRS=16)
INTEGER ISTACK(NPARMS,MAXSTACK)
C ..
C
C Small inputs will only get through insertion sort.
C
IF (N.GT.ITHRS) THEN
C
C Init stack pointer
C
ISTP = 1
ISTACK(1,ISTP) = 1
ISTACK(2,ISTP) = N
DO WHILE (ISTP.GT.0)
ILX = ISTACK(1,ISTP)
IUX = ISTACK(2,ISTP)
ISTP = ISTP - 1
c$$$ write(0,*) 'Debug 1: ',ilx,iux
C
C Choose a pivot with median-of-three heuristics, leave it
C in the LPIV location
C
I = ILX
J = IUX
LPIV = (I+J)/2
PIV = X(LPIV)
IF (PIV.LT.X(I)) THEN
IT1 = X(I)
X(I) = X(LPIV)
X(LPIV) = IT1
PIV = X(LPIV)
ENDIF
IF (PIV.GT.X(J)) THEN
IT1 = X(J)
X(J) = X(LPIV)
X(LPIV) = IT1
PIV = X(LPIV)
ENDIF
IF (PIV.LT.X(I)) THEN
IT1 = X(I)
X(I) = X(LPIV)
X(LPIV) = IT1
PIV = X(LPIV)
ENDIF
C
C Now PIV is correct; place it into first location
IT1 = X(I)
X(I) = X(LPIV)
X(LPIV) = IT1
I = ILX - 1
J = IUX + 1
130 CONTINUE
I = I + 1
XK = X(I)
IF (XK.LT.PIV) GOTO 130
C
C Ensure finite termination for next loop
C
IT1 = XK
X(I) = PIV
140 CONTINUE
J = J - 1
XK = X(J)
IF (XK.GT.PIV) GOTO 140
X(I) = IT1
150 CONTINUE
IF (J.GT.I) THEN
IT1 = X(I)
X(I) = X(J)
X(J) = IT1
GO TO 130
END IF
if (i.eq.ilx) then
if (x(i).ne.piv) then
write(0,*) 'Should never ever get here????!!!!'
stop
endif
i = i + 1
endif
N1 = (I-1)-ILX+1
N2 = IUX-(I)+1
IF (N1.GT.N2) THEN
if (n1.gt.ithrs) then
ISTP = ISTP + 1
ISTACK(1,ISTP) = ILX
ISTACK(2,ISTP) = I-1
endif
if (n2.gt.ithrs) then
ISTP = ISTP + 1
ISTACK(1,ISTP) = I
ISTACK(2,ISTP) = IUX
endif
ELSE
if (n2.gt.ithrs) then
ISTP = ISTP + 1
ISTACK(1,ISTP) = I
ISTACK(2,ISTP) = IUX
endif
if (n1.gt.ithrs) then
ISTP = ISTP + 1
ISTACK(1,ISTP) = ILX
ISTACK(2,ISTP) = I-1
endif
ENDIF
ENDDO
ENDIF
DO J=N-1,1,-1
IF (X(J+1).LT.X(J)) THEN
XX = X(J)
I=J+1
100 CONTINUE
X(I-1) = X(I)
I = I+1
IF ((I.LE.N)) then
if (X(I).LT.XX) GOTO 100
endif
X(I-1) = XX
ENDIF
ENDDO
RETURN
END

355
base/serial/aux/isr.f90
Unescape Escape View File

@ -0,0 +1,355 @@
!!$
!!$ 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.
!!$
!!$
subroutine isr(n,x,dir)
use psb_serial_mod
!
! Quicksort.
! Adapted from a number of sources, including Don Knuth's TAOCP.
!
! .. Scalar Arguments ..
integer, intent(in) :: n, dir
integer :: x(n)
! ..
! .. Local Scalars ..
integer i, j, xx, ilx, iux, istp, piv, lpiv
integer it1, n1, n2
integer, parameter :: maxstack=64,nparms=3,ithrs=16
integer :: istack(nparms,maxstack)
! ..
!
! small inputs will only get through insertion sort.
!
select case(dir)
case(psb_sort_up_)
if (n > ithrs) then
!
! Init stack pointer
!
istp = 1
istack(1,istp) = 1
istack(2,istp) = n
do
if (istp <= 0) exit
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
!
i = ilx
j = iux
lpiv = (i+j)/2
piv = x(lpiv)
if (piv < x(i)) then
it1 = x(i)
x(i) = x(lpiv)
x(lpiv) = it1
piv = x(lpiv)
endif
if (piv > x(j)) then
it1 = x(j)
x(j) = x(lpiv)
x(lpiv) = it1
piv = x(lpiv)
endif
if (piv < x(i)) then
it1 = x(i)
x(i) = x(lpiv)
x(lpiv) = it1
piv = x(lpiv)
endif
!
! now piv is correct; place it into first location
it1 = x(i)
x(i) = x(lpiv)
x(lpiv) = it1
i = ilx - 1
j = iux + 1
outer_up: do
in_up1: do
i = i + 1
xk = x(i)
if (xk >= piv) exit in_up1
end do in_up1
!
! Ensure finite termination for next loop
!
it1 = xk
x(i) = piv
in_up2:do
j = j - 1
xk = x(j)
if (xk <= piv) exit in_up2
end do in_up2
x(i) = it1
if (j > i) then
it1 = x(i)
x(i) = x(j)
x(j) = it1
else
exit outer_up
end if
end do outer_up
if (i == ilx) then
if (x(i) /= piv) then
write(0,*) 'Should never ever get here????!!!!'
stop
endif
i = i + 1
endif
n1 = (i-1)-ilx+1
n2 = iux-(i)+1
if (n1 > n2) then
if (n1 > ithrs) then
istp = istp + 1
istack(1,istp) = ilx
istack(2,istp) = i-1
else
call iisr_up(n1,x(ilx:i-1))
endif
if (n2 > ithrs) then
istp = istp + 1
istack(1,istp) = i
istack(2,istp) = iux
else
call iisr_up(n2,x(i:iux))
endif
else
if (n2 > ithrs) then
istp = istp + 1
istack(1,istp) = i
istack(2,istp) = iux
else
call iisr_up(n2,x(i:iux))
endif
if (n1 > ithrs) then
istp = istp + 1
istack(1,istp) = ilx
istack(2,istp) = i-1
else
call iisr_up(n1,x(ilx:i-1))
endif
endif
enddo
else
call iisr_up(n,x)
endif
case(psb_sort_down_)
if (n > ithrs) then
!
! Init stack pointer
!
istp = 1
istack(1,istp) = 1
istack(2,istp) = n
do
if (istp <= 0) exit
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
!
i = ilx
j = iux
lpiv = (i+j)/2
piv = x(lpiv)
if (piv > x(i)) then
it1 = x(i)
x(i) = x(lpiv)
x(lpiv) = it1
piv = x(lpiv)
endif
if (piv < x(j)) then
it1 = x(j)
x(j) = x(lpiv)
x(lpiv) = it1
piv = x(lpiv)
endif
if (piv > x(i)) then
it1 = x(i)
x(i) = x(lpiv)
x(lpiv) = it1
piv = x(lpiv)
endif
!
! now piv is correct; place it into first location
it1 = x(i)
x(i) = x(lpiv)
x(lpiv) = it1
i = ilx - 1
j = iux + 1
outer_dw: do
in_dw1: do
i = i + 1
xk = x(i)
if (xk <= piv) exit in_dw1
end do in_dw1
!
! Ensure finite termination for next loop
!
it1 = xk
x(i) = piv
in_dw2:do
j = j - 1
xk = x(j)
if (xk >= piv) exit in_dw2
end do in_dw2
x(i) = it1
if (j > i) then
it1 = x(i)
x(i) = x(j)
x(j) = it1
else
exit outer_dw
end if
end do outer_dw
if (i == ilx) then
if (x(i) /= piv) then
write(0,*) 'Should never ever get here????!!!!'
stop
endif
i = i + 1
endif
n1 = (i-1)-ilx+1
n2 = iux-(i)+1
if (n1 > n2) then
if (n1 > ithrs) then
istp = istp + 1
istack(1,istp) = ilx
istack(2,istp) = i-1
else
call iisr_dw(n1,x(ilx:i-1))
endif
if (n2 > ithrs) then
istp = istp + 1
istack(1,istp) = i
istack(2,istp) = iux
else
call iisr_dw(n2,x(i:iux))
endif
else
if (n2 > ithrs) then
istp = istp + 1
istack(1,istp) = i
istack(2,istp) = iux
else
call iisr_dw(n2,x(i:iux))
endif
if (n1 > ithrs) then
istp = istp + 1
istack(1,istp) = ilx
istack(2,istp) = i-1
else
call iisr_dw(n1,x(ilx:i-1))
endif
endif
enddo
else
call iisr_dw(n,x)
endif
case default
write(0,*) 'isr error !',dir
end select
return
contains
subroutine iisr_up(n,x)
integer :: n
integer :: x(n)
integer :: i,j
integer :: xx
do j=n-1,1,-1
if (x(j+1) < x(j)) then
xx = x(j)
i=j+1
do
x(i-1) = x(i)
i = i+1
if (i>n) exit
if (x(i) >= xx) exit
end do
x(i-1) = xx
endif
enddo
end subroutine iisr_up
subroutine iisr_dw(n,x)
integer :: n
integer :: x(n)
integer :: i,j
integer :: xx
do j=n-1,1,-1
if (x(j+1) > x(j)) then
xx = x(j)
i=j+1
do
x(i-1) = x(i)
i = i+1
if (i>n) exit
if (x(i) <= xx) exit
end do
x(i-1) = xx
endif
enddo
end subroutine iisr_dw
end subroutine isr

198
base/serial/aux/isrx.f
Unescape Escape View File

@ -1,198 +0,0 @@
C
C Parallel Sparse BLAS v2.0
C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
C Alfredo Buttari University of Rome Tor Vergata
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C 1. Redistributions of source code must retain the above copyright
C notice, this list of conditions and the following disclaimer.
C 2. Redistributions in binary form must reproduce the above copyright
C notice, this list of conditions, and the following disclaimer in the
C documentation and/or other materials provided with the distribution.
C 3. The name of the PSBLAS group or the names of its contributors may
C not be used to endorse or promote products derived from this
C software without specific written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
C POSSIBILITY OF SUCH DAMAGE.
C
C
SUBROUTINE ISRX(N,X,INDX)
C
C Quicksort with indices into original positions.
C Adapted from a number of sources, including Don Knuth's TAOCP.
C
C .. Scalar Arguments ..
INTEGER N
C ..
C .. Array Arguments ..
INTEGER INDX(N),X(N)
C ..
C .. Local Scalars ..
INTEGER I, J, II, XX, ILX, IUX, ISTP, PIV, LPIV
INTEGER IT1, IT2, N1, N2
INTEGER MAXSTACK,NPARMS,ITHRS
PARAMETER (MAXSTACK=64,NPARMS=3,ITHRS=16)
INTEGER ISTACK(NPARMS,MAXSTACK)
C ..
DO I=1, N
INDX(I) = I
ENDDO
C
C Small inputs will only get through insertion sort.
C
IF (N.GT.ITHRS) THEN
C
C Init stack pointer
C
ISTP = 1
ISTACK(1,ISTP) = 1
ISTACK(2,ISTP) = N
DO WHILE (ISTP.GT.0)
ILX = ISTACK(1,ISTP)
IUX = ISTACK(2,ISTP)
ISTP = ISTP - 1
C
C Choose a pivot with median-of-three heuristics, leave it
C in the LPIV location
C
I = ILX
J = IUX
LPIV = (I+J)/2
PIV = X(LPIV)
IF (PIV.LT.X(I)) THEN
IT1 = X(I)
IT2 = INDX(I)
X(I) = X(LPIV)
INDX(I) = INDX(LPIV)
X(LPIV) = IT1
INDX(LPIV) = IT2
PIV = X(LPIV)
ENDIF
IF (PIV.GT.X(J)) THEN
IT1 = X(J)
IT2 = INDX(J)
X(J) = X(LPIV)
INDX(J) = INDX(LPIV)
X(LPIV) = IT1
INDX(LPIV) = IT2
PIV = X(LPIV)
ENDIF
IF (PIV.LT.X(I)) THEN
IT1 = X(I)
IT2 = INDX(I)
X(I) = X(LPIV)
INDX(I) = INDX(LPIV)
X(LPIV) = IT1
INDX(LPIV) = IT2
PIV = X(LPIV)
ENDIF
C
C Now PIV is correct; place it into first location
C
IT1 = X(I)
IT2 = INDX(I)
X(I) = X(LPIV)
INDX(I) = INDX(LPIV)
X(LPIV) = IT1
INDX(LPIV) = IT2
I = ILX - 1
J = IUX + 1
130 CONTINUE
I = I + 1
XK = X(I)
IF (XK.LT.PIV) GOTO 130
C
C Ensure finite termination for next loop
C
IT1 = XK
X(I) = PIV
140 CONTINUE
J = J - 1
XK = X(J)
IF (XK.GT.PIV) GOTO 140
X(I) = IT1
150 CONTINUE
IF (J.GT.I) THEN
IT1 = X(I)
IT2 = INDX(I)
X(I) = X(J)
INDX(I) = INDX(J)
X(J) = IT1
INDX(J) = IT2
GO TO 130
END IF
if (i.eq.ilx) then
if (x(i).ne.piv) then
write(0,*)
+ 'ISRX:: Should never ever get here????!!!!'
stop
endif
i = i + 1
endif
N1 = (I-1)-ILX+1
N2 = IUX-(I)+1
IF (N1.GT.N2) THEN
if (n1.gt.ithrs) then
ISTP = ISTP + 1
ISTACK(1,ISTP) = ILX
ISTACK(2,ISTP) = I-1
endif
if (n2.gt.ithrs) then
ISTP = ISTP + 1
ISTACK(1,ISTP) = I
ISTACK(2,ISTP) = IUX
endif
ELSE
if (n2.gt.ithrs) then
ISTP = ISTP + 1
ISTACK(1,ISTP) = I
ISTACK(2,ISTP) = IUX
endif
if (n1.gt.ithrs) then
ISTP = ISTP + 1
ISTACK(1,ISTP) = ILX
ISTACK(2,ISTP) = I-1
endif
ENDIF
ENDDO
ENDIF
DO J=N-1,1,-1
IF (X(J+1).LT.X(J)) THEN
XX = X(J)
II = INDX(J)
I=J+1
100 CONTINUE
X(I-1) = X(I)
INDX(I-1) = INDX(I)
I = I+1
IF ((I.LE.N)) then
if (X(I).LT.XX) GOTO 100
endif
X(I-1) = XX
INDX(I-1) = II
ENDIF
ENDDO
RETURN
END

405
base/serial/aux/isrx.f90
Unescape Escape View File

@ -0,0 +1,405 @@
!!$
!!$ 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.
!!$
!!$
subroutine isrx(n,x,indx,dir,flag)
use psb_serial_mod
!
! Quicksort with indices into original positions.
! Adapted from a number of sources, including Don Knuth's TAOCP.
!
! .. Scalar Arguments ..
integer, intent(in) :: n, dir, flag
integer :: x(n), indx(n)
! ..
! .. Local Scalars ..
integer i, j, ii, xx, ilx, iux, istp, piv, lpiv
integer it1, it2, n1, n2
integer, parameter :: maxstack=64,nparms=3,ithrs=16
integer :: istack(nparms,maxstack)
! ..
select case(flag)
case(psb_sort_ovw_idx_)
do i=1, n
indx(i) = i
enddo
case(psb_sort_keep_idx_)
! do nothing
case default
write(0,*) 'Error in isrx: invalid flag',flag
end select
!
!
! small inputs will only get through insertion sort.
!
select case(dir)
case(psb_sort_up_)
if (n > ithrs) then
!
! Init stack pointer
!
istp = 1
istack(1,istp) = 1
istack(2,istp) = n
do
if (istp <= 0) exit
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
!
i = ilx
j = iux
lpiv = (i+j)/2
piv = x(lpiv)
if (piv < x(i)) then
it1 = x(i)
it2 = indx(i)
x(i) = x(lpiv)
indx(i) = indx(lpiv)
x(lpiv) = it1
indx(lpiv) = it2
piv = x(lpiv)
endif
if (piv > x(j)) then
it1 = x(j)
it2 = indx(j)
x(j) = x(lpiv)
indx(j) = indx(lpiv)
x(lpiv) = it1
indx(lpiv) = it2
piv = x(lpiv)
endif
if (piv < x(i)) then
it1 = x(i)
it2 = indx(i)
x(i) = x(lpiv)
indx(i) = indx(lpiv)
x(lpiv) = it1
indx(lpiv) = it2
piv = x(lpiv)
endif
!
! now piv is correct; place it into first location
it1 = x(i)
it2 = indx(i)
x(i) = x(lpiv)
indx(i) = indx(lpiv)
x(lpiv) = it1
indx(lpiv) = it2
piv = x(lpiv)
i = ilx - 1
j = iux + 1
outer_up: do
in_up1: do
i = i + 1
xk = x(i)
if (xk >= piv) exit in_up1
end do in_up1
!
! Ensure finite termination for next loop
!
it1 = xk
x(i) = piv
in_up2:do
j = j - 1
xk = x(j)
if (xk <= piv) exit in_up2
end do in_up2
x(i) = it1
if (j > i) then
it1 = x(i)
it2 = indx(i)
x(i) = x(j)
indx(i) = indx(j)
x(j) = it1
indx(j) = it2
else
exit outer_up
end if
end do outer_up
if (i == ilx) then
if (x(i) /= piv) then
write(0,*) 'Should never ever get here????!!!!'
stop
endif
i = i + 1
endif
n1 = (i-1)-ilx+1
n2 = iux-(i)+1
if (n1 > n2) then
if (n1 > ithrs) then
istp = istp + 1
istack(1,istp) = ilx
istack(2,istp) = i-1
else
call iisrx_up(n1,x(ilx:i-1),indx(ilx:i-1))
endif
if (n2 > ithrs) then
istp = istp + 1
istack(1,istp) = i
istack(2,istp) = iux
else
call iisrx_up(n2,x(i:iux),indx(i:iux))
endif
else
if (n2 > ithrs) then
istp = istp + 1
istack(1,istp) = i
istack(2,istp) = iux
else
call iisrx_up(n2,x(i:iux),indx(i:iux))
endif
if (n1 > ithrs) then
istp = istp + 1
istack(1,istp) = ilx
istack(2,istp) = i-1
else
call iisrx_up(n1,x(ilx:i-1),indx(ilx:i-1))
endif
endif
enddo
else
call iisrx_up(n,x,indx)
endif
case(psb_sort_down_)
if (n > ithrs) then
!
! Init stack pointer
!
istp = 1
istack(1,istp) = 1
istack(2,istp) = n
do
if (istp <= 0) exit
ilx = istack(1,istp)
iux = istack(2,istp)
istp = istp - 1
!$$$ write(0,*) 'Debug 1: ',ilx,iux
!
! Choose a pivot with median-of-three heuristics, leave it
! in the LPIV location
!
i = ilx
j = iux
lpiv = (i+j)/2
piv = x(lpiv)
if (piv > x(i)) then
it1 = x(i)
it2 = indx(i)
x(i) = x(lpiv)
indx(i) = indx(lpiv)
x(lpiv) = it1
indx(lpiv) = it2
piv = x(lpiv)
endif
if (piv < x(j)) then
it1 = x(j)
it2 = indx(j)
x(j) = x(lpiv)
indx(j) = indx(lpiv)
x(lpiv) = it1
indx(lpiv) = it2
piv = x(lpiv)
endif
if (piv > x(i)) then
it1 = x(i)
it2 = indx(i)
x(i) = x(lpiv)
indx(i) = indx(lpiv)
x(lpiv) = it1
indx(lpiv) = it2
piv = x(lpiv)
endif
!
! now piv is correct; place it into first location
it1 = x(i)
it2 = indx(i)
x(i) = x(lpiv)
indx(i) = indx(lpiv)
x(lpiv) = it1
indx(lpiv) = it2
piv = x(lpiv)
i = ilx - 1
j = iux + 1
outer_dw: do
in_dw1: do
i = i + 1
xk = x(i)
if (xk <= piv) exit in_dw1
end do in_dw1
!
! Ensure finite termination for next loop
!
it1 = xk
x(i) = piv
in_dw2:do
j = j - 1
xk = x(j)
if (xk >= piv) exit in_dw2
end do in_dw2
x(i) = it1
if (j > i) then
it1 = x(i)
it2 = indx(i)
x(i) = x(j)
indx(i) = indx(j)
x(j) = it1
indx(j) = it2
else
exit outer_dw
end if
end do outer_dw
if (i == ilx) then
if (x(i) /= piv) then
write(0,*) 'Should never ever get here????!!!!'
stop
endif
i = i + 1
endif
n1 = (i-1)-ilx+1
n2 = iux-(i)+1
if (n1 > n2) then
if (n1 > ithrs) then
istp = istp + 1
istack(1,istp) = ilx
istack(2,istp) = i-1
else
call iisrx_dw(n1,x(ilx:i-1),indx(ilx:i-1))
endif
if (n2 > ithrs) then
istp = istp + 1
istack(1,istp) = i
istack(2,istp) = iux
else
call iisrx_dw(n2,x(i:iux),indx(i:iux))
endif
else
if (n2 > ithrs) then
istp = istp + 1
istack(1,istp) = i
istack(2,istp) = iux
else
call iisrx_dw(n2,x(i:iux),indx(i:iux))
endif
if (n1 > ithrs) then
istp = istp + 1
istack(1,istp) = ilx
istack(2,istp) = i-1
else
call iisrx_dw(n1,x(ilx:i-1),indx(ilx:i-1))
endif
endif
enddo
else
call iisrx_dw(n,x,indx)
endif
case default
write(0,*) 'isrx error dir ',dir
end select
return
contains
subroutine iisrx_up(n,x,indx)
integer :: n
integer :: x(n)
integer :: indx(n)
integer :: i,j,ix
integer :: xx
do j=n-1,1,-1
if (x(j+1) < x(j)) then
xx = x(j)
ix = indx(j)
i=j+1
do
x(i-1) = x(i)
indx(i-1) = indx(i)
i = i+1
if (i>n) exit
if (x(i) >= xx) exit
end do
x(i-1) = xx
indx(i-1) = ix
endif
enddo
end subroutine iisrx_up
subroutine iisrx_dw(n,x,indx)
integer :: n
integer :: x(n)
integer :: indx(n)
integer :: i,j,ix
integer :: xx
do j=n-1,1,-1
if (x(j+1) > x(j)) then
xx = x(j)
ix = indx(j)
i=j+1
do
x(i-1) = x(i)
indx(i-1) = indx(i)
i = i+1
if (i>n) exit
if (x(i) <= xx) exit
end do
x(i-1) = xx
indx(i-1) = ix
endif
enddo
end subroutine iisrx_dw
end subroutine isrx

234
base/serial/aux/mrgsrtd.f
Unescape Escape View File

@ -0,0 +1,234 @@
C
C Parallel Sparse BLAS v2.0
C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
C Alfredo Buttari University of Rome Tor Vergata
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C 1. Redistributions of source code must retain the above copyright
C notice, this list of conditions and the following disclaimer.
C 2. Redistributions in binary form must reproduce the above copyright
C notice, this list of conditions, and the following disclaimer in the
C documentation and/or other materials provided with the distribution.
C 3. The name of the PSBLAS group or the names of its contributors may
C not be used to endorse or promote products derived from this
C software without specific written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
C POSSIBILITY OF SUCH DAMAGE.
C
C
***********************************************************************
* *
* FUNCTION = This subroutine returns an array of pointers, L, *
* to be used to sort the integer input vector K; *
* the routine implements a list merge-sort *
* *
***********************************************************************
* *
* CALL MRGSRTD(N,K,L,IRET) *
* *
* INPUT = *
* *
* SYMBOLIC NAME: N *
* POSITION: First parameter. *
* ATTRIBUTES: INTEGER *
* VALUES: >= 0 *
* DESCRIPTION: Dimension of the array to be sorted *
* *
* SYMBOLIC NAME: K *
* POSITION: Second parameter *
* ATTRIBUTES: INTEGER ARRAY(N) *
* VALUES: Any *
* DESCRIPTION: Input array containing the keys, i.e., values *
* to be sorted *
* *
* *
* *
* OUTPUT = *
* *
* SYMBOLIC NAME: L *
* POSITION: Third parameter *
* ATTRIBUTES: INTEGER ARRAY(N+2) *
* VALUES: >= 0 *
* DESCRIPTION: On exit, this array contains pointers to *
* the keys array. *
* *
***********************************************************************
***********************************************************************
* *
***********************************************************************
***********************************************************************
* ALGORITHM DESCRIPTION *
* *
* REFERENCES = (1) D. E. Knuth *
* The Art of Computer Programming, *
* vol.3: Sorting and Searching *
* Addison-Wesley, 1973 *
* *
* FUNCTION = This subroutine is based on the well-known merge-sort *
* algorithm; according to (1) we are sorting 'records' *
* R(I) with respect to keys K(I), and to this purpose *
* we use 'links' L(I); at the end of the subroutine, *
* L(0) is the index of the first record in the sorted *
* sequence, then for every record R(I), we have into *
* L(I) the index of the next one in the sequence. A *
* value L(I)=0 signals the end of the sequence. *
* The sorting is stable, i.e., if K(I)=K(J) and I<J, *
* then in the sorted sequence R(I) precedes R(J); many *
* sorting algorithms, e.g. quicksort, are not stable. *
* The list merge-sort is one of the fastest stable *
* sortings available; it is guaranteed to run in *
* O(N log N) time on both the average and worst cases. *
* *
* *
***********************************************************************
***********************************************************************
* ALGORITHM EXAMPLE(S) *
* *
* EXAMPLE: Construct a sorted array of records RS from a vector R *
* according to the keys stored in K *
* *
* CALL MRGSRTD(N,K,L,*100) *
* I = L(0) *
* DO 100 J = 1, N *
* RS(J) = R(I) *
* I = L(I) *
* 100 CONTINUE ! RETURN POINT IF ARRAY ALREADY SORTED *
* *
* *
* EXAMPLE: Sort in place array R *
* *
* CALL MRGSRTD(N,K,L,*400) *
* LP = L(0) *
* KK = 1 *
* 100 CONTINUE *
* IF ((LP.EQ.0).OR.(KK.GT.N)) GOTO 400 *
* 200 CONTINUE *
* IF (LP.GE.KK) GOTO 300 *
* LP = L(LP) *
* GOTO 200 *
* 300 CONTINUE *
* SWAP = R(KK) *
* R(KK) = R(LP) *
* R(LP) = SWAP *
* LSWAP = L(LP) *
* L(LP) = L(KK) *
* L(KK) = LP *
* LP = LSWAP *
* KK = KK+1 *
* GOTO 100 *
* 400 CONTINUE *
* *
* *
***********************************************************************
SUBROUTINE MRGSRTD(N,K,L,IRET)
C .. Scalar Arguments ..
INTEGER N, IRET
C ..
C .. Array Arguments ..
INTEGER K(N),L(0:N+1)
C ..
C .. Local Scalars ..
INTEGER P,Q,S,T
C ..
C .. Intrinsic Functions ..
INTRINSIC IABS,ISIGN
C ..
IRET = 0
C First step: we are preparing ordered sublists, exploiting
C what order was already in the input data; negative links
C mark the end of the sublists
L(0) = 1
T = N + 1
DO P = 1,N - 1
IF (K(P).GE.K(P+1)) THEN
L(P) = P + 1
ELSE
L(T) = - (P+1)
T = P
END IF
END DO
L(T) = 0
L(N) = 0
C See if the input was already sorted
IF (L(N+1).EQ.0) THEN
IRET = 1
RETURN
ELSE
L(N+1) = IABS(L(N+1))
END IF
200 CONTINUE
C Otherwise, begin a pass through the list.
C Throughout all the subroutine we have:
C P, Q: pointing to the sublists being merged
C S: pointing to the most recently processed record
C T: pointing to the end of previously completed sublist
S = 0
T = N + 1
P = L(S)
Q = L(T)
IF (Q.EQ.0) RETURN
300 CONTINUE
IF (K(P).LT.K(Q)) GO TO 600
400 CONTINUE
L(S) = ISIGN(P,L(S))
S = P
P = L(P)
IF (P.GT.0) GO TO 3100
C Otherwise, one sublist ended, and we append to it the rest
C of the other one.
500 CONTINUE
L(S) = Q
S = T
550 CONTINUE
T = Q
Q = L(Q)
IF (Q.GT.0) GO TO 550
GO TO 800
600 CONTINUE
L(S) = ISIGN(Q,L(S))
S = Q
Q = L(Q)
IF (Q.GT.0) GO TO 3200
700 CONTINUE
L(S) = P
S = T
750 CONTINUE
T = P
P = L(P)
IF (P.GT.0) GO TO 750
800 CONTINUE
P = -P
Q = -Q
IF (Q.EQ.0) THEN
L(S) = ISIGN(P,L(S))
L(T) = 0
GO TO 200
ELSE
GO TO 300
END IF
3100 CONTINUE
IF (K(P).LT.K(Q)) GO TO 600
S = P
P = L(P)
IF (P.GT.0) GO TO 3100
GO TO 500
3200 CONTINUE
IF (K(P).GE.K(Q)) GO TO 400
S = Q
Q = L(Q)
IF (Q.GT.0) GO TO 3200
GO TO 700
END

3
base/serial/f77/smmp.f
Unescape Escape View File

@ -17,6 +17,7 @@ c=======================================================================
* ic, jc, diagc, * ic, jc, diagc,
* index) * index)
use psb_realloc_mod use psb_realloc_mod
use psb_serial_mod, only: psb_msort
c c
integer ia(*), ja(*), diaga, integer ia(*), ja(*), diaga,
* ib(*), jb(*), diagb, * ib(*), jb(*), diagb,
@ -120,7 +121,7 @@ c
istart=index(istart) istart=index(istart)
index(jc(j))=0 index(jc(j))=0
40 continue 40 continue
call isr(length,jc(ic(i))) call psb_msort(jc(ic(i):ic(i)+length -1))
c$$$ write(iunit,*) length,' : ',jc(ic(i):ic(i)+length-1) c$$$ write(iunit,*) length,' : ',jc(ic(i):ic(i)+length-1)
index(i) = 0 index(i) = 0
50 continue 50 continue

30
base/serial/lsame.f90
Unescape Escape View File

@ -1,3 +1,33 @@
!!$
!!$ 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.
!!$
!!$
function lsame(a,b) function lsame(a,b)
use psb_string_mod use psb_string_mod
logical :: lsame logical :: lsame

3
base/serial/psb_dsymbmm.f90
Unescape Escape View File

@ -41,6 +41,7 @@
subroutine psb_dsymbmm(a,b,c,info) subroutine psb_dsymbmm(a,b,c,info)
use psb_spmat_type use psb_spmat_type
use psb_string_mod use psb_string_mod
use psb_serial_mod, only : psb_msort
implicit none implicit none
type(psb_dspmat_type) :: a,b,c type(psb_dspmat_type) :: a,b,c
@ -209,7 +210,7 @@ contains
istart=index(istart) istart=index(istart)
index(c%ia1(j))=0 index(c%ia1(j))=0
end do end do
call isr(length,c%ia1(c%ia2(i))) call psb_msort(c%ia1(c%ia2(i):c%ia2(i)+length-1))
index(i) = 0 index(i) = 0
end do main end do main

2
base/serial/psb_zsymbmm.f90
Unescape Escape View File

@ -204,7 +204,7 @@ contains
istart=index(istart) istart=index(istart)
index(c%ia1(j))=0 index(c%ia1(j))=0
end do end do
call isr(length,c%ia1(c%ia2(i))) call psb_msort(c%ia1(c%ia2(i):c%ia2(i)+length-1))
index(i) = 0 index(i) = 0
end do main end do main

2
base/tools/psb_dcdovr.f90
Unescape Escape View File

@ -436,7 +436,7 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype)
if (i_ovr <= novr) then if (i_ovr <= novr) then
if (tot_elem > 1) then if (tot_elem > 1) then
call imsr(tot_elem,works(idxs+1)) call psb_msort(works(idxs+1:idxs+tot_elem))
lx = works(idxs+1) lx = works(idxs+1)
i = 1 i = 1
j = 1 j = 1

2
base/tools/psb_zcdovr.f90
Unescape Escape View File

@ -436,7 +436,7 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype)
if (i_ovr <= novr) then if (i_ovr <= novr) then
if (tot_elem > 1) then if (tot_elem > 1) then
call imsr(tot_elem,works(idxs+1)) call psb_msort(works(idxs+1:idxs+tot_elem))
lx = works(idxs+1) lx = works(idxs+1)
i = 1 i = 1
j = 1 j = 1

6
baseprec/psb_dsp_renum.f90
Unescape Escape View File

@ -91,7 +91,7 @@ subroutine psb_dsp_renum(a,desc_a,p,atmp,info)
! !
! We want: NEW(I) = OLD(PERM(I)) ! We want: NEW(I) = OLD(PERM(I))
! !
call isrx(nnr,itmp2,p%perm) call psb_msort(itmp2(1:nnr),ix=p%perm)
do k=1, nnr do k=1, nnr
p%invperm(p%perm(k)) = k p%invperm(p%perm(k)) = k
@ -132,7 +132,7 @@ subroutine psb_dsp_renum(a,desc_a,p,atmp,info)
atmp%ia1(j+k-1) = p%invperm(a%ia1(jj+kk-1)) atmp%ia1(j+k-1) = p%invperm(a%ia1(jj+kk-1))
endif endif
enddo enddo
call isrx(k,atmp%ia1(j:j+k-1),itmp2) call psb_msort(atmp%ia1(j:j+k-1),ix=itmp2)
do kk=1,k do kk=1,k
atmp%aspk(j+kk-1) = rtmp(itmp2(kk)) atmp%aspk(j+kk-1) = rtmp(itmp2(kk))
enddo enddo
@ -256,7 +256,7 @@ subroutine psb_dsp_renum(a,desc_a,p,atmp,info)
atmp%ia1(j+k-1) = p%invperm(a%ia1(jj+kk-1)) atmp%ia1(j+k-1) = p%invperm(a%ia1(jj+kk-1))
endif endif
enddo enddo
call isrx(k,atmp%ia1(j:j+k-1),itmp2) call psb_msort(atmp%ia1(j:j+k-1),ix=itmp2)
do kk=1,k do kk=1,k
atmp%aspk(j+kk-1) = rtmp(itmp2(kk)) atmp%aspk(j+kk-1) = rtmp(itmp2(kk))
enddo enddo

6
baseprec/psb_zsp_renum.f90
Unescape Escape View File

@ -89,7 +89,7 @@ subroutine psb_zsp_renum(a,desc_a,p,atmp,info)
! !
! We want: NEW(I) = OLD(PERM(I)) ! We want: NEW(I) = OLD(PERM(I))
! !
call isrx(nnr,itmp2,p%perm) call psb_msort(itmp2(1:nnr),ix=p%perm)
do k=1, nnr do k=1, nnr
p%invperm(p%perm(k)) = k p%invperm(p%perm(k)) = k
@ -130,7 +130,7 @@ subroutine psb_zsp_renum(a,desc_a,p,atmp,info)
atmp%ia1(j+k-1) = p%invperm(a%ia1(jj+kk-1)) atmp%ia1(j+k-1) = p%invperm(a%ia1(jj+kk-1))
endif endif
enddo enddo
call isrx(k,atmp%ia1(j:j+k-1),itmp2) call psb_msort(atmp%ia1(j:j+k-1),ix=itmp2)
do kk=1,k do kk=1,k
atmp%aspk(j+kk-1) = ztmp(itmp2(kk)) atmp%aspk(j+kk-1) = ztmp(itmp2(kk))
enddo enddo
@ -254,7 +254,7 @@ subroutine psb_zsp_renum(a,desc_a,p,atmp,info)
atmp%ia1(j+k-1) = p%invperm(a%ia1(jj+kk-1)) atmp%ia1(j+k-1) = p%invperm(a%ia1(jj+kk-1))
endif endif
enddo enddo
call isrx(k,atmp%ia1(j:j+k-1),itmp2) call psb_msort(atmp%ia1(j:j+k-1),ix=itmp2)
do kk=1,k do kk=1,k
atmp%aspk(j+kk-1) = ztmp(itmp2(kk)) atmp%aspk(j+kk-1) = ztmp(itmp2(kk))
enddo enddo

Write Preview
Loading…
Cancel
Save