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

Removed old merge sort from sparker, duplicate BLAS 1 code.

Browse Source
psblas3-type-indexed
Salvatore Filippone 19 years ago
parent 8ea7a2fddf
commit e6f0038390
8 changed files with 3 additions and 352 deletions
Show Stats Download Patch File Download Diff File

4
src/serial/aux/mrgsrt.f
Unescape Escape View File

@ -152,14 +152,14 @@ 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 100 P = 1,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
100 CONTINUE
END DO
L(T) = 0
L(N) = 0
C See if the input was already sorted

107
src/serial/aux/msrtrw.f
Unescape Escape View File

@ -1,107 +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
C Original version from Sparker
C
C msrtrw.f
C Author: Carlo Vittoli
C Date: May 19, 1994
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C Subroutine msrtrw C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C Purpose: Sort rows by column indices (merge sorting) C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C Parameters: C
C Input: C
C C
C Output: C
C C
C Others: C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C Algorithm: C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C References: C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE msrtrw(m,a,ia1,ia2,work,lwork,awork,lawork,ierrv)
IMPLICIT NONE
C .. Scalar Arguments ..
INTEGER m, lwork, lawork
C .. Array Arguments ..
DOUBLE PRECISION a(*), awork(m)
INTEGER ia1(*), ia2(*), work(2*m), ierrv(*)
C .. Local Scalars ..
INTEGER i1, i, j2, lrow, jend, jstart, irow, iret
C .. External Subroutines ..
EXTERNAL xsperr
C .. Executable Statements ..
if(lawork.lt.m) then
call xsperr('LWORK ',lawork,8,'MSRTRW',IERRV)
goto 9999
endif
if(lwork.lt.2*m) then
call xsperr('LWORK ',lwork,6,'MSRTRW',IERRV)
goto 9999
endif
C Start
DO 160 IROW = 1, M
JSTART = ia2(IROW)
JEND = ia2(IROW+1) - 1
LROW = JEND - JSTART + 1
call mrgsrt(lrow,ia1(jstart),work,iret)
if (iret.eq.0) then
I1 = WORK(1)
DO 20 I = 1, LROW
work(m+I) = I1 + JSTART - 1
I1 = WORK(I1+1)
20 CONTINUE
DO 40 I = 1, LROW
WORK(I) = ia1(work(m+I))
AWORK(I) = A(work(m+I))
40 CONTINUE
DO 60 I = 1, LROW
J2 = I + JSTART - 1
A(J2) = AWORK(I)
ia1(J2) = WORK(I)
60 CONTINUE
endif
160 continue
9999 continue
RETURN
END

3
src/serial/f77/Makefile
Unescape Escape View File

@ -5,8 +5,7 @@ include ../../../Make.inc
#
FOBJS = daxpby.o dcsmm.o dcsnmi.o dcsrp.o dcssm.o \
dcsupd.o dgelp.o dlpupd.o dswmm.o dswprt.o \
dswsm.o smmp.o ddot.o dscal.o dcsrws.o idamax.o \
dnrm2.o dcopy.o
dswsm.o smmp.o dcsrws.o
#zcsrck.o zcrnrmi.o zcsrmm.o zsrmv.o zcsrsm.o zsrsv.o

50
src/serial/f77/dcopy.f
Unescape Escape View File

@ -1,50 +0,0 @@
subroutine dcopy(n,dx,incx,dy,incy)
c
c copies a vector, x, to a vector, y.
c uses unrolled loops for increments equal to one.
c jack dongarra, linpack, 3/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double precision dx(*),dy(*)
integer i,incx,incy,ix,iy,m,mp1,n
c
if(n.le.0)return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
dy(iy) = dx(ix)
ix = ix + incx
iy = iy + incy
10 continue
return
c
c code for both increments equal to 1
c
c
c clean-up loop
c
20 m = mod(n,7)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dy(i) = dx(i)
30 continue
if( n .lt. 7 ) return
40 mp1 = m + 1
do 50 i = mp1,n,7
dy(i) = dx(i)
dy(i + 1) = dx(i + 1)
dy(i + 2) = dx(i + 2)
dy(i + 3) = dx(i + 3)
dy(i + 4) = dx(i + 4)
dy(i + 5) = dx(i + 5)
dy(i + 6) = dx(i + 6)
50 continue
return
end

49
src/serial/f77/ddot.f
Unescape Escape View File

@ -1,49 +0,0 @@
double precision function ddot(n,dx,incx,dy,incy)
c
c forms the dot product of two vectors.
c uses unrolled loops for increments equal to one.
c jack dongarra, linpack, 3/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double precision dx(*),dy(*),dtemp
integer i,incx,incy,ix,iy,m,mp1,n
c
ddot = 0.0d0
dtemp = 0.0d0
if(n.le.0)return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
dtemp = dtemp + dx(ix)*dy(iy)
ix = ix + incx
iy = iy + incy
10 continue
ddot = dtemp
return
c
c code for both increments equal to 1
c
c
c clean-up loop
c
20 m = mod(n,5)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dtemp = dtemp + dx(i)*dy(i)
30 continue
if( n .lt. 5 ) go to 60
40 mp1 = m + 1
do 50 i = mp1,n,5
dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
* dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
50 continue
60 ddot = dtemp
return
end

60
src/serial/f77/dnrm2.f
Unescape Escape View File

@ -1,60 +0,0 @@
DOUBLE PRECISION FUNCTION DNRM2 ( N, X, INCX )
* .. Scalar Arguments ..
INTEGER INCX, N
* .. Array Arguments ..
DOUBLE PRECISION X( * )
* ..
*
* DNRM2 returns the euclidean norm of a vector via the function
* name, so that
*
* DNRM2 := sqrt( x'*x )
*
*
*
* -- This version written on 25-October-1982.
* Modified on 14-October-1993 to inline the call to DLASSQ.
* Sven Hammarling, Nag Ltd.
*
*
* .. Parameters ..
DOUBLE PRECISION ONE , ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
* .. Local Scalars ..
INTEGER IX
DOUBLE PRECISION ABSXI, NORM, SCALE, SSQ
* .. Intrinsic Functions ..
INTRINSIC ABS, SQRT
* ..
* .. Executable Statements ..
IF( N.LT.1 .OR. INCX.LT.1 )THEN
NORM = ZERO
ELSE IF( N.EQ.1 )THEN
NORM = ABS( X( 1 ) )
ELSE
SCALE = ZERO
SSQ = ONE
* The following loop is equivalent to this call to the LAPACK
* auxiliary routine:
* CALL DLASSQ( N, X, INCX, SCALE, SSQ )
*
DO 10, IX = 1, 1 + ( N - 1 )*INCX, INCX
IF( X( IX ).NE.ZERO )THEN
ABSXI = ABS( X( IX ) )
IF( SCALE.LT.ABSXI )THEN
SSQ = ONE + SSQ*( SCALE/ABSXI )**2
SCALE = ABSXI
ELSE
SSQ = SSQ + ( ABSXI/SCALE )**2
END IF
END IF
10 CONTINUE
NORM = SCALE * SQRT( SSQ )
END IF
*
DNRM2 = NORM
RETURN
*
* End of DNRM2.
*
END

43
src/serial/f77/dscal.f
Unescape Escape View File

@ -1,43 +0,0 @@
subroutine dscal(n,da,dx,incx)
c
c scales a vector by a constant.
c uses unrolled loops for increment equal to one.
c jack dongarra, linpack, 3/11/78.
c modified 3/93 to return if incx .le. 0.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double precision da,dx(*)
integer i,incx,m,mp1,n,nincx
c
if( n.le.0 .or. incx.le.0 )return
if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c
nincx = n*incx
do 10 i = 1,nincx,incx
dx(i) = da*dx(i)
10 continue
return
c
c code for increment equal to 1
c
c
c clean-up loop
c
20 m = mod(n,5)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dx(i) = da*dx(i)
30 continue
if( n .lt. 5 ) return
40 mp1 = m + 1
do 50 i = mp1,n,5
dx(i) = da*dx(i)
dx(i + 1) = da*dx(i + 1)
dx(i + 2) = da*dx(i + 2)
dx(i + 3) = da*dx(i + 3)
dx(i + 4) = da*dx(i + 4)
50 continue
return
end

39
src/serial/f77/idamax.f
Unescape Escape View File

@ -1,39 +0,0 @@
integer function idamax(n,dx,incx)
c
c finds the index of element having max. absolute value.
c jack dongarra, linpack, 3/11/78.
c modified 3/93 to return if incx .le. 0.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double precision dx(*),dmax
integer i,incx,ix,n
c
idamax = 0
if( n.lt.1 .or. incx.le.0 ) return
idamax = 1
if(n.eq.1)return
if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c
ix = 1
dmax = dabs(dx(1))
ix = ix + incx
do 10 i = 2,n
if(dabs(dx(ix)).le.dmax) go to 5
idamax = i
dmax = dabs(dx(ix))
5 ix = ix + incx
10 continue
return
c
c code for increment equal to 1
c
20 dmax = dabs(dx(1))
do 30 i = 2,n
if(dabs(dx(i)).le.dmax) go to 30
idamax = i
dmax = dabs(dx(i))
30 continue
return
end
Write Preview
Loading…
Cancel
Save