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 42
Archive-name: applications/matlab/src.2
# 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-2
# This archive created: Wed Nov 2 16:20:14 1988
cat << \SHAR_EOF > src-2
DOUBLE PRECISION SYV,S,FLOP
INTEGER BLANK,Z,DOT,D,E,PLUS,MINUS,NAME,NUM,SIGN,CHCNT,EOL
INTEGER STAR,SLASH,BSLASH,SS
DATA BLANK/36/,Z/35/,DOT/47/,D/13/,E/14/,EOL/99/,PLUS/41/
DATA MINUS/42/,NAME/1/,NUM/0/,STAR/43/,SLASH/44/,BSLASH/45/
10 IF (CHAR .NE. BLANK) GO TO 20
CALL GETCH
GO TO 10
20 LPT(2) = LPT(3)
LPT(3) = LPT(4)
IF (CHAR .LE. 9) GO TO 50
IF (CHAR .LE. Z) GO TO 30
C
C SPECIAL CHARACTER
SS = SYM
SYM = CHAR
CALL GETCH
IF (SYM .NE. DOT) GO TO 90
C
C IS DOT PART OF NUMBER OR OPERATOR
SYV = 0.0D0
IF (CHAR .LE. 9) GO TO 55
IF (CHAR.EQ.STAR .OR. CHAR.EQ.SLASH .OR. CHAR.EQ.BSLASH) GO TO 90
IF (SS.EQ.STAR .OR. SS.EQ.SLASH .OR. SS.EQ.BSLASH) GO TO 90
GO TO 55
C
C NAME
30 SYM = NAME
SYN(1) = CHAR
CHCNT = 1
40 CALL GETCH
CHCNT = CHCNT+1
IF (CHAR .GT. Z) GO TO 45
IF (CHCNT .LE. 4) SYN(CHCNT) = CHAR
GO TO 40
45 IF (CHCNT .GT. 4) GO TO 47
DO 46 I = CHCNT, 4
46 SYN(I) = BLANK
47 CONTINUE
GO TO 90
C
C NUMBER
50 CALL GETVAL(SYV)
IF (CHAR .NE. DOT) GO TO 60
CALL GETCH
55 CHCNT = LPT(4)
CALL GETVAL(S)
CHCNT = LPT(4) - CHCNT
IF (CHAR .EQ. EOL) CHCNT = CHCNT+1
SYV = SYV + S/10.0D0**CHCNT
60 IF (CHAR.NE.D .AND. CHAR.NE.E) GO TO 70
CALL GETCH
SIGN = CHAR
IF (SIGN.EQ.MINUS .OR. SIGN.EQ.PLUS) CALL GETCH
CALL GETVAL(S)
IF (SIGN .NE. MINUS) SYV = SYV*10.0D0**S
IF (SIGN .EQ. MINUS) SYV = SYV/10.0D0**S
70 STKI(VSIZE) = FLOP(SYV)
SYM = NUM
C
90 IF (CHAR .NE. BLANK) GO TO 99
CALL GETCH
GO TO 90
99 IF (DDT .NE. 1) RETURN
IF (SYM.GT.NAME .AND. SYM.LT.ALFL) WRITE(WTE,197) ALFA(SYM+1)
IF (SYM .GE. ALFL) WRITE(WTE,198)
IF (SYM .EQ. NAME) CALL PRNTID(SYN,1)
IF (SYM .EQ. NUM) WRITE(WTE,199) SYV
197 FORMAT(1X,A1)
198 FORMAT(1X,'EOL')
199 FORMAT(1X,G8.2)
RETURN
END
SUBROUTINE GETVAL(S)
DOUBLE PRECISION S
C FORM NUMERICAL VALUE
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
S = 0.0D0
10 IF (CHAR .GT. 9) RETURN
S = 10.0D0*S + CHAR
CALL GETCH
GO TO 10
END
SUBROUTINE MATFN1
C
C EVALUATE FUNCTIONS INVOLVING GAUSSIAN ELIMINATION
C
DOUBLE PRECISION STKR(5005),STKI(5005)
INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE
INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE
COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
DOUBLE PRECISION DTR(2),DTI(2),SR,SI,RCOND,T,T0,T1,FLOP,EPS,WASUM
C
IF (DDT .EQ. 1) WRITE(WTE,100) FIN
100 FORMAT(1X,'MATFN1',I4)
C
L = LSTK(TOP)
M = MSTK(TOP)
N = NSTK(TOP)
IF (FIN .EQ. -1) GO TO 10
IF (FIN .EQ. -2) GO TO 20
GO TO (30,40,50,60,70,80,85),FIN
C
C MATRIX RIGHT DIVISION, A/A2
10 L2 = LSTK(TOP+1)
M2 = MSTK(TOP+1)
N2 = NSTK(TOP+1)
IF (M2 .NE. N2) CALL ERROR(20)
IF (ERR .GT. 0) RETURN
IF (M*N .EQ. 1) GO TO 16
IF (N .NE. N2) CALL ERROR(11)
IF (ERR .GT. 0) RETURN
L3 = L2 + M2*N2
ERR = L3+N2 - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
CALL WGECO(STKR(L2),STKI(L2),M2,N2,BUF,RCOND,STKR(L3),STKI(L3))
IF (RCOND .EQ. 0.0D0) CALL ERROR(19)
IF (ERR .GT. 0) RETURN
T = FLOP(1.0D0 + RCOND)
IF (T.EQ.1.0D0 .AND. FUN.NE.21) WRITE(WTE,11) RCOND
IF (T.EQ.1.0D0 .AND. FUN.NE.21 .AND. WIO.NE.0) WRITE(WIO,11) RCOND
11 FORMAT(1X,'WARNING.'
$ /1X,'MATRIX IS CLOSE TO SINGULAR OR BADLY SCALED.'
$ /1X,'RESULTS MAY BE INACCURATE. RCOND =', 1PD13.4/)
IF (T.EQ.1.0D0 .AND. FUN.EQ.21) WRITE(WTE,12) RCOND
IF (T.EQ.1.0D0 .AND. FUN.EQ.21 .AND. WIO.NE.0) WRITE(WIO,12) RCOND
12 FORMAT(1X,'WARNING.'
$ /1X,'EIGENVECTORS ARE BADLY CONDITIONED.'
$ /1X,'RESULTS MAY BE INACCURATE. RCOND =', 1PD13.4/)
DO 15 I = 1, M
DO 13 J = 1, N
LS = L+I-1+(J-1)*M
LL = L3+J-1
STKR(LL) = STKR(LS)
STKI(LL) = -STKI(LS)
13 CONTINUE
CALL WGESL(STKR(L2),STKI(L2),M2,N2,BUF,STKR(L3),STKI(L3),1)
DO 14 J = 1, N
LL = L+I-1+(J-1)*M
LS = L3+J-1
STKR(LL) = STKR(LS)
STKI(LL) = -STKI(LS)
14 CONTINUE
15 CONTINUE
IF (FUN .NE. 21) GO TO 99
C
C CHECK FOR IMAGINARY ROUNDOFF IN MATRIX FUNCTIONS
SR = WASUM(N*N,STKR(L),STKR(L),1)
SI = WASUM(N*N,STKI(L),STKI(L),1)
EPS = STKR(VSIZE-4)
T = EPS*SR
IF (DDT .EQ. 18) WRITE(WTE,115) SR,SI,EPS,T
115 FORMAT(1X,'SR,SI,EPS,T',1P4D13.4)
IF (SI .LE. EPS*SR) CALL RSET(N*N,0.0D0,STKI(L),1)
GO TO 99
C
16 SR = STKR(L)
SI = STKI(L)
N = N2
M = N
MSTK(TOP) = N
NSTK(TOP) = N
CALL WCOPY(N*N,STKR(L2),STKI(L2),1,STKR(L),STKI(L),1)
GO TO 30
C
C MATRIX LEFT DIVISION A BACKSLASH A2
20 L2 = LSTK(TOP+1)
M2 = MSTK(TOP+1)
N2 = NSTK(TOP+1)
IF (M .NE. N) CALL ERROR(20)
IF (ERR .GT. 0) RETURN
IF (M2*N2 .EQ. 1) GO TO 26
L3 = L2 + M2*N2
ERR = L3+N - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
CALL WGECO(STKR(L),STKI(L),M,N,BUF,RCOND,STKR(L3),STKI(L3))
IF (RCOND .EQ. 0.0D0) CALL ERROR(19)
IF (ERR .GT. 0) RETURN
T = FLOP(1.0D0 + RCOND)
IF (T .EQ. 1.0D0) WRITE(WTE,11) RCOND
IF (T.EQ.1.0D0 .AND. WIO.NE.0) WRITE(WIO,11) RCOND
IF (M2 .NE. N) CALL ERROR(12)
IF (ERR .GT. 0) RETURN
DO 23 J = 1, N2
LJ = L2+(J-1)*M2
CALL WGESL(STKR(L),STKI(L),M,N,BUF,STKR(LJ),STKI(LJ),0)
23 CONTINUE
NSTK(TOP) = N2
CALL WCOPY(M2*N2,STKR(L2),STKI(L2),1,STKR(L),STKI(L),1)
GO TO 99
26 SR = STKR(L2)
SI = STKI(L2)
GO TO 30
C
C INV
C
30 IF (M .NE. N) CALL ERROR(20)
IF (ERR .GT. 0) RETURN
IF (DDT .EQ. 17) GO TO 32
DO 31 J = 1, N
DO 31 I = 1, N
LS = L+I-1+(J-1)*N
T0 = STKR(LS)
T1 = FLOP(1.0D0/(DFLOAT(I+J-1)))
IF (T0 .NE. T1) GO TO 32
31 CONTINUE
GO TO 72
32 L3 = L + N*N
ERR = L3+N - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
CALL WGECO(STKR(L),STKI(L),M,N,BUF,RCOND,STKR(L3),STKI(L3))
IF (RCOND .EQ. 0.0D0) CALL ERROR(19)
IF (ERR .GT. 0) RETURN
T = FLOP(1.0D0 + RCOND)
IF (T .EQ. 1.0D0) WRITE(WTE,11) RCOND
IF (T.EQ.1.0D0 .AND. WIO.NE.0) WRITE(WIO,11) RCOND
CALL WGEDI(STKR(L),STKI(L),M,N,BUF,DTR,DTI,STKR(L3),STKI(L3),1)
IF (FIN .LT. 0) CALL WSCAL(N*N,SR,SI,STKR(L),STKI(L),1)
GO TO 99
C
C DET
C
40 IF (M .NE. N) CALL ERROR(20)
IF (ERR .GT. 0) RETURN
CALL WGEFA(STKR(L),STKI(L),M,N,BUF,INFO)
CALL WGEDI(STKR(L),STKI(L),M,N,BUF,DTR,DTI,SR,SI,10)
K = IDINT(DTR(2))
KA = IABS(K)+2
T = 1.0D0
DO 41 I = 1, KA
T = T/10.0D0
IF (T .EQ. 0.0D0) GO TO 42
41 CONTINUE
STKR(L) = DTR(1)*10.D0**K
STKI(L) = DTI(1)*10.D0**K
MSTK(TOP) = 1
NSTK(TOP) = 1
GO TO 99
42 IF (DTI(1) .EQ. 0.0D0) WRITE(WTE,43) DTR(1),K
IF (DTI(1) .NE. 0.0D0) WRITE(WTE,44) DTR(1),DTI(1),K
43 FORMAT(1X,'DET = ',F7.4,7H * 10**,I4)
44 FORMAT(1X,'DET = ',F7.4,' + ',F7.4,' i ',7H * 10**,I4)
STKR(L) = DTR(1)
STKI(L) = DTI(1)
STKR(L+1) = DTR(2)
STKI(L+1) = 0.0D0
MSTK(TOP) = 1
NSTK(TOP) = 2
GO TO 99
C
C RCOND
C
50 IF (M .NE. N) CALL ERROR(20)
IF (ERR .GT. 0) RETURN
L3 = L + N*N
ERR = L3+N - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
CALL WGECO(STKR(L),STKI(L),M,N,BUF,RCOND,STKR(L3),STKI(L3))
STKR(L) = RCOND
STKI(L) = 0.0D0
MSTK(TOP) = 1
NSTK(TOP) = 1
IF (LHS .EQ. 1) GO TO 99
L = L + 1
CALL WCOPY(N,STKR(L3),STKI(L3),1,STKR(L),STKI(L),1)
TOP = TOP + 1
LSTK(TOP) = L
MSTK(TOP) = N
NSTK(TOP) = 1
GO TO 99
C
C LU
C
60 IF (M .NE. N) CALL ERROR(20)
IF (ERR .GT. 0) RETURN
CALL WGEFA(STKR(L),STKI(L),M,N,BUF,INFO)
IF (LHS .NE. 2) GO TO 99
NN = N*N
IF (TOP+1 .GE. BOT) CALL ERROR(18)
IF (ERR .GT. 0) RETURN
TOP = TOP+1
LSTK(TOP) = L + NN
MSTK(TOP) = N
NSTK(TOP) = N
ERR = L+NN+NN - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
DO 64 KB = 1, N
K = N+1-KB
DO 61 I = 1, N
LL = L+I-1+(K-1)*N
LU = LL + NN
IF (I .LE. K) STKR(LU) = STKR(LL)
IF (I .LE. K) STKI(LU) = STKI(LL)
IF (I .GT. K) STKR(LU) = 0.0D0
IF (I .GT. K) STKI(LU) = 0.0D0
IF (I .LT. K) STKR(LL) = 0.0D0
IF (I .LT. K) STKI(LL) = 0.0D0
IF (I .EQ. K) STKR(LL) = 1.0D0
IF (I .EQ. K) STKI(LL) = 0.0D0
IF (I .GT. K) STKR(LL) = -STKR(LL)
IF (I .GT. K) STKI(LL) = -STKI(LL)
61 CONTINUE
I = BUF(K)
IF (I .EQ. K) GO TO 64
LI = L+I-1+(K-1)*N
LK = L+K-1+(K-1)*N
CALL WSWAP(N-K+1,STKR(LI),STKI(LI),N,STKR(LK),STKI(LK),N)
64 CONTINUE
GO TO 99
C
C HILBERT
70 N = IDINT(STKR(L))
MSTK(TOP) = N
NSTK(TOP) = N
72 CALL HILBER(STKR(L),N,N)
CALL RSET(N*N,0.0D0,STKI(L),1)
IF (FIN .LT. 0) CALL WSCAL(N*N,SR,SI,STKR(L),STKI(L),1)
GO TO 99
C
C CHOLESKY
80 IF (M .NE. N) CALL ERROR(20)
IF (ERR .GT. 0) RETURN
CALL WPOFA(STKR(L),STKI(L),M,N,ERR)
IF (ERR .NE. 0) CALL ERROR(29)
IF (ERR .GT. 0) RETURN
DO 81 J = 1, N
LL = L+J+(J-1)*M
CALL WSET(M-J,0.0D0,0.0D0,STKR(LL),STKI(LL),1)
81 CONTINUE
GO TO 99
C
C RREF
85 IF (RHS .LT. 2) GO TO 86
TOP = TOP-1
L = LSTK(TOP)
IF (MSTK(TOP) .NE. M) CALL ERROR(5)
IF (ERR .GT. 0) RETURN
N = N + NSTK(TOP)
86 CALL RREF(STKR(L),STKI(L),M,M,N,STKR(VSIZE-4))
NSTK(TOP) = N
GO TO 99
C
99 RETURN
END
SUBROUTINE MATFN2
C
C EVALUATE ELEMENTARY FUNCTIONS AND FUNCTIONS INVOLVING
C EIGENVALUES AND EIGENVECTORS
C
DOUBLE PRECISION STKR(5005),STKI(5005)
INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE
INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE
COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
DOUBLE PRECISION PYTHAG,ROUND,TR,TI,SR,SI,POWR,POWI,FLOP
LOGICAL HERM,SCHUR,VECT,HESS
C
IF (DDT .EQ. 1) WRITE(WTE,100) FIN
100 FORMAT(1X,'MATFN2',I4)
C
C FUNCTIONS/FIN
C ** SIN COS ATAN EXP SQRT LOG
C 0 1 2 3 4 5 6
C EIG SCHU HESS POLY ROOT
C 11 12 13 14 15
C ABS ROUN REAL IMAG CONJ
C 21 22 23 24 25
IF (FIN .NE. 0) GO TO 05
L = LSTK(TOP+1)
POWR = STKR(L)
POWI = STKI(L)
05 L = LSTK(TOP)
M = MSTK(TOP)
N = NSTK(TOP)
IF (FIN .GE. 11 .AND. FIN .LE. 13) GO TO 10
IF (FIN .EQ. 14 .AND. (M.EQ.1 .OR. N.EQ.1)) GO TO 50
IF (FIN .EQ. 14) GO TO 10
IF (FIN .EQ. 15) GO TO 60
IF (FIN .GT. 20) GO TO 40
IF (M .EQ. 1 .OR. N .EQ. 1) GO TO 40
C
C EIGENVALUES AND VECTORS
10 IF (M .NE. N) CALL ERROR(20)
IF (ERR .GT. 0) RETURN
SCHUR = FIN .EQ. 12
HESS = FIN .EQ. 13
VECT = LHS.EQ.2 .OR. FIN.LT.10
NN = N*N
L2 = L + NN
LD = L2 + NN
LE = LD + N
LW = LE + N
ERR = LW+N - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
CALL WCOPY(NN,STKR(L),STKI(L),1,STKR(L2),STKI(L2),1)
C
C CHECK IF HERMITIAN
DO 15 J = 1, N
DO 15 I = 1, J
LS = L+I-1+(J-1)*N
LL = L+(I-1)*N+J-1
HERM = STKR(LL).EQ.STKR(LS) .AND. STKI(LL).EQ.-STKI(LS)
IF (.NOT. HERM) GO TO 30
15 CONTINUE
C
C HERMITIAN EIGENVALUE PROBLEM
CALL WSET(NN,0.0D0,0.0D0,STKR(L),STKI(L),1)
CALL WSET(N,1.0D0,0.0D0,STKR(L),STKI(L),N+1)
CALL WSET(N,0.0D0,0.0D0,STKI(LD),STKI(LE),1)
JOB = 0
IF (VECT) JOB = 1
CALL HTRIDI(N,N,STKR(L2),STKI(L2),STKR(LD),STKR(LE),
$ STKR(LE),STKR(LW))
IF (.NOT.HESS) CALL IMTQL2(N,N,STKR(LD),STKR(LE),STKR(L),ERR,JOB)
IF (ERR .GT. 0) CALL ERROR(24)
IF (ERR .GT. 0) RETURN
IF (JOB .NE. 0)
$ CALL HTRIBK(N,N,STKR(L2),STKI(L2),STKR(LW),N,STKR(L),STKI(L))
GO TO 31
C
C NON-HERMITIAN EIGENVALUE PROBLEM
30 CALL CORTH(N,N,1,N,STKR(L2),STKI(L2),STKR(LW),STKI(LW))
IF (.NOT.VECT .AND. HESS) GO TO 31
JOB = 0
IF (VECT) JOB = 2
IF (VECT .AND. SCHUR) JOB = 1
IF (HESS) JOB = 3
CALL COMQR3(N,N,1,N,STKR(LW),STKI(LW),STKR(L2),STKI(L2),
$ STKR(LD),STKI(LD),STKR(L),STKI(L),ERR,JOB)
IF (ERR .GT. 0) CALL ERROR(24)
IF (ERR .GT. 0) RETURN
C
C VECTORS
31 IF (.NOT.VECT) GO TO 34
IF (TOP+1 .GE. BOT) CALL ERROR(18)
IF (ERR .GT. 0) RETURN
TOP = TOP+1
LSTK(TOP) = L2
MSTK(TOP) = N
NSTK(TOP) = N
C
C DIAGONAL OF VALUES OR CANONICAL FORMS
34 IF (.NOT.VECT .AND. .NOT.SCHUR .AND. .NOT.HESS) GO TO 37
DO 36 J = 1, N
LJ = L2+(J-1)*N
IF (SCHUR .AND. (.NOT.HERM)) LJ = LJ+J
IF (HESS .AND. (.NOT.HERM)) LJ = LJ+J+1
LL = L2+J*N-LJ
CALL WSET(LL,0.0D0,0.0D0,STKR(LJ),STKI(LJ),1)
36 CONTINUE
IF (.NOT.HESS .OR. HERM)
$ CALL WCOPY(N,STKR(LD),STKI(LD),1,STKR(L2),STKI(L2),N+1)
LL = L2+1
IF (HESS .AND. HERM)
$ CALL WCOPY(N-1,STKR(LE+1),STKI(LE+1),1,STKR(LL),STKI(LL),N+1)
LL = L2+N
IF (HESS .AND. HERM)
$ CALL WCOPY(N-1,STKR(LE+1),STKI(LE+1),1,STKR(LL),STKI(LL),N+1)
IF (FIN .LT. 10) GO TO 42
IF (VECT .OR. .NOT.(SCHUR.OR.HESS)) GO TO 99
CALL WCOPY(NN,STKR(L2),STKI(L2),1,STKR(L),STKI(L),1)
GO TO 99
C
C VECTOR OF EIGENVALUES
37 IF (FIN .EQ. 14) GO TO 52
CALL WCOPY(N,STKR(LD),STKI(LD),1,STKR(L),STKI(L),1)
NSTK(TOP) = 1
GO TO 99
C
C ELEMENTARY FUNCTIONS
C FOR MATRICES.. X,D = EIG(A), FUN(A) = X*FUN(D)/X
40 INC = 1
N = M*N
L2 = L
GO TO 44
42 INC = N+1
44 DO 46 J = 1, N
LS = L2+(J-1)*INC
SR = STKR(LS)
SI = STKI(LS)
TI = 0.0D0
IF (FIN .NE. 0) GO TO 45
CALL WLOG(SR,SI,SR,SI)
CALL WMUL(SR,SI,POWR,POWI,SR,SI)
TR = DEXP(SR)*DCOS(SI)
TI = DEXP(SR)*DSIN(SI)
45 IF (FIN .EQ. 1) TR = DSIN(SR)*DCOSH(SI)
IF (FIN .EQ. 1) TI = DCOS(SR)*DSINH(SI)
IF (FIN .EQ. 2) TR = DCOS(SR)*DCOSH(SI)
IF (FIN .EQ. 2) TI = -DSIN(SR)*DSINH(SI)
IF (FIN .EQ. 3) CALL WATAN(SR,SI,TR,TI)
IF (FIN .EQ. 4) TR = DEXP(SR)*DCOS(SI)
IF (FIN .EQ. 4) TI = DEXP(SR)*DSIN(SI)
IF (FIN .EQ. 5) CALL WSQRT(SR,SI,TR,TI)
IF (FIN .EQ. 6) CALL WLOG(SR,SI,TR,TI)
IF (FIN .EQ. 21) TR = PYTHAG(SR,SI)
IF (FIN .EQ. 22) TR = ROUND(SR)
IF (FIN .EQ. 23) TR = SR
IF (FIN .EQ. 24) TR = SI
IF (FIN .EQ. 25) TR = SR
IF (FIN .EQ. 25) TI = -SI
IF (ERR .GT. 0) RETURN
STKR(LS) = FLOP(TR)
STKI(LS) = 0.0D0
IF (TI .NE. 0.0D0) STKI(LS) = FLOP(TI)
46 CONTINUE
IF (INC .EQ. 1) GO TO 99
DO 48 J = 1, N
LS = L2+(J-1)*INC
SR = STKR(LS)
SI = STKI(LS)
LS = L+(J-1)*N
LL = L2+(J-1)*N
CALL WCOPY(N,STKR(LS),STKI(LS),1,STKR(LL),STKI(LL),1)
CALL WSCAL(N,SR,SI,STKR(LS),STKI(LS),1)
48 CONTINUE
C SIGNAL MATFN1 TO DIVIDE BY EIGENVECTORS
FUN = 21
FIN = -1
TOP = TOP-1
GO TO 99
C
C POLY
C FORM POLYNOMIAL WITH GIVEN VECTOR AS ROOTS
50 N = MAX0(M,N)
LD = L+N+1
CALL WCOPY(N,STKR(L),STKI(L),1,STKR(LD),STKI(LD),1)
C
C FORM CHARACTERISTIC POLYNOMIAL
52 CALL WSET(N+1,0.0D0,0.0D0,STKR(L),STKI(L),1)
STKR(L) = 1.0D0
DO 56 J = 1, N
CALL WAXPY(J,-STKR(LD),-STKI(LD),STKR(L),STKI(L),-1,
$ STKR(L+1),STKI(L+1),-1)
LD = LD+1
56 CONTINUE
MSTK(TOP) = N+1
NSTK(TOP) = 1
GO TO 99
C
C ROOTS
60 LL = L+M*N
STKR(LL) = -1.0D0
STKI(LL) = 0.0D0
K = -1
61 K = K+1
L1 = L+K
IF (DABS(STKR(L1))+DABS(STKI(L1)) .EQ. 0.0D0) GO TO 61
N = MAX0(M*N - K-1, 0)
IF (N .LE. 0) GO TO 65
L2 = L1+N+1
LW = L2+N*N
ERR = LW+N - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
CALL WSET(N*N+N,0.0D0,0.0D0,STKR(L2),STKI(L2),1)
DO 64 J = 1, N
LL = L2+J+(J-1)*N
STKR(LL) = 1.0D0
LS = L1+J
LL = L2+(J-1)*N
CALL WDIV(-STKR(LS),-STKI(LS),STKR(L1),STKI(L1),
$ STKR(LL),STKI(LL))
IF (ERR .GT. 0) RETURN
64 CONTINUE
CALL COMQR3(N,N,1,N,STKR(LW),STKI(LW),STKR(L2),STKI(L2),
$ STKR(L),STKI(L),TR,TI,ERR,0)
IF (ERR .GT. 0) CALL ERROR(24)
IF (ERR .GT. 0) RETURN
65 MSTK(TOP) = N
NSTK(TOP) = 1
GO TO 99
99 RETURN
END
SUBROUTINE MATFN3
C
C EVALUATE FUNCTIONS INVOLVING SINGULAR VALUE DECOMPOSITION
C
DOUBLE PRECISION STKR(5005),STKI(5005)
INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE
INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE
COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
LOGICAL FRO,INF
DOUBLE PRECISION P,S,T,TOL,EPS
DOUBLE PRECISION WDOTCR,WDOTCI,PYTHAG,WNRM2,WASUM,FLOP
C
IF (DDT .EQ. 1) WRITE(WTE,100) FIN
100 FORMAT(1X,'MATFN3',I4)
C
IF (FIN.EQ.1 .AND. RHS.EQ.2) TOP = TOP-1
L = LSTK(TOP)
M = MSTK(TOP)
N = NSTK(TOP)
MN = M*N
GO TO (50,70,10,30,70), FIN
C
C COND
C
10 LD = L + M*N
L1 = LD + MIN0(M+1,N)
L2 = L1 + N
ERR = L2+MIN0(M,N) - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
CALL WSVDC(STKR(L),STKI(L),M,M,N,STKR(LD),STKI(LD),
$ STKR(L1),STKI(L1),T,T,1,T,T,1,STKR(L2),STKI(L2),
$ 0,ERR)
IF (ERR .NE. 0) CALL ERROR(24)
IF (ERR .GT. 0) RETURN
S = STKR(LD)
LD = LD + MIN0(M,N) - 1
T = STKR(LD)
IF (T .EQ. 0.0D0) GO TO 13
STKR(L) = FLOP(S/T)
STKI(L) = 0.0D0
MSTK(TOP) = 1
NSTK(TOP) = 1
GO TO 99
13 WRITE(WTE,14)
IF (WIO .NE. 0) WRITE(WIO,14)
14 FORMAT(1X,'CONDITION IS INFINITE')
MSTK(TOP) = 0
GO TO 99
C
C NORM
C
30 P = 2.0D0
INF = .FALSE.
IF (RHS .NE. 2) GO TO 31
FRO = IDINT(STKR(L)).EQ.15 .AND. MN.GT.1
INF = IDINT(STKR(L)).EQ.18 .AND. MN.GT.1
IF (.NOT. FRO) P = STKR(L)
TOP = TOP-1
L = LSTK(TOP)
M = MSTK(TOP)
N = NSTK(TOP)
MN = M*N
IF (FRO) M = MN
IF (FRO) N = 1
31 IF (M .GT. 1 .AND. N .GT. 1) GO TO 40
IF (P .EQ. 1.0D0) GO TO 36
IF (P .EQ. 2.0D0) GO TO 38
I = IWAMAX(MN,STKR(L),STKI(L),1) + L - 1
S = DABS(STKR(I)) + DABS(STKI(I))
IF (INF .OR. S .EQ. 0.0D0) GO TO 49
T = 0.0D0
DO 33 I = 1, MN
LS = L+I-1
T = FLOP(T + (PYTHAG(STKR(LS),STKI(LS))/S)**P)
33 CONTINUE
IF (P .NE. 0.0D0) P = 1.0D0/P
S = FLOP(S*T**P)
GO TO 49
36 S = WASUM(MN,STKR(L),STKI(L),1)
GO TO 49
38 S = WNRM2(MN,STKR(L),STKI(L),1)
GO TO 49
C
C MATRIX NORM
C
40 IF (INF) GO TO 43
IF (P .EQ. 1.0D0) GO TO 46
IF (P .NE. 2.0D0) CALL ERROR(23)
IF (ERR .GT. 0) RETURN
LD = L + M*N
L1 = LD + MIN0(M+1,N)
L2 = L1 + N
ERR = L2+MIN0(M,N) - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
CALL WSVDC(STKR(L),STKI(L),M,M,N,STKR(LD),STKI(LD),
$ STKR(L1),STKI(L1),T,T,1,T,T,1,STKR(L2),STKI(L2),
$ 0,ERR)
IF (ERR .NE. 0) CALL ERROR(24)
IF (ERR .GT. 0) RETURN
S = STKR(LD)
GO TO 49
43 S = 0.0D0
DO 45 I = 1, M
LI = L+I-1
T = WASUM(N,STKR(LI),STKI(LI),M)
S = DMAX1(S,T)
45 CONTINUE
GO TO 49
46 S = 0.0D0
DO 48 J = 1, N
LJ = L+(J-1)*M
T = WASUM(M,STKR(LJ),STKI(LJ),1)
S = DMAX1(S,T)
48 CONTINUE
GO TO 49
49 STKR(L) = S
STKI(L) = 0.0D0
MSTK(TOP) = 1
NSTK(TOP) = 1
GO TO 99
C
C SVD
C
50 IF (LHS .NE. 3) GO TO 52
K = M
IF (RHS .EQ. 2) K = MIN0(M,N)
LU = L + M*N
LD = LU + M*K
LV = LD + K*N
L1 = LV + N*N
L2 = L1 + N
ERR = L2+MIN0(M,N) - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
JOB = 11
IF (RHS .EQ. 2) JOB = 21
CALL WSVDC(STKR(L),STKI(L),M,M,N,STKR(LD),STKI(LD),
$ STKR(L1),STKI(L1),STKR(LU),STKI(LU),M,STKR(LV),STKI(LV),
$ N,STKR(L2),STKI(L2),JOB,ERR)
DO 51 JB = 1, N
DO 51 I = 1, K
J = N+1-JB
LL = LD+I-1+(J-1)*K
IF (I.NE.J) STKR(LL) = 0.0D0
STKI(LL) = 0.0D0
LS = LD+I-1
IF (I.EQ.J) STKR(LL) = STKR(LS)
LS = L1+I-1
IF (ERR.NE.0 .AND. I.EQ.J-1) STKR(LL) = STKR(LS)
51 CONTINUE
IF (ERR .NE. 0) CALL ERROR(24)
ERR = 0
CALL WCOPY(M*K+K*N+N*N,STKR(LU),STKI(LU),1,STKR(L),STKI(L),1)
MSTK(TOP) = M
NSTK(TOP) = K
IF (TOP+1 .GE. BOT) CALL ERROR(18)
IF (ERR .GT. 0) RETURN
TOP = TOP+1
LSTK(TOP) = L + M*K
MSTK(TOP) = K
NSTK(TOP) = N
IF (TOP+1 .GE. BOT) CALL ERROR(18)
IF (ERR .GT. 0) RETURN
TOP = TOP+1
LSTK(TOP) = L + M*K + K*N
MSTK(TOP) = N
NSTK(TOP) = N
GO TO 99
C
52 LD = L + M*N
L1 = LD + MIN0(M+1,N)
L2 = L1 + N
ERR = L2+MIN0(M,N) - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
CALL WSVDC(STKR(L),STKI(L),M,M,N,STKR(LD),STKI(LD),
$ STKR(L1),STKI(L1),T,T,1,T,T,1,STKR(L2),STKI(L2),
$ 0,ERR)
IF (ERR .NE. 0) CALL ERROR(24)
IF (ERR .GT. 0) RETURN
K = MIN0(M,N)
CALL WCOPY(K,STKR(LD),STKI(LD),1,STKR(L),STKI(L),1)
MSTK(TOP) = K
NSTK(TOP) = 1
GO TO 99
C
C PINV AND RANK
C
70 TOL = -1.0D0
IF (RHS .NE. 2) GO TO 71
TOL = STKR(L)
TOP = TOP-1
L = LSTK(TOP)
M = MSTK(TOP)
N = NSTK(TOP)
71 LU = L + M*N
LD = LU + M*M
IF (FIN .EQ. 5) LD = L + M*N
LV = LD + M*N
L1 = LV + N*N
IF (FIN .EQ. 5) L1 = LD + N
L2 = L1 + N
ERR = L2+MIN0(M,N) - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
IF (FIN .EQ. 2) JOB = 11
IF (FIN .EQ. 5) JOB = 0
CALL WSVDC(STKR(L),STKI(L),M,M,N,STKR(LD),STKI(LD),
$ STKR(L1),STKI(L1),STKR(LU),STKI(LU),M,STKR(LV),STKI(LV),
$ N,STKR(L2),STKI(L2),JOB,ERR)
IF (ERR .NE. 0) CALL ERROR(24)
IF (ERR .GT. 0) RETURN
EPS = STKR(VSIZE-4)
IF (TOL .LT. 0.0D0) TOL = FLOP(DFLOAT(MAX0(M,N))*EPS*STKR(LD))
MN = MIN0(M,N)
K = 0
DO 72 J = 1, MN
LS = LD+J-1
S = STKR(LS)
IF (S .LE. TOL) GO TO 73
K = J
LL = LV+(J-1)*N
IF (FIN .EQ. 2) CALL WRSCAL(N,1.0D0/S,STKR(LL),STKI(LL),1)
72 CONTINUE
73 IF (FIN .EQ. 5) GO TO 78
DO 76 J = 1, M
DO 76 I = 1, N
LL = L+I-1+(J-1)*N
L1 = LV+I-1
L2 = LU+J-1
STKR(LL) = WDOTCR(K,STKR(L2),STKI(L2),M,STKR(L1),STKI(L1),N)
STKI(LL) = WDOTCI(K,STKR(L2),STKI(L2),M,STKR(L1),STKI(L1),N)
76 CONTINUE
MSTK(TOP) = N
NSTK(TOP) = M
GO TO 99
78 STKR(L) = DFLOAT(K)
STKI(L) = 0.0D0
MSTK(TOP) = 1
NSTK(TOP) = 1
GO TO 99
C
99 RETURN
END
SUBROUTINE MATFN4
C
C EVALUATE FUNCTIONS INVOLVING QR DECOMPOSITION (LEAST SQUARES)
C
DOUBLE PRECISION STKR(5005),STKI(5005)
INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE
INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE
COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
DOUBLE PRECISION T,TOL,EPS,FLOP
INTEGER QUOTE
DATA QUOTE/49/
C
IF (DDT .EQ. 1) WRITE(WTE,100) FIN
100 FORMAT(1X,'MATFN4',I4)
C
L = LSTK(TOP)
M = MSTK(TOP)
N = NSTK(TOP)
IF (FIN .EQ. -1) GO TO 10
IF (FIN .EQ. -2) GO TO 20
GO TO 40
C
C RECTANGULAR MATRIX RIGHT DIVISION, A/A2
10 L2 = LSTK(TOP+1)
M2 = MSTK(TOP+1)
N2 = NSTK(TOP+1)
TOP = TOP + 1
IF (N.GT.1 .AND. N.NE.N2) CALL ERROR(11)
IF (ERR .GT. 0) RETURN
CALL STACK1(QUOTE)
IF (ERR .GT. 0) RETURN
LL = L2+M2*N2
CALL WCOPY(M*N,STKR(L),STKI(L),1,STKR(LL),STKI(LL),1)
CALL WCOPY(M*N+M2*N2,STKR(L2),STKI(L2),1,STKR(L),STKI(L),1)
LSTK(TOP) = L+M2*N2
MSTK(TOP) = M
NSTK(TOP) = N
CALL STACK1(QUOTE)
IF (ERR .GT. 0) RETURN
TOP = TOP - 1
M = N2
N = M2
GO TO 20
C
C RECTANGULAR MATRIX LEFT DIVISION A BACKSLASH A2
C
20 L2 = LSTK(TOP+1)
M2 = MSTK(TOP+1)
N2 = NSTK(TOP+1)
IF (M2*N2 .GT. 1) GO TO 21
M2 = M
N2 = M
ERR = L2+M*M - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
CALL WSET(M*M-1,0.0D0,0.0D0,STKR(L2+1),STKI(L2+1),1)
CALL WCOPY(M,STKR(L2),STKI(L2),0,STKR(L2),STKI(L2),M+1)
21 IF (M2 .NE. M) CALL ERROR(12)
IF (ERR .GT. 0) RETURN
L3 = L2 + MAX0(M,N)*N2
L4 = L3 + N
ERR = L4 + N - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
IF (M .GT. N) GO TO 23
DO 22 JB = 1, N2
J = N+1-JB
LS = L2 + (J-1)*M
LL = L2 + (J-1)*N
CALL WCOPY(M,STKR(LS),STKI(LS),-1,STKR(LL),STKI(LL),-1)
22 CONTINUE
23 DO 24 J = 1, N
BUF(J) = 0
24 CONTINUE
CALL WQRDC(STKR(L),STKI(L),M,M,N,STKR(L4),STKI(L4),
$ BUF,STKR(L3),STKI(L3),1)
K = 0
EPS = STKR(VSIZE-4)
T = DABS(STKR(L))+DABS(STKI(L))
TOL = FLOP(DFLOAT(MAX0(M,N))*EPS*T)
MN = MIN0(M,N)
DO 27 J = 1, MN
LS = L+J-1+(J-1)*M
T = DABS(STKR(LS)) + DABS(STKI(LS))
IF (T .GT. TOL) K = J
27 CONTINUE
IF (K .LT. MN) WRITE(WTE,28) K,TOL
IF (K.LT.MN .AND. WIO.NE.0) WRITE(WIO,28) K,TOL
28 FORMAT(1X,'RANK DEFICIENT, RANK =',I4,', TOL =',1PD13.4)
MN = MAX0(M,N)
DO 29 J = 1, N2
LS = L2+(J-1)*MN
CALL WQRSL(STKR(L),STKI(L),M,M,K,STKR(L4),STKI(L4),
$ STKR(LS),STKI(LS),T,T,STKR(LS),STKI(LS),
$ STKR(LS),STKI(LS),T,T,T,T,100,INFO)
LL = LS+K
CALL WSET(N-K,0.0D0,0.0D0,STKR(LL),STKI(LL),1)
29 CONTINUE
DO 31 J = 1, N
BUF(J) = -BUF(J)
31 CONTINUE
DO 35 J = 1, N
IF (BUF(J) .GT. 0) GO TO 35
K = -BUF(J)
BUF(J) = K
33 CONTINUE
IF (K .EQ. J) GO TO 34
LS = L2+J-1
LL = L2+K-1
CALL WSWAP(N2,STKR(LS),STKI(LS),MN,STKR(LL),STKI(LL),MN)
BUF(K) = -BUF(K)
K = BUF(K)
GO TO 33
34 CONTINUE
35 CONTINUE
DO 36 J = 1, N2
LS = L2+(J-1)*MN
LL = L+(J-1)*N
CALL WCOPY(N,STKR(LS),STKI(LS),1,STKR(LL),STKI(LL),1)
36 CONTINUE
MSTK(TOP) = N
NSTK(TOP) = N2
IF (FIN .EQ. -1) CALL STACK1(QUOTE)
IF (ERR .GT. 0) RETURN
GO TO 99
C
C QR
C
40 MM = MAX0(M,N)
LS = L + MM*MM
IF (LHS.EQ.1 .AND. FIN.EQ.1) LS = L
LE = LS + M*N
L4 = LE + MM
ERR = L4+MM - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
IF (LS.NE.L) CALL WCOPY(M*N,STKR(L),STKI(L),1,STKR(LS),STKI(LS),1)
JOB = 1
IF (LHS.LT.3) JOB = 0
DO 42 J = 1, N
BUF(J) = 0
42 CONTINUE
CALL WQRDC(STKR(LS),STKI(LS),M,M,N,STKR(L4),STKI(L4),
$ BUF,STKR(LE),STKI(LE),JOB)
IF (LHS.EQ.1 .AND. FIN.EQ.1) GO TO 99
CALL WSET(M*M,0.0D0,0.0D0,STKR(L),STKI(L),1)
CALL WSET(M,1.0D0,0.0D0,STKR(L),STKI(L),M+1)
DO 43 J = 1, M
LL = L+(J-1)*M
CALL WQRSL(STKR(LS),STKI(LS),M,M,N,STKR(L4),STKI(L4),
$ STKR(LL),STKI(LL),STKR(LL),STKI(LL),T,T,
$ T,T,T,T,T,T,10000,INFO)
43 CONTINUE
IF (FIN .EQ. 2) GO TO 99
NSTK(TOP) = M
DO 45 J = 1, N
LL = LS+J+(J-1)*M
CALL WSET(M-J,0.0D0,0.0D0,STKR(LL),STKI(LL),1)
45 CONTINUE
IF (TOP+1 .GE. BOT) CALL ERROR(18)
IF (ERR .GT. 0) RETURN
TOP = TOP+1
LSTK(TOP) = LS
MSTK(TOP) = M
NSTK(TOP) = N
IF (LHS .EQ. 2) GO TO 99
CALL WSET(N*N,0.0D0,0.0D0,STKR(LE),STKI(LE),1)
DO 47 J = 1, N
LL = LE+BUF(J)-1+(J-1)*N
STKR(LL) = 1.0D0
47 CONTINUE
IF (TOP+1 .GE. BOT) CALL ERROR(18)
IF (ERR .GT. 0) RETURN
TOP = TOP+1
LSTK(TOP) = LE
MSTK(TOP) = N
NSTK(TOP) = N
GO TO 99
C
99 RETURN
END
SUBROUTINE MATFN5
C
C FILE HANDLING AND OTHER I/O
C
DOUBLE PRECISION STKR(5005),STKI(5005)
INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
INTEGER ALFA(52),ALFB(52),ALFL,CASE
INTEGER IDS(4,32),PSTK(32),RSTK(32),PSIZE,PT,PTZ
INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE
INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
COMMON /ALFS/ ALFA,ALFB,ALFL,CASE
COMMON /RECU/ IDS,PSTK,RSTK,PSIZE,PT,PTZ
COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE
COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
INTEGER EOL,CH,BLANK,FLAG,TOP2,PLUS,MINUS,QUOTE,SEMI,LRAT,MRAT
INTEGER ID(4)
DOUBLE PRECISION EPS,B,S,T,FLOP,WASUM
LOGICAL TEXT
DATA EOL/99/,BLANK/36/,PLUS/41/,MINUS/42/,QUOTE/49/,SEMI/39/
DATA LRAT/5/,MRAT/100/
C
IF (DDT .EQ. 1) WRITE(WTE,100) FIN
100 FORMAT(1X,'MATFN5',I4)
C FUNCTIONS/FIN
C EXEC SAVE LOAD PRIN DIAR DISP BASE LINE CHAR PLOT RAT DEBU
C 1 2 3 4 5 6 7 8 9 10 11 12
L = LSTK(TOP)
M = MSTK(TOP)
N = NSTK(TOP)
IF (FIN .GT. 5) GO TO 15
C
C CONVERT FILE NAME
MN = M*N
FLAG = 3
IF (SYM .EQ. SEMI) FLAG = 0
IF (RHS .LT. 2) GO TO 12
FLAG = IDINT(STKR(L))
TOP2 = TOP
TOP = TOP-1
L = LSTK(TOP)
MN = MSTK(TOP)*NSTK(TOP)
12 LUN = -1
IF (MN.EQ.1 .AND. STKR(L).LT.10.0D0) LUN = IDINT(STKR(L))
IF (LUN .GE. 0) GO TO 15
DO 14 J = 1, 32
LS = L+J-1
IF (J .LE. MN) CH = IDINT(STKR(LS))
IF (J .GT. MN) CH = BLANK
IF (CH.LT.0 .OR. CH.GE.ALFL) CALL ERROR(38)
IF (ERR .GT. 0) RETURN
IF (CASE .EQ. 0) BUF(J) = ALFA(CH+1)
IF (CASE .EQ. 1) BUF(J) = ALFB(CH+1)
14 CONTINUE
C
15 GO TO (20,30,35,25,27,60,65,70,50,80,40,95),FIN
C
C EXEC
20 IF (LUN .EQ. 0) GO TO 23
K = LPT(6)
LIN(K+1) = LPT(1)
LIN(K+2) = LPT(3)
LIN(K+3) = LPT(6)
LIN(K+4) = PTZ
LIN(K+5) = RIO
LIN(K+6) = LCT(4)
LPT(1) = K + 7
LCT(4) = FLAG
PTZ = PT - 4
IF (RIO .EQ. RTE) RIO = 12
RIO = RIO + 1
IF (LUN .GT. 0) RIO = LUN
IF (LUN .LT. 0) CALL FILES(RIO,BUF)
IF (FLAG .GE. 4) WRITE(WTE,22)
22 FORMAT(1X,'PAUSE MODE. ENTER BLANK LINES.')
SYM = EOL
MSTK(TOP) = 0
GO TO 99
C
C EXEC(0)
23 RIO = RTE
ERR = 99
GO TO 99
C
C PRINT
25 K = WTE
WTE = LUN
IF (LUN .LT. 0) WTE = 7
IF (LUN .LT. 0) CALL FILES(WTE,BUF)
L = LCT(2)
LCT(2) = 9999
IF (RHS .GT. 1) CALL PRINT(SYN,TOP2)
LCT(2) = L
WTE = K
MSTK(TOP) = 0
GO TO 99
C
C DIARY
27 WIO = LUN
IF (LUN .LT. 0) WIO = 8
IF (LUN .LT. 0) CALL FILES(WIO,BUF)
MSTK(TOP) = 0
GO TO 99
C
C SAVE
30 IF (LUN .LT. 0) LUNIT = 1
IF (LUN .LT. 0) CALL FILES(LUNIT,BUF)
IF (LUN .GT. 0) LUNIT = LUN
K = LSIZE-4
IF (K .LT. BOT) K = LSIZE
IF (RHS .EQ. 2) K = TOP2
IF (RHS .EQ. 2) CALL PUTID(IDSTK(1,K),SYN)
32 L = LSTK(K)
M = MSTK(K)
N = NSTK(K)
DO 34 I = 1, 4
J = IDSTK(I,K)+1
BUF(I) = ALFA(J)
34 CONTINUE
IMG = 0
IF (WASUM(M*N,STKI(L),STKI(L),1) .NE. 0.0D0) IMG = 1
IF(FE .EQ. 0)CALL SAVLOD(LUNIT,BUF,M,N,IMG,0,STKR(L),STKI(L))
K = K-1
IF (K .GE. BOT) GO TO 32
CALL FILES(-LUNIT,BUF)
MSTK(TOP) = 0
GO TO 99
C
C LOAD
35 IF (LUN .LT. 0) LUNIT = 2
IF (LUN .LT. 0) CALL FILES(LUNIT,BUF)
IF (LUN .GT. 0) LUNIT = LUN
36 JOB = LSTK(BOT) - L
IF(FE .EQ. 0)
+CALL SAVLOD(LUNIT,ID,MSTK(TOP),NSTK(TOP),IMG,JOB,STKR(L),STKI(L))
MN = MSTK(TOP)*NSTK(TOP)
IF (MN .EQ. 0) GO TO 39
IF (IMG .EQ. 0) CALL RSET(MN,0.0D0,STKI(L),1)
DO 38 I = 1, 4
J = 0
37 J = J+1
IF (ID(I).NE.ALFA(J) .AND. J.LE.BLANK) GO TO 37
ID(I) = J-1
38 CONTINUE
SYM = SEMI
RHS = 0
CALL STACKP(ID)
TOP = TOP + 1
GO TO 36
39 CALL FILES(-LUNIT,BUF)
MSTK(TOP) = 0
GO TO 99
C
C RAT
40 IF (RHS .EQ. 2) GO TO 44
MN = M*N
L2 = L
IF (LHS .EQ. 2) L2 = L + MN
LW = L2 + MN
ERR = LW + LRAT - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
IF (LHS .EQ. 2) TOP = TOP + 1
LSTK(TOP) = L2
MSTK(TOP) = M
NSTK(TOP) = N
CALL RSET(LHS*MN,0.0D0,STKI(L),1)
DO 42 I = 1, MN
CALL RAT(STKR(L),LRAT,MRAT,S,T,STKR(LW))
STKR(L) = S
STKR(L2) = T
IF (LHS .EQ. 1) STKR(L) = FLOP(S/T)
L = L + 1
L2 = L2 + 1
42 CONTINUE
GO TO 99
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.