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 44
Archive-name: applications/matlab/src.4
# 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-4
# This archive created: Wed Nov 2 16:20:28 1988
cat << \SHAR_EOF > src-4
IF (OP .GT. DOT) GO TO 70
C
C ADDITION
01 IF (M .LT. 0) GO TO 50
IF (M2 .LT. 0) GO TO 52
IF (M .NE. M2) CALL ERROR(8)
IF (ERR .GT. 0) RETURN
IF (N .NE. N2) CALL ERROR(8)
IF (ERR .GT. 0) RETURN
CALL WAXPY(M*N,1.0D0,0.0D0,STKR(L2),STKI(L2),1,
$ STKR(L),STKI(L),1)
GO TO 99
C
C SUBTRACTION
03 IF (M .LT. 0) GO TO 54
IF (M2 .LT. 0) GO TO 56
IF (M .NE. M2) CALL ERROR(9)
IF (ERR .GT. 0) RETURN
IF (N .NE. N2) CALL ERROR(9)
IF (ERR .GT. 0) RETURN
CALL WAXPY(M*N,-1.0D0,0.0D0,STKR(L2),STKI(L2),1,
$ STKR(L),STKI(L),1)
GO TO 99
C
C MULTIPLICATION
05 IF (M2*M2*N2 .EQ. 1) GO TO 10
IF (M*N .EQ. 1) GO TO 11
IF (M2*N2 .EQ. 1) GO TO 10
IF (N .NE. M2) CALL ERROR(10)
IF (ERR .GT. 0) RETURN
MN = M*N2
LL = L + MN
ERR = LL+M*N+M2*N2 - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
CALL WCOPY(M*N+M2*N2,STKR(L),STKI(L),-1,STKR(LL),STKI(LL),-1)
DO 08 J = 1, N2
DO 08 I = 1, M
K1 = L + MN + (I-1)
K2 = L2 + MN + (J-1)*M2
K = L + (I-1) + (J-1)*M
STKR(K) = WDOTUR(N,STKR(K1),STKI(K1),M,STKR(K2),STKI(K2),1)
STKI(K) = WDOTUI(N,STKR(K1),STKI(K1),M,STKR(K2),STKI(K2),1)
08 CONTINUE
NSTK(TOP) = N2
GO TO 99
C
C MULTIPLICATION BY SCALAR
10 SR = STKR(L2)
SI = STKI(L2)
L1 = L
GO TO 13
11 SR = STKR(L)
SI = STKI(L)
L1 = L+1
MSTK(TOP) = M2
NSTK(TOP) = N2
13 MN = MSTK(TOP)*NSTK(TOP)
CALL WSCAL(MN,SR,SI,STKR(L1),STKI(L1),1)
IF (L1.NE.L)
$ CALL WCOPY(MN,STKR(L1),STKI(L1),1,STKR(L),STKI(L),1)
GO TO 99
C
C RIGHT DIVISION
20 IF (M2*N2 .EQ. 1) GO TO 21
IF (M2 .EQ. N2) FUN = 1
IF (M2 .NE. N2) FUN = 4
FIN = -1
RHS = 2
GO TO 99
21 SR = STKR(L2)
SI = STKI(L2)
MN = M*N
DO 22 I = 1, MN
LL = L+I-1
CALL WDIV(STKR(LL),STKI(LL),SR,SI,STKR(LL),STKI(LL))
IF (ERR .GT. 0) RETURN
22 CONTINUE
GO TO 99
C
C LEFT DIVISION
25 IF (M*N .EQ. 1) GO TO 26
IF (M .EQ. N) FUN = 1
IF (M .NE. N) FUN = 4
FIN = -2
RHS = 2
GO TO 99
26 SR = STKR(L)
SI = STKI(L)
MSTK(TOP) = M2
NSTK(TOP) = N2
MN = M2*N2
DO 27 I = 1, MN
LL = L+I-1
CALL WDIV(STKR(LL+1),STKI(LL+1),SR,SI,STKR(LL),STKI(LL))
IF (ERR .GT. 0) RETURN
27 CONTINUE
GO TO 99
C
C POWER
30 IF (M2*N2 .NE. 1) CALL ERROR(30)
IF (ERR .GT. 0) RETURN
IF (M .NE. N) CALL ERROR(20)
IF (ERR .GT. 0) RETURN
NEXP = IDINT(STKR(L2))
IF (STKR(L2) .NE. DFLOAT(NEXP)) GO TO 39
IF (STKI(L2) .NE. 0.0D0) GO TO 39
IF (NEXP .LT. 2) GO TO 39
MN = M*N
ERR = L2+MN+N - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
CALL WCOPY(MN,STKR(L),STKI(L),1,STKR(L2),STKI(L2),1)
L3 = L2+MN
DO 36 KEXP = 2, NEXP
DO 35 J = 1, N
LS = L+(J-1)*N
CALL WCOPY(N,STKR(LS),STKI(LS),1,STKR(L3),STKI(L3),1)
DO 34 I = 1, N
LS = L2+I-1
LL = L+I-1+(J-1)*N
STKR(LL) = WDOTUR(N,STKR(LS),STKI(LS),N,STKR(L3),STKI(L3),1)
STKI(LL) = WDOTUI(N,STKR(LS),STKI(LS),N,STKR(L3),STKI(L3),1)
34 CONTINUE
35 CONTINUE
36 CONTINUE
GO TO 99
C
C NONINTEGER OR NONPOSITIVE POWER, USE EIGENVECTORS
39 FUN = 2
FIN = 0
GO TO 99
C
C ADD OR SUBTRACT SCALAR
50 IF (M2 .NE. N2) CALL ERROR(8)
IF (ERR .GT. 0) RETURN
M = M2
N = N2
MSTK(TOP) = M
NSTK(TOP) = N
SR = STKR(L)
SI = STKI(L)
CALL WCOPY(M*N,STKR(L+1),STKI(L+1),1,STKR(L),STKI(L),1)
GO TO 58
52 IF (M .NE. N) CALL ERROR(8)
IF (ERR .GT. 0) RETURN
SR = STKR(L2)
SI = STKI(L2)
GO TO 58
54 IF (M2 .NE. N2) CALL ERROR(9)
IF (ERR .GT. 0) RETURN
M = M2
N = N2
MSTK(TOP) = M
NSTK(TOP) = N
SR = STKR(L)
SI = STKI(L)
CALL WCOPY(M*N,STKR(L+1),STKI(L+1),1,STKR(L),STKI(L),1)
CALL WRSCAL(M*N,-1.0D0,STKR(L),STKI(L),1)
GO TO 58
56 IF (M .NE. N) CALL ERROR(9)
IF (ERR .GT. 0) RETURN
SR = -STKR(L2)
SI = -STKI(L2)
GO TO 58
58 DO 59 I = 1, N
LL = L + (I-1)*(N+1)
STKR(LL) = FLOP(STKR(LL)+SR)
STKI(LL) = FLOP(STKI(LL)+SI)
59 CONTINUE
GO TO 99
C
C COLON
60 E2 = STKR(L2)
ST = 1.0D0
N = 0
IF (RHS .LT. 3) GO TO 61
ST = STKR(L)
TOP = TOP-1
L = LSTK(TOP)
IF (ST .EQ. 0.0D0) GO TO 63
61 E1 = STKR(L)
C CHECK FOR CLAUSE
IF (RSTK(PT) .EQ. 3) GO TO 64
ERR = L + MAX0(3,IDINT((E2-E1)/ST)) - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
62 IF (ST .GT. 0.0D0 .AND. STKR(L) .GT. E2) GO TO 63
IF (ST .LT. 0.0D0 .AND. STKR(L) .LT. E2) GO TO 63
N = N+1
L = L+1
STKR(L) = E1 + DFLOAT(N)*ST
STKI(L) = 0.0D0
GO TO 62
63 NSTK(TOP) = N
MSTK(TOP) = 1
IF (N .EQ. 0) MSTK(TOP) = 0
GO TO 99
C
C FOR CLAUSE
64 STKR(L) = E1
STKR(L+1) = ST
STKR(L+2) = E2
MSTK(TOP) = -3
NSTK(TOP) = -1
GO TO 99
C
C ELEMENTWISE OPERATIONS
70 OP = OP - DOT
IF (M.NE.M2 .OR. N.NE.N2) CALL ERROR(10)
IF (ERR .GT. 0) RETURN
MN = M*N
DO 72 I = 1, MN
J = L+I-1
K = L2+I-1
IF (OP .EQ. STAR)
$ CALL WMUL(STKR(J),STKI(J),STKR(K),STKI(K),STKR(J),STKI(J))
IF (OP .EQ. SLASH)
$ CALL WDIV(STKR(J),STKI(J),STKR(K),STKI(K),STKR(J),STKI(J))
IF (OP .EQ. BSLASH)
$ CALL WDIV(STKR(K),STKI(K),STKR(J),STKI(J),STKR(J),STKI(J))
IF (ERR .GT. 0) RETURN
72 CONTINUE
GO TO 99
C
C KRONECKER
80 FIN = OP - 2*DOT - STAR + 11
FUN = 6
TOP = TOP + 1
RHS = 2
GO TO 99
C
99 RETURN
END
SUBROUTINE STACKG(ID)
INTEGER ID(4)
C
C GET VARIABLES FROM STORAGE
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 EQID
IF (DDT .EQ. 1) WRITE(WTE,100) ID
100 FORMAT(1X,'STACKG',4I4)
CALL PUTID(IDSTK(1,BOT-1), ID)
K = LSIZE+1
10 K = K-1
IF (.NOT.EQID(IDSTK(1,K), ID)) GO TO 10
IF (K .GE. LSIZE-1 .AND. RHS .GT. 0) GO TO 98
IF (K .EQ. BOT-1) GO TO 98
LK = LSTK(K)
IF (RHS .EQ. 1) GO TO 40
IF (RHS .EQ. 2) GO TO 60
IF (RHS .GT. 2) CALL ERROR(21)
IF (ERR .GT. 0) RETURN
L = 1
IF (TOP .GT. 0) L = LSTK(TOP) + MSTK(TOP)*NSTK(TOP)
IF (TOP+1 .GE. BOT) CALL ERROR(18)
IF (ERR .GT. 0) RETURN
TOP = TOP+1
C
C LOAD VARIABLE TO TOP OF STACK
LSTK(TOP) = L
MSTK(TOP) = MSTK(K)
NSTK(TOP) = NSTK(K)
MN = MSTK(K)*NSTK(K)
ERR = L+MN - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
C IF RAND, MATFN6 GENERATES RANDOM NUMBER
IF (K .EQ. LSIZE) GO TO 97
CALL WCOPY(MN,STKR(LK),STKI(LK),1,STKR(L),STKI(L),1)
GO TO 99
C
C VECT(ARG)
40 IF (MSTK(TOP) .EQ. 0) GO TO 99
L = LSTK(TOP)
MN = MSTK(TOP)*NSTK(TOP)
MNK = MSTK(K)*NSTK(K)
IF (MSTK(TOP) .LT. 0) MN = MNK
DO 50 I = 1, MN
LL = L+I-1
LS = LK+I-1
IF (MSTK(TOP) .GT. 0) LS = LK + IDINT(STKR(LL)) - 1
IF (LS .LT. LK .OR. LS .GE. LK+MNK) CALL ERROR(21)
IF (ERR .GT. 0) RETURN
STKR(LL) = STKR(LS)
STKI(LL) = STKI(LS)
50 CONTINUE
MSTK(TOP) = 1
NSTK(TOP) = 1
IF (MSTK(K) .GT. 1) MSTK(TOP) = MN
IF (MSTK(K) .EQ. 1) NSTK(TOP) = MN
GO TO 99
C
C MATRIX(ARG,ARG)
60 TOP = TOP-1
L = LSTK(TOP)
IF (MSTK(TOP+1) .EQ. 0) MSTK(TOP) = 0
IF (MSTK(TOP) .EQ. 0) GO TO 99
L2 = LSTK(TOP+1)
M = MSTK(TOP)*NSTK(TOP)
IF (MSTK(TOP) .LT. 0) M = MSTK(K)
N = MSTK(TOP+1)*NSTK(TOP+1)
IF (MSTK(TOP+1) .LT. 0) N = NSTK(K)
L3 = L2 + N
MK = MSTK(K)
MNK = MSTK(K)*NSTK(K)
DO 70 J = 1, N
DO 70 I = 1, M
LI = L+I-1
IF (MSTK(TOP) .GT. 0) LI = L + IDINT(STKR(LI)) - 1
LJ = L2+J-1
IF (MSTK(TOP+1) .GT. 0) LJ = L2 + IDINT(STKR(LJ)) - 1
LS = LK + LI-L + (LJ-L2)*MK
IF (LS.LT.LK .OR. LS.GE.LK+MNK) CALL ERROR(21)
IF (ERR .GT. 0) RETURN
LL = L3 + I-1 + (J-1)*M
STKR(LL) = STKR(LS)
STKI(LL) = STKI(LS)
70 CONTINUE
MN = M*N
CALL WCOPY(MN,STKR(L3),STKI(L3),1,STKR(L),STKI(L),1)
MSTK(TOP) = M
NSTK(TOP) = N
GO TO 99
97 FIN = 7
FUN = 6
RETURN
98 FIN = 0
RETURN
99 FIN = -1
FUN = 0
RETURN
END
SUBROUTINE STACKP(ID)
INTEGER ID(4)
C
C PUT VARIABLES INTO STORAGE
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 EQID
INTEGER SEMI
DATA SEMI/39/
IF (DDT .EQ. 1) WRITE(WTE,100) ID
100 FORMAT(1X,'STACKP',4I4)
IF (TOP .LE. 0) CALL ERROR(1)
IF (ERR .GT. 0) RETURN
CALL FUNS(ID)
IF (FIN .NE. 0) CALL ERROR(25)
IF (ERR .GT. 0) RETURN
M = MSTK(TOP)
N = NSTK(TOP)
IF (M .GT. 0) L = LSTK(TOP)
IF (M .LT. 0) CALL ERROR(14)
IF (ERR .GT. 0) RETURN
IF (M .EQ. 0 .AND. N .NE. 0) GO TO 99
MN = M*N
LK = 0
MK = 1
NK = 0
LT = 0
MT = 0
NT = 0
C
C DOES VARIABLE ALREADY EXIST
CALL PUTID(IDSTK(1,BOT-1),ID)
K = LSIZE+1
05 K = K-1
IF (.NOT.EQID(IDSTK(1,K),ID)) GO TO 05
IF (K .EQ. BOT-1) GO TO 30
LK = LSTK(K)
MK = MSTK(K)
NK = NSTK(K)
MNK = MK*NK
IF (RHS .EQ. 0) GO TO 20
IF (RHS .GT. 2) CALL ERROR(15)
IF (ERR .GT. 0) RETURN
MT = MK
NT = NK
LT = L + MN
ERR = LT + MNK - LSTK(BOT)
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
CALL WCOPY(MNK,STKR(LK),STKI(LK),1,STKR(LT),STKI(LT),1)
C
C DOES IT FIT
20 IF (RHS.EQ.0 .AND. MN.EQ.MNK) GO TO 40
IF (K .GE. LSIZE-3) CALL ERROR(13)
IF (ERR .GT. 0) RETURN
C
C SHIFT STORAGE
IF (K .EQ. BOT) GO TO 25
LS = LSTK(BOT)
LL = LS + MNK
CALL WCOPY(LK-LS,STKR(LS),STKI(LS),-1,STKR(LL),STKI(LL),-1)
KM1 = K-1
DO 24 IB = BOT, KM1
I = BOT+KM1-IB
CALL PUTID(IDSTK(1,I+1),IDSTK(1,I))
MSTK(I+1) = MSTK(I)
NSTK(I+1) = NSTK(I)
LSTK(I+1) = LSTK(I)+MNK
24 CONTINUE
C
C DESTROY OLD VARIABLE
25 BOT = BOT+1
C
C CREATE NEW VARIABLE
30 IF (MN .EQ. 0) GO TO 99
IF (BOT-2 .LE. TOP) CALL ERROR(18)
IF (ERR .GT. 0) RETURN
K = BOT-1
CALL PUTID(IDSTK(1,K), ID)
IF (RHS .EQ. 1) GO TO 50
IF (RHS .EQ. 2) GO TO 55
C
C STORE
40 IF (K .LT. LSIZE) LSTK(K) = LSTK(K+1) - MN
MSTK(K) = M
NSTK(K) = N
LK = LSTK(K)
CALL WCOPY(MN,STKR(L),STKI(L),-1,STKR(LK),STKI(LK),-1)
GO TO 90
C
C VECT(ARG)
50 IF (MSTK(TOP-1) .LT. 0) GO TO 59
MN1 = 1
MN2 = 1
L1 = 0
L2 = 0
IF (N.NE.1 .OR. NK.NE.1) GO TO 52
L1 = LSTK(TOP-1)
M1 = MSTK(TOP-1)
MN1 = M1*NSTK(TOP-1)
M2 = -1
GO TO 60
52 IF (M.NE.1 .OR. MK.NE.1) CALL ERROR(15)
IF (ERR .GT. 0) RETURN
L2 = LSTK(TOP-1)
M2 = MSTK(TOP-1)
MN2 = M2*NSTK(TOP-1)
M1 = -1
GO TO 60
C
C MATRIX(ARG,ARG)
55 IF (MSTK(TOP-1).LT.0 .AND. MSTK(TOP-2).LT.0) GO TO 59
L2 = LSTK(TOP-1)
M2 = MSTK(TOP-1)
MN2 = M2*NSTK(TOP-1)
IF (M2 .LT. 0) MN2 = N
L1 = LSTK(TOP-2)
M1 = MSTK(TOP-2)
MN1 = M1*NSTK(TOP-2)
IF (M1 .LT. 0) MN1 = M
GO TO 60
C
59 IF (MN .NE. MNK) CALL ERROR(15)
IF (ERR .GT. 0) RETURN
LK = LSTK(K)
CALL WCOPY(MN,STKR(L),STKI(L),-1,STKR(LK),STKI(LK),-1)
GO TO 90
C
60 IF (MN1.NE.M .OR. MN2.NE.N) CALL ERROR(15)
IF (ERR .GT. 0) RETURN
LL = 1
IF (M1 .LT. 0) GO TO 62
DO 61 I = 1, MN1
LS = L1+I-1
MK = MAX0(MK,IDINT(STKR(LS)))
LL = MIN0(LL,IDINT(STKR(LS)))
61 CONTINUE
62 MK = MAX0(MK,M)
IF (M2 .LT. 0) GO TO 64
DO 63 I = 1, MN2
LS = L2+I-1
NK = MAX0(NK,IDINT(STKR(LS)))
LL = MIN0(LL,IDINT(STKR(LS)))
63 CONTINUE
64 NK = MAX0(NK,N)
IF (LL .LT. 1) CALL ERROR(21)
IF (ERR .GT. 0) RETURN
MNK = MK*NK
LK = LSTK(K+1) - MNK
ERR = LT + MT*NT - LK
IF (ERR .GT. 0) CALL ERROR(17)
IF (ERR .GT. 0) RETURN
LSTK(K) = LK
MSTK(K) = MK
NSTK(K) = NK
CALL WSET(MNK,0.0D0,0.0D0,STKR(LK),STKI(LK),1)
IF (NT .LT. 1) GO TO 67
DO 66 J = 1, NT
LS = LT+(J-1)*MT
LL = LK+(J-1)*MK
CALL WCOPY(MT,STKR(LS),STKI(LS),-1,STKR(LL),STKI(LL),-1)
66 CONTINUE
67 DO 68 J = 1, N
DO 68 I = 1, M
LI = L1+I-1
IF (M1 .GT. 0) LI = L1 + IDINT(STKR(LI)) - 1
LJ = L2+J-1
IF (M2 .GT. 0) LJ = L2 + IDINT(STKR(LJ)) - 1
LL = LK+LI-L1+(LJ-L2)*MK
LS = L+I-1+(J-1)*M
STKR(LL) = STKR(LS)
STKI(LL) = STKI(LS)
68 CONTINUE
GO TO 90
C
C PRINT IF DESIRED AND POP STACK
90 IF (SYM.NE.SEMI .AND. LCT(3).EQ.0) CALL PRINT(ID,K)
IF (SYM.EQ.SEMI .AND. LCT(3).EQ.1) CALL PRINT(ID,K)
IF (K .EQ. BOT-1) BOT = BOT-1
99 IF (M .NE. 0) TOP = TOP - 1 - RHS
IF (M .EQ. 0) TOP = TOP - 1
RETURN
END
SUBROUTINE TERM
DOUBLE PRECISION STKR(5005),STKI(5005)
INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
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 /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 R,OP,BSLASH,STAR,SLASH,DOT
DATA BSLASH/45/,STAR/43/,SLASH/44/,DOT/47/
IF (DDT .EQ. 1) WRITE(WTE,100) PT,RSTK(PT)
100 FORMAT(1X,'TERM ',2I4)
R = RSTK(PT)
GO TO (99,99,99,99,99,01,01,05,25,99,99,99,99,99,35,99,99,99,99),R
01 PT = PT+1
RSTK(PT) = 8
C *CALL* FACTOR
RETURN
05 PT = PT-1
10 OP = 0
IF (SYM .EQ. DOT) OP = DOT
IF (SYM .EQ. DOT) CALL GETSYM
IF (SYM.EQ.STAR .OR. SYM.EQ.SLASH .OR. SYM.EQ.BSLASH) GO TO 20
RETURN
20 OP = OP + SYM
CALL GETSYM
IF (SYM .EQ. DOT) OP = OP + SYM
IF (SYM .EQ. DOT) CALL GETSYM
PT = PT+1
PSTK(PT) = OP
RSTK(PT) = 9
C *CALL* FACTOR
RETURN
25 OP = PSTK(PT)
PT = PT-1
CALL STACK2(OP)
IF (ERR .GT. 0) RETURN
C SOME BINARY OPS DONE IN MATFNS
IF (FUN .EQ. 0) GO TO 10
PT = PT+1
RSTK(PT) = 15
C *CALL* MATFN
RETURN
35 PT = PT-1
GO TO 10
99 CALL ERROR(22)
IF (ERR .GT. 0) RETURN
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) = DFLOAT(A3(I))
10 CONTINUE
M = 3
N = 3
RETURN
END
SUBROUTINE XCHAR(BUF,K)
INTEGER BUF(1),K
C
C SYSTEM DEPENDENT ROUTINE TO HANDLE SPECIAL CHARACTERS
C
C
INTEGER BACK,MASK
DATA BACK/Z'20202008'/,MASK/Z'000000FF'/
C
IF (BUF(1) .EQ. BACK) K = -1
L = BUF(1) .AND. MASK
IF (K .NE. -1) WRITE(6,10) BUF(1),L
10 FORMAT(1X,1H',A1,4H' = ,Z2,' hex is not a MATLAB character.')
RETURN
END
SUBROUTINE WGECO(AR,AI,LDA,N,IPVT,RCOND,ZR,ZI)
INTEGER LDA,N,IPVT(1)
DOUBLE PRECISION AR(LDA,1),AI(LDA,1),ZR(1),ZI(1)
DOUBLE PRECISION RCOND
C
C WGECO FACTORS A DOUBLE-COMPLEX MATRIX BY GAUSSIAN ELIMINATION
C AND ESTIMATES THE CONDITION OF THE MATRIX.
C
C IF RCOND IS NOT NEEDED, WGEFA IS SLIGHTLY FASTER.
C TO SOLVE A*X = B , FOLLOW WGECO BY WGESL.
C TO COMPUTE INVERSE(A)*C , FOLLOW WGECO BY WGESL.
C TO COMPUTE DETERMINANT(A) , FOLLOW WGECO BY WGEDI.
C TO COMPUTE INVERSE(A) , FOLLOW WGECO BY WGEDI.
C
C ON ENTRY
C
C A DOUBLE-COMPLEX(LDA, N)
C THE MATRIX TO BE FACTORED.
C
C LDA INTEGER
C THE LEADING DIMENSION OF THE ARRAY A .
C
C N INTEGER
C THE ORDER OF THE MATRIX A .
C
C ON RETURN
C
C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS
C WHICH WERE USED TO OBTAIN IT.
C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE
C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER
C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR.
C
C IPVT INTEGER(N)
C AN INTEGER VECTOR OF PIVOT INDICES.
C
C RCOND DOUBLE PRECISION
C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A .
C FOR THE SYSTEM A*X = B , RELATIVE PERTURBATIONS
C IN A AND B OF SIZE EPSILON MAY CAUSE
C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND .
C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION
C 1.0 + RCOND .EQ. 1.0
C IS TRUE, THEN A MAY BE SINGULAR TO WORKING
C PRECISION. IN PARTICULAR, RCOND IS ZERO IF
C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE
C UNDERFLOWS.
C
C Z DOUBLE-COMPLEX(N)
C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT.
C IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS
C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT
C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) .
C
C LINPACK. THIS VERSION DATED 07/01/79 .
C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C SUBROUTINES AND FUNCTIONS
C
C LINPACK WGEFA
C BLAS WAXPY,WDOTC,WASUM
C FORTRAN DABS,DMAX1
C
C INTERNAL VARIABLES
C
DOUBLE PRECISION WDOTCR,WDOTCI,EKR,EKI,TR,TI,WKR,WKI,WKMR,WKMI
DOUBLE PRECISION ANORM,S,WASUM,SM,YNORM,FLOP
INTEGER INFO,J,K,KB,KP1,L
C
DOUBLE PRECISION ZDUMR,ZDUMI
DOUBLE PRECISION CABS1
CABS1(ZDUMR,ZDUMI) = DABS(ZDUMR) + DABS(ZDUMI)
C
C COMPUTE 1-NORM OF A
C
ANORM = 0.0D0
DO 10 J = 1, N
ANORM = DMAX1(ANORM,WASUM(N,AR(1,J),AI(1,J),1))
10 CONTINUE
C
C FACTOR
C
CALL WGEFA(AR,AI,LDA,N,IPVT,INFO)
C
C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) .
C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND CTRANS(A)*Y = E .
C CTRANS(A) IS THE CONJUGATE TRANSPOSE OF A .
C THE COMPONENTS OF E ARE CHOSEN TO CAUSE MAXIMUM LOCAL
C GROWTH IN THE ELEMENTS OF W WHERE CTRANS(U)*W = E .
C THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW.
C
C SOLVE CTRANS(U)*W = E
C
EKR = 1.0D0
EKI = 0.0D0
DO 20 J = 1, N
ZR(J) = 0.0D0
ZI(J) = 0.0D0
20 CONTINUE
DO 110 K = 1, N
CALL WSIGN(EKR,EKI,-ZR(K),-ZI(K),EKR,EKI)
IF (CABS1(EKR-ZR(K),EKI-ZI(K))
* .LE. CABS1(AR(K,K),AI(K,K))) GO TO 40
S = CABS1(AR(K,K),AI(K,K))
* /CABS1(EKR-ZR(K),EKI-ZI(K))
CALL WRSCAL(N,S,ZR,ZI,1)
EKR = S*EKR
EKI = S*EKI
40 CONTINUE
WKR = EKR - ZR(K)
WKI = EKI - ZI(K)
WKMR = -EKR - ZR(K)
WKMI = -EKI - ZI(K)
S = CABS1(WKR,WKI)
SM = CABS1(WKMR,WKMI)
IF (CABS1(AR(K,K),AI(K,K)) .EQ. 0.0D0) GO TO 50
CALL WDIV(WKR,WKI,AR(K,K),-AI(K,K),WKR,WKI)
CALL WDIV(WKMR,WKMI,AR(K,K),-AI(K,K),WKMR,WKMI)
GO TO 60
50 CONTINUE
WKR = 1.0D0
WKI = 0.0D0
WKMR = 1.0D0
WKMI = 0.0D0
60 CONTINUE
KP1 = K + 1
IF (KP1 .GT. N) GO TO 100
DO 70 J = KP1, N
CALL WMUL(WKMR,WKMI,AR(K,J),-AI(K,J),TR,TI)
SM = FLOP(SM + CABS1(ZR(J)+TR,ZI(J)+TI))
CALL WAXPY(1,WKR,WKI,AR(K,J),-AI(K,J),1,
$ ZR(J),ZI(J),1)
S = FLOP(S + CABS1(ZR(J),ZI(J)))
70 CONTINUE
IF (S .GE. SM) GO TO 90
TR = WKMR - WKR
TI = WKMI - WKI
WKR = WKMR
WKI = WKMI
DO 80 J = KP1, N
CALL WAXPY(1,TR,TI,AR(K,J),-AI(K,J),1,
$ ZR(J),ZI(J),1)
80 CONTINUE
90 CONTINUE
100 CONTINUE
ZR(K) = WKR
ZI(K) = WKI
110 CONTINUE
S = 1.0D0/WASUM(N,ZR,ZI,1)
CALL WRSCAL(N,S,ZR,ZI,1)
C
C SOLVE CTRANS(L)*Y = W
C
DO 140 KB = 1, N
K = N + 1 - KB
IF (K .GE. N) GO TO 120
ZR(K) = ZR(K)
* + WDOTCR(N-K,AR(K+1,K),AI(K+1,K),1,ZR(K+1),ZI(K+1),1)
ZI(K) = ZI(K)
* + WDOTCI(N-K,AR(K+1,K),AI(K+1,K),1,ZR(K+1),ZI(K+1),1)
120 CONTINUE
IF (CABS1(ZR(K),ZI(K)) .LE. 1.0D0) GO TO 130
S = 1.0D0/CABS1(ZR(K),ZI(K))
CALL WRSCAL(N,S,ZR,ZI,1)
130 CONTINUE
L = IPVT(K)
TR = ZR(L)
TI = ZI(L)
ZR(L) = ZR(K)
ZI(L) = ZI(K)
ZR(K) = TR
ZI(K) = TI
140 CONTINUE
S = 1.0D0/WASUM(N,ZR,ZI,1)
CALL WRSCAL(N,S,ZR,ZI,1)
C
YNORM = 1.0D0
C
C SOLVE L*V = Y
C
DO 160 K = 1, N
L = IPVT(K)
TR = ZR(L)
TI = ZI(L)
ZR(L) = ZR(K)
ZI(L) = ZI(K)
ZR(K) = TR
ZI(K) = TI
IF (K .LT. N)
* CALL WAXPY(N-K,TR,TI,AR(K+1,K),AI(K+1,K),1,ZR(K+1),ZI(K+1),
* 1)
IF (CABS1(ZR(K),ZI(K)) .LE. 1.0D0) GO TO 150
S = 1.0D0/CABS1(ZR(K),ZI(K))
CALL WRSCAL(N,S,ZR,ZI,1)
YNORM = S*YNORM
150 CONTINUE
160 CONTINUE
S = 1.0D0/WASUM(N,ZR,ZI,1)
CALL WRSCAL(N,S,ZR,ZI,1)
YNORM = S*YNORM
C
C SOLVE U*Z = V
C
DO 200 KB = 1, N
K = N + 1 - KB
IF (CABS1(ZR(K),ZI(K))
* .LE. CABS1(AR(K,K),AI(K,K))) GO TO 170
S = CABS1(AR(K,K),AI(K,K))
* /CABS1(ZR(K),ZI(K))
CALL WRSCAL(N,S,ZR,ZI,1)
YNORM = S*YNORM
170 CONTINUE
IF (CABS1(AR(K,K),AI(K,K)) .EQ. 0.0D0) GO TO 180
CALL WDIV(ZR(K),ZI(K),AR(K,K),AI(K,K),ZR(K),ZI(K))
180 CONTINUE
IF (CABS1(AR(K,K),AI(K,K)) .NE. 0.0D0) GO TO 190
ZR(K) = 1.0D0
ZI(K) = 0.0D0
190 CONTINUE
TR = -ZR(K)
TI = -ZI(K)
CALL WAXPY(K-1,TR,TI,AR(1,K),AI(1,K),1,ZR(1),ZI(1),1)
200 CONTINUE
C MAKE ZNORM = 1.0
S = 1.0D0/WASUM(N,ZR,ZI,1)
CALL WRSCAL(N,S,ZR,ZI,1)
YNORM = S*YNORM
C
IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM
IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0
RETURN
END
SUBROUTINE WGEFA(AR,AI,LDA,N,IPVT,INFO)
INTEGER LDA,N,IPVT(1),INFO
DOUBLE PRECISION AR(LDA,1),AI(LDA,1)
C
C WGEFA FACTORS A DOUBLE-COMPLEX MATRIX BY GAUSSIAN ELIMINATION.
C
C WGEFA IS USUALLY CALLED BY WGECO, BUT IT CAN BE CALLED
C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED.
C (TIME FOR WGECO) = (1 + 9/N)*(TIME FOR WGEFA) .
C
C ON ENTRY
C
C A DOUBLE-COMPLEX(LDA, N)
C THE MATRIX TO BE FACTORED.
C
C LDA INTEGER
C THE LEADING DIMENSION OF THE ARRAY A .
C
C N INTEGER
C THE ORDER OF THE MATRIX A .
C
C ON RETURN
C
C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS
C WHICH WERE USED TO OBTAIN IT.
C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE
C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER
C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR.
C
C IPVT INTEGER(N)
C AN INTEGER VECTOR OF PIVOT INDICES.
C
C INFO INTEGER
C = 0 NORMAL VALUE.
C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR
C CONDITION FOR THIS SUBROUTINE, BUT IT DOES
C INDICATE THAT WGESL OR WGEDI WILL DIVIDE BY ZERO
C IF CALLED. USE RCOND IN WGECO FOR A RELIABLE
C INDICATION OF SINGULARITY.
C
C LINPACK. THIS VERSION DATED 07/01/79 .
C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C SUBROUTINES AND FUNCTIONS
C
C BLAS WAXPY,WSCAL,IWAMAX
C FORTRAN DABS
C
C INTERNAL VARIABLES
C
DOUBLE PRECISION TR,TI
INTEGER IWAMAX,J,K,KP1,L,NM1
C
DOUBLE PRECISION ZDUMR,ZDUMI
DOUBLE PRECISION CABS1
CABS1(ZDUMR,ZDUMI) = DABS(ZDUMR) + DABS(ZDUMI)
C
C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
C
INFO = 0
NM1 = N - 1
IF (NM1 .LT. 1) GO TO 70
DO 60 K = 1, NM1
KP1 = K + 1
C
C FIND L = PIVOT INDEX
C
L = IWAMAX(N-K+1,AR(K,K),AI(K,K),1) + K - 1
IPVT(K) = L
C
C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
C
IF (CABS1(AR(L,K),AI(L,K)) .EQ. 0.0D0) GO TO 40
C
C INTERCHANGE IF NECESSARY
C
IF (L .EQ. K) GO TO 10
TR = AR(L,K)
TI = AI(L,K)
AR(L,K) = AR(K,K)
AI(L,K) = AI(K,K)
AR(K,K) = TR
AI(K,K) = TI
10 CONTINUE
C
C COMPUTE MULTIPLIERS
C
CALL WDIV(-1.0D0,0.0D0,AR(K,K),AI(K,K),TR,TI)
CALL WSCAL(N-K,TR,TI,AR(K+1,K),AI(K+1,K),1)
C
C ROW ELIMINATION WITH COLUMN INDEXING
C
DO 30 J = KP1, N
TR = AR(L,J)
TI = AI(L,J)
IF (L .EQ. K) GO TO 20
AR(L,J) = AR(K,J)
AI(L,J) = AI(K,J)
AR(K,J) = TR
AI(K,J) = TI
20 CONTINUE
CALL WAXPY(N-K,TR,TI,AR(K+1,K),AI(K+1,K),1,AR(K+1,J),
* AI(K+1,J),1)
30 CONTINUE
GO TO 50
40 CONTINUE
INFO = K
50 CONTINUE
60 CONTINUE
70 CONTINUE
IPVT(N) = N
IF (CABS1(AR(N,N),AI(N,N)) .EQ. 0.0D0) INFO = N
RETURN
END
SUBROUTINE WGESL(AR,AI,LDA,N,IPVT,BR,BI,JOB)
INTEGER LDA,N,IPVT(1),JOB
DOUBLE PRECISION AR(LDA,1),AI(LDA,1),BR(1),BI(1)
C
C WGESL SOLVES THE DOUBLE-COMPLEX SYSTEM
C A * X = B OR CTRANS(A) * X = B
C USING THE FACTORS COMPUTED BY WGECO OR WGEFA.
C
C ON ENTRY
C
C A DOUBLE-COMPLEX(LDA, N)
C THE OUTPUT FROM WGECO OR WGEFA.
C
C LDA INTEGER
C THE LEADING DIMENSION OF THE ARRAY A .
C
C N INTEGER
C THE ORDER OF THE MATRIX A .
C
C IPVT INTEGER(N)
C THE PIVOT VECTOR FROM WGECO OR WGEFA.
C
C B DOUBLE-COMPLEX(N)
C THE RIGHT HAND SIDE VECTOR.
C
C JOB INTEGER
C = 0 TO SOLVE A*X = B ,
C = NONZERO TO SOLVE CTRANS(A)*X = B WHERE
C CTRANS(A) IS THE CONJUGATE TRANSPOSE.
C
C ON RETURN
C
C B THE SOLUTION VECTOR X .
C
C ERROR CONDITION
C
C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A
C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY
C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER
C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE
C CALLED CORRECTLY AND IF WGECO HAS SET RCOND .GT. 0.0
C OR WGEFA HAS SET INFO .EQ. 0 .
C
C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX
C WITH P COLUMNS
C CALL WGECO(A,LDA,N,IPVT,RCOND,Z)
C IF (RCOND IS TOO SMALL) GO TO ...
C DO 10 J = 1, P
C CALL WGESL(A,LDA,N,IPVT,C(1,J),0)
C 10 CONTINUE
C
C LINPACK. THIS VERSION DATED 07/01/79 .
C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C SUBROUTINES AND FUNCTIONS
C
C BLAS WAXPY,WDOTC
C
C INTERNAL VARIABLES
C
DOUBLE PRECISION WDOTCR,WDOTCI,TR,TI
INTEGER K,KB,L,NM1
C
NM1 = N - 1
IF (JOB .NE. 0) GO TO 50
C
C JOB = 0 , SOLVE A * X = B
C FIRST SOLVE L*Y = B
C
IF (NM1 .LT. 1) GO TO 30
DO 20 K = 1, NM1
L = IPVT(K)
TR = BR(L)
TI = BI(L)
IF (L .EQ. K) GO TO 10
BR(L) = BR(K)
BI(L) = BI(K)
BR(K) = TR
BI(K) = TI
10 CONTINUE
CALL WAXPY(N-K,TR,TI,AR(K+1,K),AI(K+1,K),1,BR(K+1),BI(K+1),
* 1)
20 CONTINUE
30 CONTINUE
C
C NOW SOLVE U*X = Y
C
DO 40 KB = 1, N
K = N + 1 - KB
CALL WDIV(BR(K),BI(K),AR(K,K),AI(K,K),BR(K),BI(K))
TR = -BR(K)
TI = -BI(K)
CALL WAXPY(K-1,TR,TI,AR(1,K),AI(1,K),1,BR(1),BI(1),1)
40 CONTINUE
GO TO 100
50 CONTINUE
C
C JOB = NONZERO, SOLVE CTRANS(A) * X = B
C FIRST SOLVE CTRANS(U)*Y = B
C
DO 60 K = 1, N
TR = BR(K) - WDOTCR(K-1,AR(1,K),AI(1,K),1,BR(1),BI(1),1)
TI = BI(K) - WDOTCI(K-1,AR(1,K),AI(1,K),1,BR(1),BI(1),1)
CALL WDIV(TR,TI,AR(K,K),-AI(K,K),BR(K),BI(K))
60 CONTINUE
C
C NOW SOLVE CTRANS(L)*X = Y
C
IF (NM1 .LT. 1) GO TO 90
DO 80 KB = 1, NM1
K = N - KB
BR(K) = BR(K)
* + WDOTCR(N-K,AR(K+1,K),AI(K+1,K),1,BR(K+1),BI(K+1),1)
BI(K) = BI(K)
* + WDOTCI(N-K,AR(K+1,K),AI(K+1,K),1,BR(K+1),BI(K+1),1)
L = IPVT(K)
IF (L .EQ. K) GO TO 70
TR = BR(L)
TI = BI(L)
BR(L) = BR(K)
BI(L) = BI(K)
BR(K) = TR
BI(K) = TI
70 CONTINUE
80 CONTINUE
90 CONTINUE
100 CONTINUE
RETURN
END
SUBROUTINE WGEDI(AR,AI,LDA,N,IPVT,DETR,DETI,WORKR,WORKI,JOB)
INTEGER LDA,N,IPVT(1),JOB
DOUBLE PRECISION AR(LDA,1),AI(LDA,1),DETR(2),DETI(2),WORKR(1),
* WORKI(1)
C
C WGEDI COMPUTES THE DETERMINANT AND INVERSE OF A MATRIX
C USING THE FACTORS COMPUTED BY WGECO OR WGEFA.
C
C ON ENTRY
C
C A DOUBLE-COMPLEX(LDA, N)
C THE OUTPUT FROM WGECO OR WGEFA.
C
C LDA INTEGER
C THE LEADING DIMENSION OF THE ARRAY A .
C
C N INTEGER
C THE ORDER OF THE MATRIX A .
C
C IPVT INTEGER(N)
C THE PIVOT VECTOR FROM WGECO OR WGEFA.
C
C WORK DOUBLE-COMPLEX(N)
C WORK VECTOR. CONTENTS DESTROYED.
C
C JOB INTEGER
C = 11 BOTH DETERMINANT AND INVERSE.
C = 01 INVERSE ONLY.
C = 10 DETERMINANT ONLY.
C
C ON RETURN
C
C A INVERSE OF ORIGINAL MATRIX IF REQUESTED.
C OTHERWISE UNCHANGED.
C
C DET DOUBLE-COMPLEX(2)
C DETERMINANT OF ORIGINAL MATRIX IF REQUESTED.
C OTHERWISE NOT REFERENCED.
C DETERMINANT = DET(1) * 10.0**DET(2)
C WITH 1.0 .LE. CABS1(DET(1) .LT. 10.0
C OR DET(1) .EQ. 0.0 .
C
C ERROR CONDITION
C
C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS
C A ZERO ON THE DIAGONAL AND THE INVERSE IS REQUESTED.
C IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED CORRECTLY
C AND IF WGECO HAS SET RCOND .GT. 0.0 OR WGEFA HAS SET
C INFO .EQ. 0 .
C
C LINPACK. THIS VERSION DATED 07/01/79 .
C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C SUBROUTINES AND FUNCTIONS
C
C BLAS WAXPY,WSCAL,WSWAP
C FORTRAN DABS,MOD
C
C INTERNAL VARIABLES
C
DOUBLE PRECISION TR,TI
DOUBLE PRECISION TEN
INTEGER I,J,K,KB,KP1,L,NM1
C
DOUBLE PRECISION ZDUMR,ZDUMI
DOUBLE PRECISION CABS1
CABS1(ZDUMR,ZDUMI) = DABS(ZDUMR) + DABS(ZDUMI)
C
C COMPUTE DETERMINANT
C
IF (JOB/10 .EQ. 0) GO TO 80
DETR(1) = 1.0D0
DETI(1) = 0.0D0
DETR(2) = 0.0D0
DETI(2) = 0.0D0
TEN = 10.0D0
DO 60 I = 1, N
IF (IPVT(I) .EQ. I) GO TO 10
DETR(1) = -DETR(1)
DETI(1) = -DETI(1)
10 CONTINUE
CALL WMUL(AR(I,I),AI(I,I),DETR(1),DETI(1),DETR(1),DETI(1))
C ...EXIT
C ...EXIT
IF (CABS1(DETR(1),DETI(1)) .EQ. 0.0D0) GO TO 70
20 IF (CABS1(DETR(1),DETI(1)) .GE. 1.0D0) GO TO 30
DETR(1) = TEN*DETR(1)
DETI(1) = TEN*DETI(1)
DETR(2) = DETR(2) - 1.0D0
DETI(2) = DETI(2) - 0.0D0
GO TO 20
30 CONTINUE
40 IF (CABS1(DETR(1),DETI(1)) .LT. TEN) GO TO 50
DETR(1) = DETR(1)/TEN
DETI(1) = DETI(1)/TEN
DETR(2) = DETR(2) + 1.0D0
DETI(2) = DETI(2) + 0.0D0
GO TO 40
50 CONTINUE
60 CONTINUE
70 CONTINUE
80 CONTINUE
C
C COMPUTE INVERSE(U)
C
IF (MOD(JOB,10) .EQ. 0) GO TO 160
DO 110 K = 1, N
CALL WDIV(1.0D0,0.0D0,AR(K,K),AI(K,K),AR(K,K),AI(K,K))
TR = -AR(K,K)
TI = -AI(K,K)
CALL WSCAL(K-1,TR,TI,AR(1,K),AI(1,K),1)
KP1 = K + 1
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.