SUBROUTINE DAXPY( N, DA, DX, INCX, DY, INCY ) * * constant times a vector plus a vector. * uses unrolled loops for increments equal to one. * jack dongarra, linpack, 3/11/78. * * .. Scalar Arguments .. INTEGER INCX, INCY, N DOUBLE PRECISION DA * .. * .. Array Arguments .. DOUBLE PRECISION DX( 1 ), DY( 1 ) * .. * .. Local Scalars .. INTEGER I, IX, IY, M, MP1 * .. * .. Intrinsic Functions .. INTRINSIC MOD * .. * .. Executable Statements .. * IF( N.LE.0 ) \$ RETURN IF( DA.EQ.0.0D0 ) \$ RETURN IF( INCX.EQ.1 .AND. INCY.EQ.1 ) \$ GO TO 20 * * code for unequal increments or equal increments * not equal to 1 * 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 * * code for both increments equal to 1 * * * clean-up loop * 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 SUBROUTINE DCOPY( N, DX, INCX, DY, INCY ) * * copies a vector, x, to a vector, y. * uses unrolled loops for increments equal to one. * jack dongarra, linpack, 3/11/78. * * .. Scalar Arguments .. INTEGER INCX, INCY, N * .. * .. Array Arguments .. DOUBLE PRECISION DX( 1 ), DY( 1 ) * .. * .. Local Scalars .. INTEGER I, IX, IY, M, MP1 * .. * .. Intrinsic Functions .. INTRINSIC MOD * .. * .. Executable Statements .. * IF( N.LE.0 ) \$ RETURN IF( INCX.EQ.1 .AND. INCY.EQ.1 ) \$ GO TO 20 * * code for unequal increments or equal increments * not equal to 1 * 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 * * code for both increments equal to 1 * * * clean-up loop * 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 DOUBLE PRECISION FUNCTION DDOT( N, DX, INCX, DY, INCY ) * * forms the dot product of two vectors. * uses unrolled loops for increments equal to one. * jack dongarra, linpack, 3/11/78. * * .. Scalar Arguments .. INTEGER INCX, INCY, N * .. * .. Array Arguments .. DOUBLE PRECISION DX( 1 ), DY( 1 ) * .. * .. Local Scalars .. INTEGER I, IX, IY, M, MP1 DOUBLE PRECISION DTEMP * .. * .. Intrinsic Functions .. INTRINSIC MOD * .. * .. Executable Statements .. * DDOT = 0.0D0 DTEMP = 0.0D0 IF( N.LE.0 ) \$ RETURN IF( INCX.EQ.1 .AND. INCY.EQ.1 ) \$ GO TO 20 * * code for unequal increments or equal increments * not equal to 1 * 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 * * code for both increments equal to 1 * * * clean-up loop * 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 DOUBLE PRECISION FUNCTION DNRM2( N, DX, INCX ) * * euclidean norm of the n-vector stored in dx() with storage * increment incx . * if n .le. 0 return with result = 0. * if n .ge. 1 then incx must be .ge. 1 * * c.l.lawson, 1978 jan 08 * * four phase method using two built-in constants that are * hopefully applicable to all machines. * cutlo = maximum of dsqrt(u/eps) over all known machines. * cuthi = minimum of dsqrt(v) over all known machines. * where * eps = smallest no. such that eps + 1. .gt. 1. * u = smallest positive no. (underflow limit) * v = largest no. (overflow limit) * * brief outline of algorithm.. * * phase 1 scans zero components. * move to phase 2 when a component is nonzero and .le. cutlo * move to phase 3 when a component is .gt. cutlo * move to phase 4 when a component is .ge. cuthi/m * where m = n for x() real and m = 2*n for complex. * * values for cutlo and cuthi.. * from the environmental parameters listed in the imsl converter * document the limiting values are as follows.. * cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are * univac and dec at 2**(-103) * thus cutlo = 2**(-51) = 4.44089e-16 * cuthi, s.p. v = 2**127 for univac, honeywell, and dec. * thus cuthi = 2**(63.5) = 1.30438e19 * cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. * thus cutlo = 2**(-33.5) = 8.23181d-11 * cuthi, d.p. same as s.p. cuthi = 1.30438d19 * data cutlo, cuthi / 8.232d-11, 1.304d19 / * data cutlo, cuthi / 4.441e-16, 1.304e19 / * * .. Scalar Arguments .. INTEGER INCX, N * .. * .. Array Arguments .. DOUBLE PRECISION DX( 1 ) * .. * .. Local Scalars .. INTEGER I, IX, J, NEXT, NN DOUBLE PRECISION CUTHI, CUTLO, HITEST, ONE, SUM, \$ XMAX, ZERO * .. * .. Intrinsic Functions .. INTRINSIC DABS, DSQRT, FLOAT * .. * .. Data statements .. DATA ZERO, ONE / 0.0D0, 1.0D0 / DATA CUTLO, CUTHI / 8.232D-11, \$ 1.304D19 / * .. * .. Executable Statements .. * IF( N.GT.0 ) \$ GO TO 10 DNRM2 = ZERO GO TO 140 * 10 ASSIGN 30 TO NEXT SUM = ZERO * * begin main loop * IX = 1 IF( INCX.LT.0 ) \$ IX = 1 - ( N-1 )*INCX I = IX NN = IX + ( N-1 )*INCX 20 GO TO NEXT( 30, 40, 70, 80 ) 30 IF( DABS( DX( I ) ).GT.CUTLO ) \$ GO TO 110 ASSIGN 40 TO NEXT XMAX = ZERO * * phase 1. sum is zero * 40 IF( DX( I ).EQ.ZERO ) \$ GO TO 130 IF( DABS( DX( I ) ).GT.CUTLO ) \$ GO TO 110 * * prepare for phase 2. * ASSIGN 70 TO NEXT GO TO 60 * * prepare for phase 4. * 50 I = J ASSIGN 80 TO NEXT SUM = ( SUM / DX( I ) ) / DX( I ) 60 XMAX = DABS( DX( I ) ) GO TO 90 * * phase 2. sum is small. * scale to avoid destructive underflow. * 70 IF( DABS( DX( I ) ).GT.CUTLO ) \$ GO TO 100 * * common code for phases 2 and 4. * in phase 4 sum is large. scale to avoid overflow. * 80 IF( DABS( DX( I ) ).LE.XMAX ) \$ GO TO 90 SUM = ONE + SUM*( XMAX / DX( I ) )**2 XMAX = DABS( DX( I ) ) GO TO 130 * 90 SUM = SUM + ( DX( I ) / XMAX )**2 GO TO 130 * * prepare for phase 3. * 100 SUM = ( SUM*XMAX )*XMAX * * for real or d.p. set hitest = cuthi/n * for complex set hitest = cuthi/(2*n) * 110 HITEST = CUTHI / FLOAT( N ) * * phase 3. sum is mid-range. no scaling. * DO 120 J = I, NN, INCX IF( DABS( DX( J ) ).GE.HITEST ) \$ GO TO 50 SUM = SUM + DX( J )**2 120 CONTINUE DNRM2 = DSQRT( SUM ) GO TO 140 * 130 CONTINUE IF( I.NE.NN ) THEN I = I + INCX GO TO 20 ENDIF * * end of main loop. * * compute square root and adjust for scaling. * DNRM2 = XMAX*DSQRT( SUM ) 140 CONTINUE RETURN END SUBROUTINE DROT( N, DX, INCX, DY, INCY, C, S ) * * applies a plane rotation. * jack dongarra, linpack, 3/11/78. * * .. Scalar Arguments .. INTEGER INCX, INCY, N DOUBLE PRECISION C, S * .. * .. Array Arguments .. DOUBLE PRECISION DX( 1 ), DY( 1 ) * .. * .. Local Scalars .. INTEGER I, IX, IY DOUBLE PRECISION DTEMP * .. * .. Executable Statements .. * IF( N.LE.0 ) \$ RETURN IF( INCX.EQ.1 .AND. INCY.EQ.1 ) \$ GO TO 20 * * code for unequal increments or equal increments * not equal to 1 * 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 * * code for both increments equal to 1 * 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 SUBROUTINE DROTG( DA, DB, C, S ) * * construct givens plane rotation. * jack dongarra, linpack, 3/11/78. * modified 9/27/86. * * .. Scalar Arguments .. DOUBLE PRECISION C, DA, DB, S * .. * .. Local Scalars .. DOUBLE PRECISION R, ROE, SCALE, Z * .. * .. Intrinsic Functions .. INTRINSIC DABS, DSIGN, DSQRT * .. * .. Executable Statements .. * ROE = DB IF( DABS( DA ).GT.DABS( DB ) ) \$ ROE = DA SCALE = DABS( DA ) + DABS( DB ) IF( SCALE.NE.0.0D0 ) \$ GO TO 10 C = 1.0D0 S = 0.0D0 R = 0.0D0 GO TO 20 10 R = SCALE*DSQRT( ( DA / SCALE )**2+( DB / SCALE )**2 ) R = DSIGN( 1.0D0, ROE )*R C = DA / R S = DB / R 20 Z = S IF( DABS( C ).GT.0.0D0 .AND. DABS( C ).LE.S ) \$ Z = 1.0D0 / C DA = R DB = Z RETURN END SUBROUTINE DSCAL( N, DA, DX, INCX ) * * scales a vector by a constant. * uses unrolled loops for increment equal to one. * jack dongarra, linpack, 3/11/78. * * .. Scalar Arguments .. INTEGER INCX, N DOUBLE PRECISION DA * .. * .. Array Arguments .. DOUBLE PRECISION DX( 1 ) * .. * .. Local Scalars .. INTEGER I, IX, M, MP1, NINCX * .. * .. Intrinsic Functions .. INTRINSIC MOD * .. * .. Executable Statements .. * IF( N.LE.0 ) \$ RETURN IF( INCX.EQ.1 ) \$ GO TO 20 * * code for increment not equal to 1 * IX = 1 IF( INCX.LT.0 ) \$ IX = 1 - ( N-1 )*INCX NINCX = IX + ( N-1 )*INCX DO 10 I = IX, NINCX, INCX DX( I ) = DA*DX( I ) 10 CONTINUE RETURN * * code for increment equal to 1 * * * clean-up loop * 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 SUBROUTINE DSWAP( N, DX, INCX, DY, INCY ) * * interchanges two vectors. * uses unrolled loops for increments equal one. * jack dongarra, linpack, 3/11/78. * * .. Scalar Arguments .. INTEGER INCX, INCY, N * .. * .. Array Arguments .. DOUBLE PRECISION DX( 1 ), DY( 1 ) * .. * .. Local Scalars .. INTEGER I, IX, IY, M, MP1 DOUBLE PRECISION DTEMP * .. * .. Intrinsic Functions .. INTRINSIC MOD * .. * .. Executable Statements .. * IF( N.LE.0 ) \$ RETURN IF( INCX.EQ.1 .AND. INCY.EQ.1 ) \$ GO TO 20 * * code for unequal increments or equal increments * not equal to 1 * 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 * * code for both increments equal to 1 * * * clean-up loop * 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 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.0D0 DTEMP = 0.0D0 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 + DABS(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 + DABS(DX(I)) 30 CONTINUE IF( N .LT. 6 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,6 DTEMP = DTEMP + DABS(DX(I)) + DABS(DX(I + 1)) + DABS(DX(I + 2)) * + DABS(DX(I + 3)) + DABS(DX(I + 4)) + DABS(DX(I + 5)) 50 CONTINUE 60 DASUM = DTEMP RETURN END