diff --git a/base/internals/psi_crea_bnd_elem.f90 b/base/internals/psi_crea_bnd_elem.f90 index 2a79ab8a..cb5355ef 100644 --- a/base/internals/psi_crea_bnd_elem.f90 +++ b/base/internals/psi_crea_bnd_elem.f90 @@ -32,6 +32,7 @@ subroutine psi_crea_bnd_elem(bndel,desc_a,info) use psb_realloc_mod use psb_descriptor_type use psb_error_mod + use psb_serial_mod implicit none integer, allocatable :: bndel(:) @@ -67,7 +68,7 @@ subroutine psi_crea_bnd_elem(bndel,desc_a,info) enddo if (i>0) then - call isr(i,work) + call psb_msort(work(1:i)) j=1 irv = work(1) do k=2, i diff --git a/base/internals/psi_crea_ovr_elem.f90 b/base/internals/psi_crea_ovr_elem.f90 index 28930043..003b4f50 100644 --- a/base/internals/psi_crea_ovr_elem.f90 +++ b/base/internals/psi_crea_ovr_elem.f90 @@ -33,6 +33,7 @@ subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem,info) use psb_realloc_mod use psb_error_mod use psb_penv_mod + use psb_serial_mod implicit none ! ...parameter arrays.... @@ -154,7 +155,7 @@ subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem,info) i=i+2*desc_overlap(i)+2 enddo if (nel > 0) then - call imsr(nel,telem(:,1)) + call psb_msort(telem(1:nel,1)) iel = telem(1,1) telem(1,2) = 2 ix = 1 diff --git a/base/internals/psi_ldsc_pre_halo.f90 b/base/internals/psi_ldsc_pre_halo.f90 index af774438..57b98275 100644 --- a/base/internals/psi_ldsc_pre_halo.f90 +++ b/base/internals/psi_ldsc_pre_halo.f90 @@ -113,7 +113,8 @@ subroutine psi_ldsc_pre_halo(desc,ext_hv,info) idx = desc%hashv(i) nh = desc%hashv(i+1) - desc%hashv(i) 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 do diff --git a/base/internals/srtlist.f b/base/internals/srtlist.f index 83cf43ce..0000c27f 100644 --- a/base/internals/srtlist.f +++ b/base/internals/srtlist.f @@ -78,6 +78,7 @@ C *********************************************************************** SUBROUTINE SRTLIST(DEP_LIST,DL_LDA,LDL,NP,dg,dgp,upd, + edges,idx,ich,INFO) + use psb_serial_mod IMPLICIT NONE INTEGER NP, DL_LDA, INFO 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))) ENDDO - CALL ISRX(NEDGES-IST+1,DGP(IST),IDX(IST)) + call psb_msort(dgp(ist:nedges),ix=idx(ist:nedges)) I1 = IST NCH = 0 DO I = IST, NEDGES @@ -156,7 +157,7 @@ c$$$ write(0,*) 'SRTLIST Input :',i,ip info = 30 return ENDIF - CALL ISR(NCH,ICH) + call psb_msort(ich(1:nch)) DO I=1, NCH ISWAP(1) = EDGES(1,IST) ISWAP(2) = EDGES(2,IST) diff --git a/base/modules/Makefile b/base/modules/Makefile index 2cbf591c..e39fd6cb 100644 --- a/base/modules/Makefile +++ b/base/modules/Makefile @@ -6,9 +6,9 @@ MODULES = psb_realloc_mod.o psb_string_mod.o psb_spmat_type.o \ psb_error_mod.o \ psb_const_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) MPFOBJS=psb_penv_mod.o diff --git a/base/modules/blacs_env.F90 b/base/modules/blacs_env.F90 deleted file mode 100644 index 8aae4ad2..00000000 --- a/base/modules/blacs_env.F90 +++ /dev/null @@ -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 diff --git a/base/modules/psb_penv_mod.F90 b/base/modules/psb_penv_mod.F90 index 3ad735ee..abd9f05f 100644 --- a/base/modules/psb_penv_mod.F90 +++ b/base/modules/psb_penv_mod.F90 @@ -71,7 +71,7 @@ module psb_penv_mod module procedure psb_ibcasts, psb_ibcastv, psb_ibcastm,& & psb_dbcasts, psb_dbcastv, psb_dbcastm,& & psb_zbcasts, psb_zbcastv, psb_zbcastm,& - & psb_hbcasts, psb_lbcasts + & psb_hbcasts, psb_lbcasts, psb_lbcastv end interface @@ -544,6 +544,29 @@ contains 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) use mpi @@ -4785,6 +4808,69 @@ contains end subroutine zgamn2dm #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 diff --git a/base/modules/psb_serial_mod.f90 b/base/modules/psb_serial_mod.f90 index 6218e35d..bfd8de51 100644 --- a/base/modules/psb_serial_mod.f90 +++ b/base/modules/psb_serial_mod.f90 @@ -419,8 +419,136 @@ module psb_serial_mod end subroutine psb_zspgetrow 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 + - end module psb_serial_mod diff --git a/base/serial/aux/Makefile b/base/serial/aux/Makefile index 859a1f52..98807995 100644 --- a/base/serial/aux/Makefile +++ b/base/serial/aux/Makefile @@ -4,7 +4,7 @@ include ../../../Make.inc # 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) diff --git a/base/serial/aux/imsr.f90 b/base/serial/aux/imsr.f90 index c183ece9..7032740f 100644 --- a/base/serial/aux/imsr.f90 +++ b/base/serial/aux/imsr.f90 @@ -31,10 +31,12 @@ ! File: imsr.f90 ! Subroutine: ! Parameters: -subroutine imsr(n,x) - integer :: n +subroutine imsr(n,x,idir) + use psb_serial_mod + integer :: n, idir integer :: x(n) + integer, allocatable :: iaux(:) integer :: iswap, iret, info, lp, k @@ -53,8 +55,11 @@ subroutine imsr(n,x) return endif - - call mrgsrt(n,x,iaux,iret) + if (idir==psb_sort_up_) then + call mrgsrt(n,x,iaux,iret) + else + call mrgsrtd(n,x,iaux,iret) + end if if (iret == 0) then lp = iaux(0) diff --git a/base/serial/aux/imsrx.f90 b/base/serial/aux/imsrx.f90 index 8bc7d15a..8c7ee072 100644 --- a/base/serial/aux/imsrx.f90 +++ b/base/serial/aux/imsrx.f90 @@ -31,8 +31,9 @@ ! File: imsrx.f90 ! Subroutine: ! Parameters: -subroutine imsrx(n,x,indx,flag) - integer :: n, flag +subroutine imsrx(n,x,indx,idir,flag) + use psb_serial_mod + integer :: n,idir,flag integer :: x(n) integer :: indx(n) @@ -47,12 +48,14 @@ subroutine imsrx(n,x,indx,flag) endif if (n==0) return - if (n==1) then - if (flag == 0) then - indx(1) = 1 - endif - return - endif + + if (flag == psb_sort_ovw_idx_) then + do k=1,n + indx(k) = k + enddo + end if + + if (n==1) return allocate(iaux(0:n+1),stat=info) if (info/=0) then @@ -60,12 +63,11 @@ subroutine imsrx(n,x,indx,flag) return endif - if (flag == 0) then - do k=1,n - indx(k) = k - enddo + if (idir == psb_sort_up_) then + call mrgsrt(n,x,iaux,iret) + else + call mrgsrtd(n,x,iaux,iret) end if - call mrgsrt(n,x,iaux,iret) if (iret /= 1) then lp = iaux(0) diff --git a/base/serial/aux/isr.f b/base/serial/aux/isr.f deleted file mode 100644 index c318ae3c..00000000 --- a/base/serial/aux/isr.f +++ /dev/null @@ -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 diff --git a/base/serial/aux/isr.f90 b/base/serial/aux/isr.f90 new file mode 100644 index 00000000..7b421eb2 --- /dev/null +++ b/base/serial/aux/isr.f90 @@ -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 diff --git a/base/serial/aux/isrx.f b/base/serial/aux/isrx.f deleted file mode 100644 index e59ea7df..00000000 --- a/base/serial/aux/isrx.f +++ /dev/null @@ -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 diff --git a/base/serial/aux/isrx.f90 b/base/serial/aux/isrx.f90 new file mode 100644 index 00000000..df48e338 --- /dev/null +++ b/base/serial/aux/isrx.f90 @@ -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 diff --git a/base/serial/aux/mrgsrtd.f b/base/serial/aux/mrgsrtd.f new file mode 100644 index 00000000..b6064409 --- /dev/null +++ b/base/serial/aux/mrgsrtd.f @@ -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 1) then - call imsr(tot_elem,works(idxs+1)) + call psb_msort(works(idxs+1:idxs+tot_elem)) lx = works(idxs+1) i = 1 j = 1 diff --git a/base/tools/psb_zcdovr.f90 b/base/tools/psb_zcdovr.f90 index 9f6c5fa6..010a42db 100644 --- a/base/tools/psb_zcdovr.f90 +++ b/base/tools/psb_zcdovr.f90 @@ -436,7 +436,7 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) if (i_ovr <= novr) 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) i = 1 j = 1 diff --git a/baseprec/psb_dsp_renum.f90 b/baseprec/psb_dsp_renum.f90 index e038ac61..c51bee39 100644 --- a/baseprec/psb_dsp_renum.f90 +++ b/baseprec/psb_dsp_renum.f90 @@ -91,7 +91,7 @@ subroutine psb_dsp_renum(a,desc_a,p,atmp,info) ! ! 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 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)) endif 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 atmp%aspk(j+kk-1) = rtmp(itmp2(kk)) 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)) endif 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 atmp%aspk(j+kk-1) = rtmp(itmp2(kk)) enddo diff --git a/baseprec/psb_zsp_renum.f90 b/baseprec/psb_zsp_renum.f90 index 76107607..5ee89564 100644 --- a/baseprec/psb_zsp_renum.f90 +++ b/baseprec/psb_zsp_renum.f90 @@ -89,7 +89,7 @@ subroutine psb_zsp_renum(a,desc_a,p,atmp,info) ! ! 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 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)) endif 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 atmp%aspk(j+kk-1) = ztmp(itmp2(kk)) 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)) endif 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 atmp%aspk(j+kk-1) = ztmp(itmp2(kk)) enddo