+++ /dev/null
-C 10 NOV 94 - MWS - DNRM2: REMOVE FTNCHECK WARNINGS
-C 11 JUN 94 - MWS - INCLUDE A COPY OF DGEMV (LEVEL TWO ROUTINE)
-C 11 AUG 87 - MWS - SANITIZE FLOATING POINT CONSTANTS IN DNRM2
-C 26 MAR 87 - MWS - USE GENERIC SIGN IN DROTG
-C 28 NOV 86 - STE - SUPPLY ALL LEVEL ONE BLAS
-C 7 JUL 86 - JAB - SANITIZE FLOATING POINT CONSTANTS
-C
-C BASIC LINEAR ALGEBRA SUBPROGRAMS (BLAS) FROM LINPACK (LEVEL 1)
-C
-C THIS MODULE SHOULD BE COMPILED ONLY IF SPECIALLY CODED
-C VERSIONS OF THESE ROUTINES ARE NOT AVAILABLE ON THE TARGET MACHINE
-C
-C*MODULE BLAS1 *DECK DASUM
- DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX)
-C
-C TAKES THE SUM OF THE ABSOLUTE VALUES.
-C JACK DONGARRA, LINPACK, 3/11/78.
-C
- DOUBLE PRECISION DX(1),DTEMP
- INTEGER I,INCX,M,MP1,N,NINCX
-C
- DASUM = 0.0D+00
- DTEMP = 0.0D+00
- IF(N.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
- DTEMP = DTEMP + ABS(DX(I))
- 10 CONTINUE
- DASUM = DTEMP
- RETURN
-C
-C CODE FOR INCREMENT EQUAL TO 1
-C
-C
-C CLEAN-UP LOOP
-C
- 20 M = MOD(N,6)
- IF( M .EQ. 0 ) GO TO 40
- DO 30 I = 1,M
- DTEMP = DTEMP + ABS(DX(I))
- 30 CONTINUE
- IF( N .LT. 6 ) GO TO 60
- 40 MP1 = M + 1
- DO 50 I = MP1,N,6
- DTEMP = DTEMP + ABS(DX(I)) + ABS(DX(I + 1)) + ABS(DX(I + 2))
- * + ABS(DX(I + 3)) + ABS(DX(I + 4)) + ABS(DX(I + 5))
- 50 CONTINUE
- 60 DASUM = DTEMP
- RETURN
- END
-C*MODULE BLAS1 *DECK DAXPY
- SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
- IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- DIMENSION DX(1),DY(1)
-C
-C CONSTANT TIMES A VECTOR PLUS A VECTOR.
-C DY(I) = DY(I) + DA * DX(I)
-C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
-C JACK DONGARRA, LINPACK, 3/11/78.
-C
- IF(N.LE.0)RETURN
- IF (DA .EQ. 0.0D+00) 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) = DY(IY) + DA*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,4)
- IF( M .EQ. 0 ) GO TO 40
- DO 30 I = 1,M
- DY(I) = DY(I) + DA*DX(I)
- 30 CONTINUE
- IF( N .LT. 4 ) RETURN
- 40 MP1 = M + 1
- DO 50 I = MP1,N,4
- DY(I) = DY(I) + DA*DX(I)
- DY(I + 1) = DY(I + 1) + DA*DX(I + 1)
- DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
- DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
- 50 CONTINUE
- RETURN
- END
-C*MODULE BLAS1 *DECK DCOPY
- SUBROUTINE DCOPY(N,DX,INCX,DY,INCY)
- IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- DIMENSION DX(*),DY(*)
-C
-C COPIES A VECTOR.
-C DY(I) <== DX(I)
-C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
-C JACK DONGARRA, LINPACK, 3/11/78.
-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
-C*MODULE BLAS1 *DECK DDOT
- DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
- IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- DIMENSION DX(1),DY(1)
-C
-C FORMS THE DOT PRODUCT OF TWO VECTORS.
-C DOT = DX(I) * DY(I)
-C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
-C JACK DONGARRA, LINPACK, 3/11/78.
-C
- DDOT = 0.0D+00
- DTEMP = 0.0D+00
- 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
-C*MODULE BLAS1 *DECK DNRM2
- DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX)
- INTEGER NEXT
- DOUBLE PRECISION DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE
- DATA ZERO, ONE /0.0D+00, 1.0D+00/
-C
-C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE
-C INCREMENT INCX .
-C IF N .LE. 0 RETURN WITH RESULT = 0.
-C IF N .GE. 1 THEN INCX MUST BE .GE. 1
-C
-C C.L.LAWSON, 1978 JAN 08
-C
-C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE
-C HOPEFULLY APPLICABLE TO ALL MACHINES.
-C CUTLO = MAXIMUM OF SQRT(U/EPS) OVER ALL KNOWN MACHINES.
-C CUTHI = MINIMUM OF SQRT(V) OVER ALL KNOWN MACHINES.
-C WHERE
-C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
-C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT)
-C V = LARGEST NO. (OVERFLOW LIMIT)
-C
-C BRIEF OUTLINE OF ALGORITHM..
-C
-C PHASE 1 SCANS ZERO COMPONENTS.
-C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
-C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
-C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
-C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
-C
-C VALUES FOR CUTLO AND CUTHI..
-C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
-C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
-C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE
-C UNIVAC AND DEC AT 2**(-103)
-C THUS CUTLO = 2**(-51) = 4.44089E-16
-C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
-C THUS CUTHI = 2**(63.5) = 1.30438E19
-C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
-C THUS CUTLO = 2**(-33.5) = 8.23181D-11
-C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D+19
-C DATA CUTLO, CUTHI / 8.232D-11, 1.304D+19 /
-C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 /
- DATA CUTLO, CUTHI / 8.232D-11, 1.304D+19 /
-C
- J=0
- IF(N .GT. 0) GO TO 10
- DNRM2 = ZERO
- GO TO 300
-C
- 10 ASSIGN 30 TO NEXT
- SUM = ZERO
- NN = N * INCX
-C BEGIN MAIN LOOP
- I = 1
- 20 GO TO NEXT,(30, 50, 70, 110)
- 30 IF( ABS(DX(I)) .GT. CUTLO) GO TO 85
- ASSIGN 50 TO NEXT
- XMAX = ZERO
-C
-C PHASE 1. SUM IS ZERO
-C
- 50 IF( DX(I) .EQ. ZERO) GO TO 200
- IF( ABS(DX(I)) .GT. CUTLO) GO TO 85
-C
-C PREPARE FOR PHASE 2.
- ASSIGN 70 TO NEXT
- GO TO 105
-C
-C PREPARE FOR PHASE 4.
-C
- 100 I = J
- ASSIGN 110 TO NEXT
- SUM = (SUM / DX(I)) / DX(I)
- 105 XMAX = ABS(DX(I))
- GO TO 115
-C
-C PHASE 2. SUM IS SMALL.
-C SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
-C
- 70 IF( ABS(DX(I)) .GT. CUTLO ) GO TO 75
-C
-C COMMON CODE FOR PHASES 2 AND 4.
-C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW.
-C
- 110 IF( ABS(DX(I)) .LE. XMAX ) GO TO 115
- SUM = ONE + SUM * (XMAX / DX(I))**2
- XMAX = ABS(DX(I))
- GO TO 200
-C
- 115 SUM = SUM + (DX(I)/XMAX)**2
- GO TO 200
-C
-C
-C PREPARE FOR PHASE 3.
-C
- 75 SUM = (SUM * XMAX) * XMAX
-C
-C
-C FOR REAL OR D.P. SET HITEST = CUTHI/N
-C FOR COMPLEX SET HITEST = CUTHI/(2*N)
-C
- 85 HITEST = CUTHI/N
-C
-C PHASE 3. SUM IS MID-RANGE. NO SCALING.
-C
- DO 95 J =I,NN,INCX
- IF(ABS(DX(J)) .GE. HITEST) GO TO 100
- 95 SUM = SUM + DX(J)**2
- DNRM2 = SQRT( SUM )
- GO TO 300
-C
- 200 CONTINUE
- I = I + INCX
- IF ( I .LE. NN ) GO TO 20
-C
-C END OF MAIN LOOP.
-C
-C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
-C
- DNRM2 = XMAX * SQRT(SUM)
- 300 CONTINUE
- RETURN
- END
-C*MODULE BLAS1 *DECK DROT
- SUBROUTINE DROT (N,DX,INCX,DY,INCY,C,S)
- IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- DIMENSION DX(1),DY(1)
-C
-C APPLIES A PLANE ROTATION.
-C DX(I) = C*DX(I) + S*DY(I)
-C DY(I) = -S*DX(I) + C*DY(I)
-C JACK DONGARRA, LINPACK, 3/11/78.
-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 NOT EQUAL
-C 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 = C*DX(IX) + S*DY(IY)
- DY(IY) = C*DY(IY) - S*DX(IX)
- DX(IX) = DTEMP
- IX = IX + INCX
- IY = IY + INCY
- 10 CONTINUE
- RETURN
-C
-C CODE FOR BOTH INCREMENTS EQUAL TO 1
-C
- 20 DO 30 I = 1,N
- DTEMP = C*DX(I) + S*DY(I)
- DY(I) = C*DY(I) - S*DX(I)
- DX(I) = DTEMP
- 30 CONTINUE
- RETURN
- END
-C*MODULE BLAS1 *DECK DROTG
- SUBROUTINE DROTG(DA,DB,C,S)
-C
-C CONSTRUCT GIVENS PLANE ROTATION.
-C JACK DONGARRA, LINPACK, 3/11/78.
-C
- DOUBLE PRECISION DA,DB,C,S,ROE,SCALE,R,Z
- DOUBLE PRECISION ZERO, ONE
-C
- PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
-C
-C-----------------------------------------------------------------------
-C
-C
- ROE = DB
- IF( ABS(DA) .GT. ABS(DB) ) ROE = DA
- SCALE = ABS(DA) + ABS(DB)
- IF( SCALE .NE. ZERO ) GO TO 10
- C = ONE
- S = ZERO
- R = ZERO
- GO TO 20
-C
- 10 R = SCALE*SQRT((DA/SCALE)**2 + (DB/SCALE)**2)
- R = SIGN(ONE,ROE)*R
- C = DA/R
- S = DB/R
- 20 Z = ONE
- IF( ABS(DA) .GT. ABS(DB) ) Z = S
- IF( ABS(DB) .GE. ABS(DA) .AND. C .NE. ZERO ) Z = ONE/C
- DA = R
- DB = Z
- RETURN
- END
-C*MODULE BLAS1 *DECK DSCAL
- SUBROUTINE DSCAL(N,DA,DX,INCX)
- IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- DIMENSION DX(1)
-C
-C SCALES A VECTOR BY A CONSTANT.
-C DX(I) = DA * DX(I)
-C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
-C JACK DONGARRA, LINPACK, 3/11/78.
-C
- IF(N.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
-C*MODULE BLAS1 *DECK DSWAP
- SUBROUTINE DSWAP (N,DX,INCX,DY,INCY)
- IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- DIMENSION DX(1),DY(1)
-C
-C INTERCHANGES TWO VECTORS.
-C DX(I) <==> DY(I)
-C USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE.
-C JACK DONGARRA, LINPACK, 3/11/78.
-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 NOT EQUAL
-C 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 = DX(IX)
- DX(IX) = DY(IY)
- DY(IY) = DTEMP
- 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,3)
- IF( M .EQ. 0 ) GO TO 40
- DO 30 I = 1,M
- DTEMP = DX(I)
- DX(I) = DY(I)
- DY(I) = DTEMP
- 30 CONTINUE
- IF( N .LT. 3 ) RETURN
- 40 MP1 = M + 1
- DO 50 I = MP1,N,3
- DTEMP = DX(I)
- DX(I) = DY(I)
- DY(I) = DTEMP
- DTEMP = DX(I + 1)
- DX(I + 1) = DY(I + 1)
- DY(I + 1) = DTEMP
- DTEMP = DX(I + 2)
- DX(I + 2) = DY(I + 2)
- DY(I + 2) = DTEMP
- 50 CONTINUE
- RETURN
- END
-C*MODULE BLAS1 *DECK IDAMAX
- INTEGER FUNCTION IDAMAX(N,DX,INCX)
- IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- DIMENSION DX(1)
-C
-C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE.
-C JACK DONGARRA, LINPACK, 3/11/78.
-C
- IDAMAX = 0
- IF( N .LT. 1 ) 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
- RMAX = ABS(DX(1))
- IX = IX + INCX
- DO 10 I = 2,N
- IF(ABS(DX(IX)).LE.RMAX) GO TO 5
- IDAMAX = I
- RMAX = ABS(DX(IX))
- 5 IX = IX + INCX
- 10 CONTINUE
- RETURN
-C
-C CODE FOR INCREMENT EQUAL TO 1
-C
- 20 RMAX = ABS(DX(1))
- DO 30 I = 2,N
- IF(ABS(DX(I)).LE.RMAX) GO TO 30
- IDAMAX = I
- RMAX = ABS(DX(I))
- 30 CONTINUE
- RETURN
- END
-C*MODULE BLAS *DECK DGEMV
- SUBROUTINE DGEMV(FORMA,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
- IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- CHARACTER*1 FORMA
- DIMENSION A(LDA,*),X(*),Y(*)
- PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
-C
-C CLONE OF -DGEMV- WRITTEN BY MIKE SCHMIDT
-C
- LOCY = 1
- IF(FORMA.EQ.'T') GO TO 200
-C
-C Y = ALPHA * A * X + BETA * Y
-C
- IF(ALPHA.EQ.ONE .AND. BETA.EQ.ZERO) THEN
- DO 110 I=1,M
- Y(LOCY) = DDOT(N,A(I,1),LDA,X,INCX)
- LOCY = LOCY+INCY
- 110 CONTINUE
- ELSE
- DO 120 I=1,M
- Y(LOCY) = ALPHA*DDOT(N,A(I,1),LDA,X,INCX) + BETA*Y(LOCY)
- LOCY = LOCY+INCY
- 120 CONTINUE
- END IF
- RETURN
-C
-C Y = ALPHA * A-TRANSPOSE * X + BETA * Y
-C
- 200 CONTINUE
- IF(ALPHA.EQ.ONE .AND. BETA.EQ.ZERO) THEN
- DO 210 I=1,N
- Y(LOCY) = DDOT(M,A(1,I),1,X,INCX)
- LOCY = LOCY+INCY
- 210 CONTINUE
- ELSE
- DO 220 I=1,N
- Y(LOCY) = ALPHA*DDOT(M,A(1,I),1,X,INCX) + BETA*Y(LOCY)
- LOCY = LOCY+INCY
- 220 CONTINUE
- END IF
- RETURN
- END
projects / unres.git / blobdiff