page@swan.ulowell.edu (Bob Page) (11/03/88)
Submitted-by: strovink%galaxy-43@afit-ab.arpa (Mark A. Strovink)
Posting-number: Volume 2, Issue 47
Archive-name: applications/matlab/src.7
# This is a shell archive.
# Remove everything above and including the cut line.
# Then run the rest of the file through sh.
#----cut here-----cut here-----cut here-----cut here----#
#!/bin/sh
# shar: Shell Archiver
# Run the following text with /bin/sh to create:
# src-7
# This archive created: Wed Nov 2 16:20:51 1988
cat << \SHAR_EOF > src-7
CR = FLOP((ARS*BRS + AIS*BIS)/D)
CI = (AIS*BRS - ARS*BIS)/D
IF (CI .NE. 0.0D0) CI = FLOP(CI)
RETURN
END
SUBROUTINE WSIGN(XR,XI,YR,YI,ZR,ZI)
DOUBLE PRECISION XR,XI,YR,YI,ZR,ZI,PYTHAG,T
C IF Y .NE. 0, Z = X*Y/ABS(Y)
C IF Y .EQ. 0, Z = X
T = PYTHAG(YR,YI)
ZR = XR
ZI = XI
IF (T .NE. 0.0D0) CALL WMUL(YR/T,YI/T,ZR,ZI,ZR,ZI)
RETURN
END
SUBROUTINE WSQRT(XR,XI,YR,YI)
DOUBLE PRECISION XR,XI,YR,YI,S,TR,TI,PYTHAG,FLOP
C Y = SQRT(X) WITH YR .GE. 0.0 AND SIGN(YI) .EQ. SIGN(XI)
C
TR = XR
TI = XI
S = DSQRT(0.5D0*(PYTHAG(TR,TI) + DABS(TR)))
IF (TR .GE. 0.0D0) YR = FLOP(S)
IF (TI .LT. 0.0D0) S = -S
IF (TR .LE. 0.0D0) YI = FLOP(S)
IF (TR .LT. 0.0D0) YR = FLOP(0.5D0*(TI/YI))
IF (TR .GT. 0.0D0) YI = FLOP(0.5D0*(TI/YR))
RETURN
END
SUBROUTINE WLOG(XR,XI,YR,YI)
DOUBLE PRECISION XR,XI,YR,YI,T,R,PYTHAG
C Y = LOG(X)
R = PYTHAG(XR,XI)
IF (R .EQ. 0.0D0) CALL ERROR(32)
IF (R .EQ. 0.0D0) RETURN
T = DATAN2(XI,XR)
IF (XI.EQ.0.0D0 .AND. XR.LT.0.0D0) T = DABS(T)
YR = DLOG(R)
YI = T
RETURN
END
SUBROUTINE WATAN(XR,XI,YR,YI)
C Y = ATAN(X) = (I/2)*LOG((I+X)/(I-X))
DOUBLE PRECISION XR,XI,YR,YI,TR,TI
IF (XI .NE. 0.0D0) GO TO 10
YR = DATAN2(XR,1.0D0)
YI = 0.0D0
RETURN
10 IF (XR.NE.0.0D0 .OR. DABS(XI).NE.1.0D0) GO TO 20
CALL ERROR(32)
RETURN
20 CALL WDIV(XR,1.0D0+XI,-XR,1.0D0-XI,TR,TI)
CALL WLOG(TR,TI,TR,TI)
YR = -TI/2.0D0
YI = TR/2.0D0
RETURN
END
DOUBLE PRECISION FUNCTION WNRM2(N,XR,XI,INCX)
DOUBLE PRECISION XR(1),XI(1),PYTHAG,S
C NORM2(X)
S = 0.0D0
IF (N .LE. 0) GO TO 20
IX = 1
DO 10 I = 1, N
S = PYTHAG(S,XR(IX))
S = PYTHAG(S,XI(IX))
IX = IX + INCX
10 CONTINUE
20 WNRM2 = S
RETURN
END
DOUBLE PRECISION FUNCTION WASUM(N,XR,XI,INCX)
DOUBLE PRECISION XR(1),XI(1),S,FLOP
C NORM1(X)
S = 0.0D0
IF (N .LE. 0) GO TO 20
IX = 1
DO 10 I = 1, N
S = FLOP(S + DABS(XR(IX)) + DABS(XI(IX)))
IX = IX + INCX
10 CONTINUE
20 WASUM = S
RETURN
END
INTEGER FUNCTION IWAMAX(N,XR,XI,INCX)
DOUBLE PRECISION XR(1),XI(1),S,P
C INDEX OF NORMINF(X)
K = 0
IF (N .LE. 0) GO TO 20
K = 1
S = 0.0D0
IX = 1
DO 10 I = 1, N
P = DABS(XR(IX)) + DABS(XI(IX))
IF (P .GT. S) K = I
IF (P .GT. S) S = P
IX = IX + INCX
10 CONTINUE
20 IWAMAX = K
RETURN
END
SUBROUTINE WRSCAL(N,S,XR,XI,INCX)
DOUBLE PRECISION S,XR(1),XI(1),FLOP
IF (N .LE. 0) RETURN
IX = 1
DO 10 I = 1, N
XR(IX) = FLOP(S*XR(IX))
IF (XI(IX) .NE. 0.0D0) XI(IX) = FLOP(S*XI(IX))
IX = IX + INCX
10 CONTINUE
RETURN
END
SUBROUTINE WSCAL(N,SR,SI,XR,XI,INCX)
DOUBLE PRECISION SR,SI,XR(1),XI(1)
IF (N .LE. 0) RETURN
IX = 1
DO 10 I = 1, N
CALL WMUL(SR,SI,XR(IX),XI(IX),XR(IX),XI(IX))
IX = IX + INCX
10 CONTINUE
RETURN
END
SUBROUTINE WAXPY(N,SR,SI,XR,XI,INCX,YR,YI,INCY)
DOUBLE PRECISION SR,SI,XR(1),XI(1),YR(1),YI(1),FLOP
IF (N .LE. 0) RETURN
IF (SR .EQ. 0.0D0 .AND. SI .EQ. 0.0D0) RETURN
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
YR(IY) = FLOP(YR(IY) + SR*XR(IX) - SI*XI(IX))
YI(IY) = YI(IY) + SR*XI(IX) + SI*XR(IX)
IF (YI(IY) .NE. 0.0D0) YI(IY) = FLOP(YI(IY))
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
END
DOUBLE PRECISION FUNCTION WDOTUR(N,XR,XI,INCX,YR,YI,INCY)
DOUBLE PRECISION XR(1),XI(1),YR(1),YI(1),S,FLOP
S = 0.0D0
IF (N .LE. 0) GO TO 20
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
S = FLOP(S + XR(IX)*YR(IY) - XI(IX)*YI(IY))
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
20 WDOTUR = S
RETURN
END
DOUBLE PRECISION FUNCTION WDOTUI(N,XR,XI,INCX,YR,YI,INCY)
DOUBLE PRECISION XR(1),XI(1),YR(1),YI(1),S,FLOP
S = 0.0D0
IF (N .LE. 0) GO TO 20
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
S = S + XR(IX)*YI(IY) + XI(IX)*YR(IY)
IF (S .NE. 0.0D0) S = FLOP(S)
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
20 WDOTUI = S
RETURN
END
DOUBLE PRECISION FUNCTION WDOTCR(N,XR,XI,INCX,YR,YI,INCY)
DOUBLE PRECISION XR(1),XI(1),YR(1),YI(1),S,FLOP
S = 0.0D0
IF (N .LE. 0) GO TO 20
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
S = FLOP(S + XR(IX)*YR(IY) + XI(IX)*YI(IY))
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
20 WDOTCR = S
RETURN
END
DOUBLE PRECISION FUNCTION WDOTCI(N,XR,XI,INCX,YR,YI,INCY)
DOUBLE PRECISION XR(1),XI(1),YR(1),YI(1),S,FLOP
S = 0.0D0
IF (N .LE. 0) GO TO 20
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
S = S + XR(IX)*YI(IY) - XI(IX)*YR(IY)
IF (S .NE. 0.0D0) S = FLOP(S)
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
20 WDOTCI = S
RETURN
END
SUBROUTINE WCOPY(N,XR,XI,INCX,YR,YI,INCY)
DOUBLE PRECISION XR(1),XI(1),YR(1),YI(1)
IF (N .LE. 0) RETURN
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
YR(IY) = XR(IX)
YI(IY) = XI(IX)
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
END
SUBROUTINE WSET(N,XR,XI,YR,YI,INCY)
INTEGER N,INCY
DOUBLE PRECISION XR,XI,YR(1),YI(1)
IY = 1
IF (N .LE. 0 ) RETURN
DO 10 I = 1,N
YR(IY) = XR
YI(IY) = XI
IY = IY + INCY
10 CONTINUE
RETURN
END
SUBROUTINE WSWAP(N,XR,XI,INCX,YR,YI,INCY)
DOUBLE PRECISION XR(1),XI(1),YR(1),YI(1),T
IF (N .LE. 0) RETURN
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
T = XR(IX)
XR(IX) = YR(IY)
YR(IY) = T
T = XI(IX)
XI(IX) = YI(IY)
YI(IY) = T
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
END
SUBROUTINE RSET(N,DX,DY,INCY)
C
C COPIES A SCALAR, X, TO A SCALAR, Y.
DOUBLE PRECISION DX,DY(1)
C
IF (N.LE.0) RETURN
IY = 1
IF (INCY.LT.0) IY = (-N+1)*INCY + 1
DO 10 I = 1,N
DY(IY) = DX
IY = IY + INCY
10 CONTINUE
RETURN
END
SUBROUTINE RSWAP(N,X,INCX,Y,INCY)
DOUBLE PRECISION X(1),Y(1),T
IF (N .LE. 0) RETURN
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
T = X(IX)
X(IX) = Y(IY)
Y(IY) = T
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
END
SUBROUTINE RROT(N,DX,INCX,DY,INCY,C,S)
C
C APPLIES A PLANE ROTATION.
DOUBLE PRECISION DX(1),DY(1),DTEMP,C,S,FLOP
INTEGER I,INCX,INCY,IX,IY,N
C
IF (N.LE.0) RETURN
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 = FLOP(C*DX(IX) + S*DY(IY))
DY(IY) = FLOP(C*DY(IY) - S*DX(IX))
DX(IX) = DTEMP
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
END
SUBROUTINE RROTG(DA,DB,C,S)
C
C CONSTRUCT GIVENS PLANE ROTATION.
C
DOUBLE PRECISION DA,DB,C,S,RHO,PYTHAG,FLOP,R,Z
C
RHO = DB
IF ( DABS(DA) .GT. DABS(DB) ) RHO = DA
C = 1.0D0
S = 0.0D0
Z = 1.0D0
R = FLOP(DSIGN(PYTHAG(DA,DB),RHO))
IF (R .NE. 0.0D0) C = FLOP(DA/R)
IF (R .NE. 0.0D0) S = FLOP(DB/R)
IF ( DABS(DA) .GT. DABS(DB) ) Z = S
IF ( DABS(DB) .GE. DABS(DA) .AND. C .NE. 0.0D0 ) Z = FLOP(1.0D0/C)
DA = R
DB = Z
RETURN
END
LOGICAL FUNCTION EQID(X,Y)
C CHECK FOR EQUALITY OF TWO NAMES
INTEGER X(4),Y(4)
EQID = .TRUE.
DO 10 I = 1, 4
10 EQID = EQID .AND. (X(I).EQ.Y(I))
RETURN
END
SUBROUTINE PUTID(X,Y)
C STORE A NAME
INTEGER X(4),Y(4)
DO 10 I = 1, 4
10 X(I) = Y(I)
RETURN
END
DOUBLE PRECISION FUNCTION ROUND(X)
DOUBLE PRECISION X,Y,Z,E,H
DATA H/1.0D9/
Z = DABS(X)
Y = Z + 1.0D0
IF (Y .EQ. Z) GO TO 40
Y = 0.0D0
E = H
10 IF (E .GE. Z) GO TO 20
E = 2.0D0*E
GO TO 10
20 IF (E .LE. H) GO TO 30
IF (E .LE. Z) Y = Y + E
IF (E .LE. Z) Z = Z - E
E = E/2.0D0
GO TO 20
30 Z = IDINT(Z + 0.5D0)
Y = Y + Z
IF (X .LT. 0.0D0) Y = -Y
ROUND = Y
RETURN
40 ROUND = X
RETURN
END
FUNCTION DFLOAT(I)
C
C THIS IS THE AMIGA FUNCTION WHICH CONVERTS INTEGERS TO DOUBLE FLOATS
C
IMPLICIT NONE
DFLOAT = DBLE(I)
RETURN
END
SHAR_EOF
# End of shell archive
exit 0
--
Bob Page, U of Lowell CS Dept. page@swan.ulowell.edu ulowell!page
Have five nice days.