V10/cmd/matlab/sys.prime
C PROGRAM MAIN
C
C ******** THIS DRIVER WAS SET UP FOR THE PRIME 400 AT
C ******** THE CENTER FOR SCIENTIFIC COMPUTATION AND INTERACTIVE GRAPHICS
C ******** DREXEL UNIVERSITY
C ******** PHILADELPHIA, PA 19104
C
C THIS PROGRAM AND ASSOCIATED SUBROUTINES SHOULD BE COMPILED
C USING THE LONG INTEGER OPTION.
C
INTEGER*2 I
LOGICAL L,CLOS$A
CALL MATLAB(0)
I4=4
DO 100 I=1,12
L = CLOS$A(I)
100 CONTINUE
CALL EXIT
END
SUBROUTINE FILES(LUNIT,NAME,IOSTAT)
INTEGER LUNIT,NAME(32)
INTEGER*2 I1,I2,I3,I4,I5,I8,ITYPE,ICODE,IUNIT,I16
INTEGER*4 NAM8(2)
DATA I1,I3,I4,I5,I8,I16/1,3,4,5,8,16/
C
C SYSTEM DEPENDENT ROUTINE TO ALLOCATE FILES
C LUNIT = LOGICAL UNIT NUMBER
C NAME = FILE NAME, 1 CHARACTER PER WORD
C
C MODIFY SUBROUTINE MATLAB SO THAT RTE = 1 AND WTE = 1
C ATTACH UNIT 9 TO HELP FILE
C
L = -LUNIT
IF (LUNIT .LT. 0) REWIND L
IF (LUNIT .LT. 0) RETURN
C
ENCODE(8,10,NAM8) (NAME(I),I=1,8)
10 FORMAT(8A1)
IUNIT = LUNIT-4
CALL SRCH$$(I4,NAM8,0,IUNIT,ITYPE,ICODE)
CALL SRCH$$(I3,NAM8,I8,IUNIT,ITYPE,ICODE)
RETURN
END
SUBROUTINE SAVLOD(LUNIT,ID,M,N,IMG,JOB,XREAL,XIMAG)
INTEGER LUNIT,ID(4),M,N,IMG,JOB
DOUBLE PRECISION XREAL(1),XIMAG(1)
C
C IMPLEMENT SAVE AND LOAD
C LUNIT = LOGICAL UNIT NUMBER
C ID = NAME, FORMAT 4A1
C M, N = DIMENSIONS
C IMG = NONZERO IF XIMAG IS NONZERO
C JOB = 0 FOR SAVE
C = SPACE AVAILABLE FOR LOAD
C XREAL, XIMAG = REAL AND OPTIONAL IMAGINARY PARTS
C
C SYSTEM DEPENDENT FORMATS
101 FORMAT(4A1,3I4)
102 FORMAT(4D25.18)
C
IF (JOB .GT. 0) GO TO 20
C
C SAVE
10 WRITE(LUNIT,101) ID,M,N,IMG
DO 15 J = 1, N
K = (J-1)*M+1
L = J*M
WRITE(LUNIT,102) (XREAL(I),I=K,L)
IF (IMG .NE. 0) WRITE(LUNIT,102) (XIMAG(I),I=K,L)
15 CONTINUE
RETURN
C
C LOAD
20 READ(LUNIT,101,END=30) ID,M,N,IMG
IF (M*N .GT. JOB) GO TO 30
DO 25 J = 1, N
K = (J-1)*M+1
L = J*M
READ(LUNIT,102,END=30) (XREAL(I),I=K,L)
IF (IMG .NE. 0) READ(LUNIT,102,END=30) (XIMAG(I),I=K,L)
25 CONTINUE
RETURN
C
C END OF FILE
30 M = 0
N = 0
RETURN
END
SUBROUTINE FORMZ(LUNIT,X,Y)
DOUBLE PRECISION X,Y,XHEX(2),YHEX(2)
INTEGER*2 I4
DATA I4/4/
C
C SYSTEM DEPENDENT ROUTINE TO PRINT WITH Z FORMAT
C
CALL HEXCHR(X,I4,XHEX)
IF(Y.EQ.0.D0)GO TO 20
CALL HEXCHR(Y,I4,YHEX)
WRITE(LUNIT,10)XHEX,YHEX
RETURN
20 WRITE(LUNIT,10)XHEX
RETURN
10 FORMAT(2A8,4X,2A8)
END
SUBROUTINE HEXCHR(VALUE,VALWDS,STRING)
C* VALUE IS THE VALUE TO BE CONVERTED TO HEXADECIMAL. ANY NUMBER OF WORDS.
C* VALWDS IS THE NUMBER OF WORDS IN THE VALUE TO BE TRANSLATED:
C*
C* INTEGER*2 -- 1
C* INTEGER*4 -- 2
C* REAL*4 -- 2
C* REAL*8 -- 4
C*
C**===> ETC.
C* STRING IS AN ARRAY INTO WHICH THE CHARACTERS WILL BE PACKED.
C* (IT DOESN'T MATTER WHAT KIND OF ARRAY IT IS).
C*
C*
C* WARNING: IF THE STRING PROVIDED IS NOT LONG ENOUGH (IT HAS TO BE TWICE
C* AS LONG IN TERMS OF BITS OR WORDS AS THE VALUE ARRAY), THERE
C* IS NO WAY TO DETECT IF IT RUNS OFF AND DESTROYS OTHER MEMORY.
C* (THIS IS PRIME STANDARD, BY THE WAY.)
C* ALSO, ANYTHING BEYOND THE 2*VALWDS ELEMENT OF THE STRING ARRAY
C* IS MEANINGLESS.
C*
C*
INTEGER*2 VALUE(1),STRING(1),VALWDS,ONE,TWO,THREE,FOUR,I
DO 10 I=1,VALWDS
STRING(2*I-1) = 0
STRING(2*I) = 0
ONE = RS(LT(VALUE(I),4),12)
TWO = RS(LT(RT(VALUE(I),12),8),8)
THREE = RS(LT(RT(VALUE(I),8),12),4)
FOUR = RT(VALUE(I),4)
IF (ONE .LE. :011) STRING(2*I-1) =
+ OR(LS(OR(:260,ONE),8),STRING(2*I-1))
IF (ONE .GT. :011) STRING(2*I-1) =
+ OR(LS(ONE + :267,8),STRING(2*I-1))
IF (TWO .LE. :011) STRING(2*I-1) =
+ OR(TWO,:260,STRING(2*I-1))
IF (TWO .GT. :011) STRING(2*I-1) =
+ OR(TWO+ :267,STRING(2*I-1))
IF (THREE .LE. :011) STRING(2*I) =
+ OR(LS(OR(:260,THREE),8),STRING(2*I))
IF (THREE .GT. :011) STRING(2*I) =
+ OR(LS(THREE + :267,8),STRING(2*I))
IF (FOUR .LE. :011) STRING(2*I) =
+ OR(FOUR,:260,STRING(2*I))
IF (FOUR .GT. :011) STRING(2*I) =
+ OR(FOUR+ :267,STRING(2*I))
10 CONTINUE
RETURN
END
DOUBLE PRECISION FUNCTION FLOP(X)
DOUBLE PRECISION X
C SYSTEM DEPENDENT FUNCTION
C COUNT AND POSSIBLY CHOP EACH FLOATING POINT OPERATION
C FLP(1) IS FLOP COUNTER
C FLP(2) IS NUMBER OF PLACES TO BE CHOPPED
C
INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
C
DOUBLE PRECISION MASK(11),XX,MM
INTEGER*4 MAS(2,11),LX(2),LM(2)
C LOGICAL LX(4),LM(4)
EQUIVALENCE (LX(1),XX),(LM(1),MM),(MASK,MAS)
C
C *** THESE ARE THE MASKS TO CHOP HEX DIGITS ON THE PRIME 400 ***
C
DATA MAS /-:00000000001,-:00003600001,
$ -:00000000001,-:00077600001,
$ -:00000000001,-:01777600001,
$ -:00000000001,-:37777600001,
$ -:00000000020,:00000177777,
$ -:00000000400,:00000177777,
$ -:00000010000,:00000177777,
$ -:00000200000,:00000177777,
$ -:00004000000,:00000177777,
$ -:00100000000,:00000177777,
$ -:02000000000,:00000177777/
FLP(1) = FLP(1) + 1
K = FLP(2)
FLOP = X
IF (K .LE. 0) RETURN
FLOP = 0.0D0
IF (K .GE. 12) RETURN
XX = X
MM = MASK(K)
LX(1) = AND(LX(1),LM(1))
LX(2) = AND(LX(2),LM(2))
FLOP = XX
RETURN
END
SUBROUTINE XCHAR(BUF,K)
INTEGER BUF(1),K
C
C SYSTEM DEPENDENT ROUTINE TO HANDLE SPECIAL CHARACTERS
C
WRITE(6,10) BUF(1)
10 FORMAT(1X,A1,' is not a MATLAB character.')
RETURN
END
SUBROUTINE USER(A,M,N,S,T)
DOUBLE PRECISION A(M,N),S,T
C
INTEGER A3(9)
DATA A3 /-149,537,-27,-50,180,-9,-154,546,-25/
IF (A(1,1) .NE. 3.0D0) RETURN
DO 10 I = 1, 9
A(I,1) = A3(I)
10 CONTINUE
M = 3
N = 3
RETURN
END
SUBROUTINE PROMPT(PAUSE)
INTEGER PAUSE
C
C ISSUE MATLAB PROMPT WITH OPTIONAL PAUSE
C
INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
WRITE(WTE,10)
IF (WIO .NE. 0) WRITE(WIO,10)
10 FORMAT(/1X,'<>')
IF (PAUSE .EQ. 1) READ(RTE,20) DUMMY
20 FORMAT(A1)
RETURN
END
SUBROUTINE PLOT(LUNIT,X,Y,N,P,K,BUF)
DOUBLE PRECISION X(N),Y(N),P(1)
INTEGER BUF(79)
C
C PLOT X VS. Y ON LUNIT
C IF K IS NONZERO, THEN P(1),...,P(K) ARE EXTRA PARAMETERS
C BUF IS WORK SPACE
C
DOUBLE PRECISION XMIN,YMIN,XMAX,YMAX,DY,DX,Y1,Y0
INTEGER AST,BLANK,H,W
DATA AST/1H*/,BLANK/1H /,H/20/,W/79/
C
C H = HEIGHT, W = WIDTH
C
XMIN = X(1)
XMAX = X(1)
YMIN = Y(1)
YMAX = Y(1)
DO 10 I = 1, N
XMIN = DMIN1(XMIN,X(I))
XMAX = DMAX1(XMAX,X(I))
YMIN = DMIN1(YMIN,Y(I))
YMAX = DMAX1(YMAX,Y(I))
10 CONTINUE
DX = XMAX - XMIN
IF (DX .EQ. 0.0D0) DX = 1.0D0
DY = YMAX - YMIN
WRITE(LUNIT,35)
DO 40 L = 1, H
DO 20 J = 1, W
BUF(J) = BLANK
20 CONTINUE
Y1 = YMIN + (H-L+1)*DY/H
Y0 = YMIN + (H-L)*DY/H
JMAX = 1
DO 30 I = 1, N
IF (Y(I) .GT. Y1) GO TO 30
IF (L.NE.H .AND. Y(I).LE.Y0) GO TO 30
J = 1 + (W-1)*(X(I) - XMIN)/DX
BUF(J) = AST
JMAX = MAX0(JMAX,J)
30 CONTINUE
WRITE(LUNIT,35) (BUF(J),J=1,JMAX)
35 FORMAT(1X,79A1)
40 CONTINUE
RETURN
END
SUBROUTINE EDIT(BUF,N)
INTEGER BUF(N)
C
C CALLED AFTER INPUT OF A SINGLE BACKSLASH
C BUF CONTAINS PREVIOUS INPUT LINE, ONE CHAR PER WORD
C ENTER LOCAL EDITOR IF AVAILABLE
C OTHERWISE JUST
RETURN
END
DOUBLE PRECISION FUNCTION DSINH(X)
IMPLICIT DOUBLE PRECISION (A-Z)
DATA A0/-.139005324307231D6/,
$ A1/-.180697428204813D5/,
$ A2/-.404150536907816D3/,
$ B0/-.139005324307231D6/,
$ B1/+.509781123072247D4/,
$ B2/-.954080394016126D2/
IF (DABS(X) .LT. 0.5D0) GO TO 10
E = DEXP(X)
DSINH = 0.5D0*(E-1.0D0/E)
RETURN
10 Y = X*X
P = ((A2*Y)+A1)*Y+A0
Q = ((Y+B2)*Y+B1)*Y+B0
DSINH = X*P/Q
RETURN
END
DOUBLE PRECISION FUNCTION DCOSH(X)
IMPLICIT DOUBLE PRECISION (A-Z)
E = DEXP(X)
DCOSH = 0.5D0*(E+1.0D0/E)
RETURN
END