Reimplemented merge-sort inner routine.

psblas3-type-indexed
Salvatore Filippone 18 years ago
parent 2968c06cfb
commit 7cb47779bd

@ -3,8 +3,8 @@ include ../../../Make.inc
# The object files # The object files
# #
FOBJS = isr.o isrx.o \ FOBJS = isr.o isrx.o msort_up.o msort_dw.o\
mrgsrt.o mrgsrtd.o isaperm.o ibsrch.o issrch.o imsr.o imsrx.o isaperm.o ibsrch.o issrch.o imsr.o imsrx.o
OBJS=$(FOBJS) OBJS=$(FOBJS)

@ -56,9 +56,9 @@ subroutine imsr(n,x,idir)
endif endif
if (idir==psb_sort_up_) then if (idir==psb_sort_up_) then
call mrgsrt(n,x,iaux,iret) call msort_up(n,x,iaux,iret)
else else
call mrgsrtd(n,x,iaux,iret) call msort_dw(n,x,iaux,iret)
end if end if
if (iret == 0) then if (iret == 0) then

@ -64,9 +64,9 @@ subroutine imsrx(n,x,indx,idir,flag)
endif endif
if (idir == psb_sort_up_) then if (idir == psb_sort_up_) then
call mrgsrt(n,x,iaux,iret) call msort_up(n,x,iaux,iret)
else else
call mrgsrtd(n,x,iaux,iret) call msort_dw(n,x,iaux,iret)
end if end if
if (iret /= 1) then if (iret /= 1) then

@ -1,234 +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
***********************************************************************
* *
* 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 MRGSRT(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 MRGSRT(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 MRGSRT(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 MRGSRT(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).LE.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).GT.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).GT.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).LE.K(Q)) GO TO 400
S = Q
Q = L(Q)
IF (Q.GT.0) GO TO 3200
GO TO 700
END

@ -1,234 +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
***********************************************************************
* *
* 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

@ -0,0 +1,170 @@
!
! 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: msort_dw.f90
!
! Subroutine: msort_dw
! This subroutine sorts an integer array into ascending order.
!
! Parameters:
! n - integer Input: size of the array
! k - integer(*) input: array of keys to be sorted
! l - integer(0:n+1) output: link list
! iret - integer output: 0 Normal termination
! 1 the array was already sorted
! *
! REFERENCES = (1) D. E. Knuth *
! The Art of Computer Programming, *
! vol.3: Sorting and Searching *
! Addison-Wesley, 1973 *
! *
! call msort_dw(n,x,iaux,iret)
!
! if (iret == 0) then
! lp = iaux(0)
! k = 1
! do
! if ((lp==0).or.(k>n)) exit
! do
! if (lp >= k) exit
! lp = iaux(lp)
! end do
! iswap = x(lp)
! x(lp) = x(k)
! x(k) = iswap
! lswap = iaux(lp)
! iaux(lp) = iaux(k)
! iaux(k) = lp
! lp = lswap
! k = k + 1
! enddo
! end if
!
!
subroutine msort_dw(n,k,l,iret)
integer n, iret
integer k(n),l(0:n+1)
!
integer p,q,s,t
intrinsic iabs,isign
! ..
iret = 0
! first step: we are preparing ordered sublists, exploiting
! what order was already in the input data; negative links
! mark the end of the sublists
l(0) = 1
t = n + 1
do p = 1,n - 1
if (k(p) >= 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
! see if the input was already sorted
if (l(n+1) == 0) then
iret = 1
return
else
l(n+1) = iabs(l(n+1))
end if
mergepass: do
! otherwise, begin a pass through the list.
! throughout all the subroutine we have:
! p, q: pointing to the sublists being merged
! s: pointing to the most recently processed record
! t: pointing to the end of previously completed sublist
s = 0
t = n + 1
p = l(s)
q = l(t)
if (q == 0) exit mergepass
outer: do
if (k(p) < k(q)) then
l(s) = isign(q,l(s))
s = q
q = l(q)
if (q > 0) then
do
if (k(p) >= k(q)) cycle outer
s = q
q = l(q)
if (q <= 0) exit
end do
end if
l(s) = p
s = t
do
t = p
p = l(p)
if (p <= 0) exit
end do
else
l(s) = isign(p,l(s))
s = p
p = l(p)
if (p>0) then
do
if (k(p) < k(q)) cycle outer
s = p
p = l(p)
if (p <= 0) exit
end do
end if
! otherwise, one sublist ended, and we append to it the rest
! of the other one.
l(s) = q
s = t
do
t = q
q = l(q)
if (q <= 0) exit
end do
end if
p = -p
q = -q
if (q == 0) then
l(s) = isign(p,l(s))
l(t) = 0
exit outer
end if
end do outer
end do mergepass
end subroutine msort_dw

@ -0,0 +1,170 @@
!
! 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: msort_up.f90
!
! Subroutine: msort_up
! This subroutine sorts an integer array into ascending order.
!
! Parameters:
! n - integer Input: size of the array
! k - integer(*) input: array of keys to be sorted
! l - integer(0:n+1) output: link list
! iret - integer output: 0 Normal termination
! 1 the array was already sorted
! *
! REFERENCES = (1) D. E. Knuth *
! The Art of Computer Programming, *
! vol.3: Sorting and Searching *
! Addison-Wesley, 1973 *
! *
! call msort_up(n,x,iaux,iret)
!
! if (iret == 0) then
! lp = iaux(0)
! k = 1
! do
! if ((lp==0).or.(k>n)) exit
! do
! if (lp >= k) exit
! lp = iaux(lp)
! end do
! iswap = x(lp)
! x(lp) = x(k)
! x(k) = iswap
! lswap = iaux(lp)
! iaux(lp) = iaux(k)
! iaux(k) = lp
! lp = lswap
! k = k + 1
! enddo
! end if
!
!
subroutine msort_up(n,k,l,iret)
integer n, iret
integer k(n),l(0:n+1)
!
integer p,q,s,t
intrinsic iabs,isign
! ..
iret = 0
! first step: we are preparing ordered sublists, exploiting
! what order was already in the input data; negative links
! mark the end of the sublists
l(0) = 1
t = n + 1
do p = 1,n - 1
if (k(p) <= 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
! see if the input was already sorted
if (l(n+1) == 0) then
iret = 1
return
else
l(n+1) = iabs(l(n+1))
end if
mergepass: do
! otherwise, begin a pass through the list.
! throughout all the subroutine we have:
! p, q: pointing to the sublists being merged
! s: pointing to the most recently processed record
! t: pointing to the end of previously completed sublist
s = 0
t = n + 1
p = l(s)
q = l(t)
if (q == 0) exit mergepass
outer: do
if (k(p) > k(q)) then
l(s) = isign(q,l(s))
s = q
q = l(q)
if (q > 0) then
do
if (k(p) <= k(q)) cycle outer
s = q
q = l(q)
if (q <= 0) exit
end do
end if
l(s) = p
s = t
do
t = p
p = l(p)
if (p <= 0) exit
end do
else
l(s) = isign(p,l(s))
s = p
p = l(p)
if (p>0) then
do
if (k(p) > k(q)) cycle outer
s = p
p = l(p)
if (p <= 0) exit
end do
end if
! otherwise, one sublist ended, and we append to it the rest
! of the other one.
l(s) = q
s = t
do
t = q
q = l(q)
if (q <= 0) exit
end do
end if
p = -p
q = -q
if (q == 0) then
l(s) = isign(p,l(s))
l(t) = 0
exit outer
end if
end do outer
end do mergepass
end subroutine msort_up

@ -147,7 +147,7 @@ c
if (debug) write(0,*) 'build check :',ia2n(ip1+psb_nnzt_) if (debug) write(0,*) 'build check :',ia2n(ip1+psb_nnzt_)
c .... order with key ia1n ... c .... order with key ia1n ...
call mrgsrt(nnz,ia1n,aux,iret) call msort_up(nnz,ia1n,aux,iret)
if (iret.eq.0) call reordvn3(nnz,arn,ia1n,ia2n,aux(ipx),aux) if (iret.eq.0) call reordvn3(nnz,arn,ia1n,ia2n,aux(ipx),aux)
c .... order with key ia2n ... c .... order with key ia2n ...
@ -159,7 +159,7 @@ c .... order with key ia2n ...
j = j+1 j = j+1
enddo enddo
nzl = j - i nzl = j - i
call mrgsrt(nzl,ia2n(i),aux,iret) call msort_up(nzl,ia2n(i),aux,iret)
if (iret.eq.0) call reordvn3(nzl,arn(i),ia1n(i),ia2n(i), if (iret.eq.0) call reordvn3(nzl,arn(i),ia1n(i),ia2n(i),
+ aux(ipx+i-1),aux) + aux(ipx+i-1),aux)
i = j i = j
@ -198,7 +198,7 @@ c ... add the duplicated element ...
else else
c .... order with key ia1n ... c .... order with key ia1n ...
call mrgsrt(nnz,ia1n,aux,iret) call msort_up(nnz,ia1n,aux,iret)
if (iret.eq.0) call reordvn(nnz,arn,ia1n,ia2n,aux) if (iret.eq.0) call reordvn(nnz,arn,ia1n,ia2n,aux)
c .... order with key ia2n ... c .... order with key ia2n ...
@ -210,7 +210,7 @@ c .... order with key ia2n ...
j = j+1 j = j+1
enddo enddo
nzl = j - i nzl = j - i
call mrgsrt(nzl,ia2n(i),aux,iret) call msort_up(nzl,ia2n(i),aux,iret)
if (iret.eq.0) call reordvn(nzl,arn(i),ia1n(i),ia2n(i), if (iret.eq.0) call reordvn(nzl,arn(i),ia1n(i),ia2n(i),
+ aux) + aux)
i = j i = j

@ -157,7 +157,7 @@ c$$$ + ip1,ip2,nnz,ian2(ip1+nnzt_)
if (debug) write(0,*) 'Build check :',ian2(ip1+psb_nnzt_) if (debug) write(0,*) 'Build check :',ian2(ip1+psb_nnzt_)
C .... Order with key IA ... C .... Order with key IA ...
CALL MRGSRT(NNZ,IA,AUX,IRET) CALL MSORT_UP(NNZ,IA,AUX,IRET)
IF (IRET.EQ.0) CALL REORDVN3(NNZ,AR,IA,JA,AUX(IPX),AUX) IF (IRET.EQ.0) CALL REORDVN3(NNZ,AR,IA,JA,AUX(IPX),AUX)
if (debug) then if (debug) then
do i=1, nnz-1 do i=1, nnz-1
@ -183,7 +183,7 @@ c$$$ + (J.LE.NNZ))
J = J+1 J = J+1
ENDDO ENDDO
NZL = J - I NZL = J - I
CALL MRGSRT(NZL,JA(I),AUX,IRET) CALL MSORT_UP(NZL,JA(I),AUX,IRET)
IF (IRET.EQ.0) CALL REORDVN3(NZL,AR(I),IA(I),JA(I), IF (IRET.EQ.0) CALL REORDVN3(NZL,AR(I),IA(I),JA(I),
+ AUX(IPX+I-1),AUX) + AUX(IPX+I-1),AUX)
I = J I = J
@ -250,7 +250,7 @@ C ... Sum the duplicated element ...
ELSE ELSE
c$$$ write(0,*) 'Going for ELSE !!!?' c$$$ write(0,*) 'Going for ELSE !!!?'
C .... Order with key IA ... C .... Order with key IA ...
CALL MRGSRT(NNZ,IA,AUX,IRET) CALL MSORT_UP(NNZ,IA,AUX,IRET)
IF (IRET.EQ.0) CALL REORDVN(NNZ,AR,IA,JA,AUX) IF (IRET.EQ.0) CALL REORDVN(NNZ,AR,IA,JA,AUX)
C .... Order with key IA2N ... C .... Order with key IA2N ...
I = 1 I = 1
@ -266,7 +266,7 @@ c$$$ + (J.LE.NNZ))
J = J+1 J = J+1
ENDDO ENDDO
NZL = J - I NZL = J - I
CALL MRGSRT(NZL,JA(I),AUX,IRET) CALL MSORT_UP(NZL,JA(I),AUX,IRET)
IF (IRET.EQ.0) CALL REORDVN(NZL,AR(I),IA(I),JA(I),AUX) IF (IRET.EQ.0) CALL REORDVN(NZL,AR(I),IA(I),JA(I),AUX)
I = J I = J
ENDDO ENDDO
@ -343,7 +343,7 @@ c$$$ endif
ELSE IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'U') THEN ELSE IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'U') THEN
C .... Order with key IA ... C .... Order with key IA ...
CALL MRGSRT(NNZ,IA,AUX,IRET) CALL MSORT_UP(NNZ,IA,AUX,IRET)
IF (IRET.EQ.0) CALL REORDVN(NNZ,AR,IA,JA,AUX) IF (IRET.EQ.0) CALL REORDVN(NNZ,AR,IA,JA,AUX)
C .... Order with key IA2N ... C .... Order with key IA2N ...
I = 1 I = 1
@ -359,7 +359,7 @@ c$$$ + (J.LE.NNZ))
J = J+1 J = J+1
ENDDO ENDDO
NZL = J - I NZL = J - I
CALL MRGSRT(NZL,JA(I),AUX,IRET) CALL MSORT_UP(NZL,JA(I),AUX,IRET)
IF (IRET.EQ.0) CALL REORDVN(NZL,AR(I),IA(I),JA(I),AUX) IF (IRET.EQ.0) CALL REORDVN(NZL,AR(I),IA(I),JA(I),AUX)
I = J I = J
ENDDO ENDDO
@ -443,7 +443,7 @@ C ... Sum the duplicated element ...
ELSE IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'L') THEN ELSE IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'L') THEN
C .... Order with key IA ... C .... Order with key IA ...
CALL MRGSRT(NNZ,IA,AUX,IRET) CALL MSORT_UP(NNZ,IA,AUX,IRET)
IF (IRET.EQ.0) CALL REORDVN(NNZ,AR,IA,JA,AUX) IF (IRET.EQ.0) CALL REORDVN(NNZ,AR,IA,JA,AUX)
C .... Order with key IA2N ... C .... Order with key IA2N ...
I = 1 I = 1
@ -459,7 +459,7 @@ c$$$ + (J.LE.NNZ))
J = J+1 J = J+1
ENDDO ENDDO
NZL = J - I NZL = J - I
CALL MRGSRT(NZL,JA(I),AUX,IRET) CALL MSORT_UP(NZL,JA(I),AUX,IRET)
IF (IRET.EQ.0) CALL REORDVN(NZL,AR(I),IA(I),JA(I),AUX) IF (IRET.EQ.0) CALL REORDVN(NZL,AR(I),IA(I),JA(I),AUX)
I = J I = J
ENDDO ENDDO

@ -53,7 +53,7 @@ C Compute number of nnzero elements per row
C Sorting Array work C Sorting Array work
C ........................ C ........................
CALL MRGSRT(M,WORK,WORK(M+1),IRET) CALL MSORT_UP(M,WORK,WORK(M+1),IRET)
IF (IRET.EQ.0) THEN IF (IRET.EQ.0) THEN
C Construct IPERM Vector C Construct IPERM Vector
LP = WORK(M+1) LP = WORK(M+1)

@ -132,7 +132,7 @@ C SCALE = (UNITD.EQ.'L') ! meaningless
C .... Order with key IA1N.... C .... Order with key IA1N....
CALL MRGSRT(NNZ,IA1N,AUX,IRET) CALL MSORT_UP(NNZ,IA1N,AUX,IRET)
IF (IRET.EQ.0) CALL REORDVN(NNZ,ARN,IA1N,IA2N,AUX) IF (IRET.EQ.0) CALL REORDVN(NNZ,ARN,IA1N,IA2N,AUX)
C .... Order with key IA2N ... C .... Order with key IA2N ...
@ -144,7 +144,7 @@ C .... Order with key IA2N ...
J = J+1 J = J+1
ENDDO ENDDO
NZL = J - I NZL = J - I
CALL MRGSRT(NZL,IA2N(I),AUX,IRET) CALL MSORT_UP(NZL,IA2N(I),AUX,IRET)
IF (IRET.EQ.0) CALL REORDVN(NZL,ARN(I),IA1N(I),IA2N(I), IF (IRET.EQ.0) CALL REORDVN(NZL,ARN(I),IA1N(I),IA2N(I),
+ AUX) + AUX)
I = J I = J

@ -160,7 +160,7 @@ C FOR CURRENT GROUP
C C
L1 = IWORK2(LEV+1) - IWORK2(LEV) L1 = IWORK2(LEV+1) - IWORK2(LEV)
L0 = IWORK2(LEV) - 1 L0 = IWORK2(LEV) - 1
CALL MRGSRT(L1,IWORK1(IWORK2(LEV)),IPAT,IRET) CALL MSORT_UP(L1,IWORK1(IWORK2(LEV)),IPAT,IRET)
IF (IRET.EQ.0) THEN IF (IRET.EQ.0) THEN
NP = IPAT(1) NP = IPAT(1)
DO 280 L = 1, L1 DO 280 L = 1, L1

@ -146,7 +146,7 @@ c
if (debug) write(0,*) 'build check :',ia2n(ip1+psb_nnzt_) if (debug) write(0,*) 'build check :',ia2n(ip1+psb_nnzt_)
c .... order with key ia1n ... c .... order with key ia1n ...
call mrgsrt(nnz,ia1n,aux,iret) call msort_up(nnz,ia1n,aux,iret)
if (iret.eq.0) if (iret.eq.0)
+ call zreordvn3(nnz,arn,ia1n,ia2n,aux(ipx),aux) + call zreordvn3(nnz,arn,ia1n,ia2n,aux(ipx),aux)
c .... order with key ia2n ... c .... order with key ia2n ...
@ -159,7 +159,7 @@ c .... order with key ia2n ...
j = j+1 j = j+1
enddo enddo
nzl = j - i nzl = j - i
call mrgsrt(nzl,ia2n(i),aux,iret) call msort_up(nzl,ia2n(i),aux,iret)
if (iret.eq.0) call zreordvn3(nzl,arn(i),ia1n(i),ia2n(i), if (iret.eq.0) call zreordvn3(nzl,arn(i),ia1n(i),ia2n(i),
+ aux(ipx+i-1),aux) + aux(ipx+i-1),aux)
i = j i = j
@ -198,7 +198,7 @@ c ... sum the duplicated element ...
else else
c .... order with key ia1n ... c .... order with key ia1n ...
call mrgsrt(nnz,ia1n,aux,iret) call msort_up(nnz,ia1n,aux,iret)
if (iret.eq.0) call zreordvn(nnz,arn,ia1n,ia2n,aux) if (iret.eq.0) call zreordvn(nnz,arn,ia1n,ia2n,aux)
c .... order with key ia2n ... c .... order with key ia2n ...
@ -210,7 +210,7 @@ c .... order with key ia2n ...
j = j+1 j = j+1
enddo enddo
nzl = j - i nzl = j - i
call mrgsrt(nzl,ia2n(i),aux,iret) call msort_up(nzl,ia2n(i),aux,iret)
if (iret.eq.0) call zreordvn(nzl,arn(i),ia1n(i),ia2n(i), if (iret.eq.0) call zreordvn(nzl,arn(i),ia1n(i),ia2n(i),
+ aux) + aux)
i = j i = j

@ -157,7 +157,7 @@ c$$$ + ip1,ip2,nnz,ian2(ip1+nnzt_)
if (debug) write(0,*) 'Build check :',ian2(ip1+psb_nnzt_) if (debug) write(0,*) 'Build check :',ian2(ip1+psb_nnzt_)
C .... Order with key IA ... C .... Order with key IA ...
CALL MRGSRT(NNZ,IA,AUX,IRET) CALL MSORT_UP(NNZ,IA,AUX,IRET)
IF (IRET.EQ.0) CALL ZREORDVN3(NNZ,AR,IA,JA,AUX(IPX),AUX) IF (IRET.EQ.0) CALL ZREORDVN3(NNZ,AR,IA,JA,AUX(IPX),AUX)
if (debug) then if (debug) then
do i=1, nnz-1 do i=1, nnz-1
@ -183,7 +183,7 @@ c$$$ + (J.LE.NNZ))
J = J+1 J = J+1
ENDDO ENDDO
NZL = J - I NZL = J - I
CALL MRGSRT(NZL,JA(I),AUX,IRET) CALL MSORT_UP(NZL,JA(I),AUX,IRET)
IF (IRET.EQ.0) CALL ZREORDVN3(NZL,AR(I),IA(I),JA(I), IF (IRET.EQ.0) CALL ZREORDVN3(NZL,AR(I),IA(I),JA(I),
+ AUX(IPX+I-1),AUX) + AUX(IPX+I-1),AUX)
I = J I = J
@ -244,7 +244,7 @@ C ... Sum the duplicated element ...
ELSE ELSE
c$$$ write(0,*) 'Going for ELSE !!!?' c$$$ write(0,*) 'Going for ELSE !!!?'
C .... Order with key IA ... C .... Order with key IA ...
CALL MRGSRT(NNZ,IA,AUX,IRET) CALL MSORT_UP(NNZ,IA,AUX,IRET)
IF (IRET.EQ.0) CALL ZREORDVN(NNZ,AR,IA,JA,AUX) IF (IRET.EQ.0) CALL ZREORDVN(NNZ,AR,IA,JA,AUX)
C .... Order with key IA2N ... C .... Order with key IA2N ...
I = 1 I = 1
@ -260,7 +260,7 @@ c$$$ + (J.LE.NNZ))
J = J+1 J = J+1
ENDDO ENDDO
NZL = J - I NZL = J - I
CALL MRGSRT(NZL,JA(I),AUX,IRET) CALL MSORT_UP(NZL,JA(I),AUX,IRET)
IF (IRET.EQ.0) CALL ZREORDVN(NZL,AR(I),IA(I),JA(I),AUX) IF (IRET.EQ.0) CALL ZREORDVN(NZL,AR(I),IA(I),JA(I),AUX)
I = J I = J
ENDDO ENDDO
@ -337,7 +337,7 @@ C ... Sum the duplicated element ...
ELSE IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'U') THEN ELSE IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'U') THEN
C .... Order with key IA ... C .... Order with key IA ...
CALL MRGSRT(NNZ,IA,AUX,IRET) CALL MSORT_UP(NNZ,IA,AUX,IRET)
IF (IRET.EQ.0) CALL ZREORDVN(NNZ,AR,IA,JA,AUX) IF (IRET.EQ.0) CALL ZREORDVN(NNZ,AR,IA,JA,AUX)
C .... Order with key IA2N ... C .... Order with key IA2N ...
I = 1 I = 1
@ -353,7 +353,7 @@ c$$$ + (J.LE.NNZ))
J = J+1 J = J+1
ENDDO ENDDO
NZL = J - I NZL = J - I
CALL MRGSRT(NZL,JA(I),AUX,IRET) CALL MSORT_UP(NZL,JA(I),AUX,IRET)
IF (IRET.EQ.0) CALL ZREORDVN(NZL,AR(I),IA(I),JA(I),AUX) IF (IRET.EQ.0) CALL ZREORDVN(NZL,AR(I),IA(I),JA(I),AUX)
I = J I = J
ENDDO ENDDO
@ -437,7 +437,7 @@ C ... Sum the duplicated element ...
ELSE IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'L') THEN ELSE IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'L') THEN
C .... Order with key IA ... C .... Order with key IA ...
CALL MRGSRT(NNZ,IA,AUX,IRET) CALL MSORT_UP(NNZ,IA,AUX,IRET)
IF (IRET.EQ.0) CALL ZREORDVN(NNZ,AR,IA,JA,AUX) IF (IRET.EQ.0) CALL ZREORDVN(NNZ,AR,IA,JA,AUX)
C .... Order with key IA2N ... C .... Order with key IA2N ...
I = 1 I = 1
@ -453,7 +453,7 @@ c$$$ + (J.LE.NNZ))
J = J+1 J = J+1
ENDDO ENDDO
NZL = J - I NZL = J - I
CALL MRGSRT(NZL,JA(I),AUX,IRET) CALL MSORT_UP(NZL,JA(I),AUX,IRET)
IF (IRET.EQ.0) CALL ZREORDVN(NZL,AR(I),IA(I),JA(I),AUX) IF (IRET.EQ.0) CALL ZREORDVN(NZL,AR(I),IA(I),JA(I),AUX)
I = J I = J
ENDDO ENDDO

@ -72,7 +72,7 @@ Subroutine psb_dfixcoo(A,INFO,idir)
case(0) ! Row major order case(0) ! Row major order
call mrgsrt(nza,a%ia1(1),iaux(1),iret) call msort_up(nza,a%ia1(1),iaux(1),iret)
if (iret.eq.0) call reordvn(nza,a%aspk(1),a%ia1(1),a%ia2(1),iaux(1)) if (iret.eq.0) call reordvn(nza,a%aspk(1),a%ia1(1),a%ia2(1),iaux(1))
i = 1 i = 1
j = i j = i
@ -82,7 +82,7 @@ Subroutine psb_dfixcoo(A,INFO,idir)
if (j > nza) exit if (j > nza) exit
enddo enddo
nzl = j - i nzl = j - i
call mrgsrt(nzl,a%ia2(i),iaux(1),iret) call msort_up(nzl,a%ia2(i),iaux(1),iret)
if (iret.eq.0) & if (iret.eq.0) &
& call reordvn(nzl,a%aspk(i),a%ia1(i),a%ia2(i),iaux(1)) & call reordvn(nzl,a%aspk(i),a%ia1(i),a%ia2(i),iaux(1))
i = j i = j
@ -113,7 +113,7 @@ Subroutine psb_dfixcoo(A,INFO,idir)
case(1) ! Col major order case(1) ! Col major order
call mrgsrt(nza,a%ia2(1),iaux(1),iret) call msort_up(nza,a%ia2(1),iaux(1),iret)
if (iret.eq.0) call reordvn(nza,a%aspk(1),a%ia1(1),a%ia2(1),iaux(1)) if (iret.eq.0) call reordvn(nza,a%aspk(1),a%ia1(1),a%ia2(1),iaux(1))
i = 1 i = 1
j = i j = i
@ -123,7 +123,7 @@ Subroutine psb_dfixcoo(A,INFO,idir)
if (j > nza) exit if (j > nza) exit
enddo enddo
nzl = j - i nzl = j - i
call mrgsrt(nzl,a%ia1(i),iaux(1),iret) call msort_up(nzl,a%ia1(i),iaux(1),iret)
if (iret.eq.0) & if (iret.eq.0) &
& call reordvn(nzl,a%aspk(i),a%ia1(i),a%ia2(i),iaux(1)) & call reordvn(nzl,a%aspk(i),a%ia1(i),a%ia2(i),iaux(1))
i = j i = j

@ -72,7 +72,7 @@ Subroutine psb_zfixcoo(A,INFO,idir)
case(0) ! Row major order case(0) ! Row major order
call mrgsrt(nza,a%ia1(1),iaux(1),iret) call msort_up(nza,a%ia1(1),iaux(1),iret)
if (iret.eq.0) call zreordvn(nza,a%aspk(1),a%ia1(1),a%ia2(1),iaux(1)) if (iret.eq.0) call zreordvn(nza,a%aspk(1),a%ia1(1),a%ia2(1),iaux(1))
i = 1 i = 1
j = i j = i
@ -82,7 +82,7 @@ Subroutine psb_zfixcoo(A,INFO,idir)
if (j > nza) exit if (j > nza) exit
enddo enddo
nzl = j - i nzl = j - i
call mrgsrt(nzl,a%ia2(i),iaux(1),iret) call msort_up(nzl,a%ia2(i),iaux(1),iret)
if (iret.eq.0) & if (iret.eq.0) &
& call zreordvn(nzl,a%aspk(i),a%ia1(i),a%ia2(i),iaux(1)) & call zreordvn(nzl,a%aspk(i),a%ia1(i),a%ia2(i),iaux(1))
i = j i = j
@ -113,7 +113,7 @@ Subroutine psb_zfixcoo(A,INFO,idir)
case(1) ! Col major order case(1) ! Col major order
call mrgsrt(nza,a%ia2(1),iaux(1),iret) call msort_up(nza,a%ia2(1),iaux(1),iret)
if (iret.eq.0) call zreordvn(nza,a%aspk(1),a%ia1(1),a%ia2(1),iaux(1)) if (iret.eq.0) call zreordvn(nza,a%aspk(1),a%ia1(1),a%ia2(1),iaux(1))
i = 1 i = 1
j = i j = i
@ -123,7 +123,7 @@ Subroutine psb_zfixcoo(A,INFO,idir)
if (j > nza) exit if (j > nza) exit
enddo enddo
nzl = j - i nzl = j - i
call mrgsrt(nzl,a%ia1(i),iaux(1),iret) call msort_up(nzl,a%ia1(i),iaux(1),iret)
if (iret.eq.0) & if (iret.eq.0) &
& call zreordvn(nzl,a%aspk(i),a%ia1(i),a%ia2(i),iaux(1)) & call zreordvn(nzl,a%aspk(i),a%ia1(i),a%ia2(i),iaux(1))
i = j i = j

Loading…
Cancel
Save