--- /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