home *** CD-ROM | disk | FTP | other *** search
/ Crawly Crypt Collection 1 / crawlyvol1.bin / apps / math / matlab / src / mat3.f < prev    next >
Text File  |  1992-04-21  |  29KB  |  944 lines
  1.       SUBROUTINE GETSYM
  2. C     GET A SYMBOL
  3.       DOUBLE PRECISION STKR(50005),STKI(50005)
  4.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
  5.       INTEGER ALFA(52),ALFB(52),ALFL,CASE
  6.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  7.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  8.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
  9.       COMMON /ALFS/ ALFA,ALFB,ALFL,CASE
  10.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  11.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  12.       DOUBLE PRECISION SYV,S,FLOP
  13.       INTEGER BLANK,Z,DOT,D,E,PLUS,MINUS,NAME,NUM,SIGN,CHCNT,EOL
  14.       INTEGER STAR,SLASH,BSLASH,SS
  15.       CHARACTER MYCHAR, CALFA
  16.       DATA BLANK/36/,Z/35/,DOT/47/,D/13/,E/14/,EOL/99/,PLUS/41/
  17.       DATA MINUS/42/,NAME/1/,NUM/0/,STAR/43/,SLASH/44/,BSLASH/45/
  18.    10 IF (CHAR .NE. BLANK) GO TO 20
  19.       CALL GETCH
  20.       GO TO 10
  21.    20 LPT(2) = LPT(3)
  22.       LPT(3) = LPT(4)
  23.       IF (CHAR .LE. 9) GO TO 50
  24.       IF (CHAR .LE. Z) GO TO 30
  25. C
  26. C     SPECIAL CHARACTER
  27.       SS = SYM
  28.       SYM = CHAR
  29.       CALL GETCH
  30.       IF (SYM .NE. DOT) GO TO 90
  31. C
  32. C     IS DOT PART OF NUMBER OR OPERATOR
  33.       SYV = 0.0D0
  34.       IF (CHAR .LE. 9) GO TO 55
  35.       IF (CHAR.EQ.STAR .OR. CHAR.EQ.SLASH .OR. CHAR.EQ.BSLASH) GO TO 90
  36.       IF (SS.EQ.STAR .OR. SS.EQ.SLASH .OR. SS.EQ.BSLASH) GO TO 90
  37.       GO TO 55
  38. C
  39. C     NAME
  40.    30 SYM = NAME
  41.       SYN(1) = CHAR
  42.       CHCNT = 1
  43.    40 CALL GETCH
  44.       CHCNT = CHCNT+1
  45.       IF (CHAR .GT. Z) GO TO 45
  46.       IF (CHCNT .LE. 4) SYN(CHCNT) = CHAR
  47.       GO TO 40
  48.    45 IF (CHCNT .GT. 4) GO TO 47
  49.       DO 46 I = CHCNT, 4
  50.    46 SYN(I) = BLANK
  51.    47 CONTINUE
  52.       GO TO 90
  53. C
  54. C     NUMBER
  55.    50 CALL GETVAL(SYV)
  56.       IF (CHAR .NE. DOT) GO TO 60
  57.       CALL GETCH
  58.    55 CHCNT = LPT(4)
  59.       CALL GETVAL(S)
  60.       CHCNT = LPT(4) - CHCNT
  61.       IF (CHAR .EQ. EOL) CHCNT = CHCNT+1
  62.       SYV = SYV + S/10.0D0**CHCNT
  63.    60 IF (CHAR.NE.D .AND. CHAR.NE.E) GO TO 70
  64.       CALL GETCH
  65.       SIGN = CHAR
  66.       IF (SIGN.EQ.MINUS .OR. SIGN.EQ.PLUS) CALL GETCH
  67.       CALL GETVAL(S)
  68.       IF (SIGN .NE. MINUS) SYV = SYV*10.0D0**S
  69.       IF (SIGN .EQ. MINUS) SYV = SYV/10.0D0**S
  70.    70 STKI(VSIZE) = FLOP(SYV)
  71.       SYM = NUM
  72. C
  73.    90 IF (CHAR .NE. BLANK) GO TO 99
  74.       CALL GETCH
  75.       GO TO 90
  76.    99 IF (DDT .NE. 1) RETURN
  77.       CALFA=MYCHAR(ALFA(SYM+1))
  78.       IF (SYM.GT.NAME .AND. SYM.LT.ALFL) WRITE(WTE,197) CALFA
  79.       IF (SYM .GE. ALFL) WRITE(WTE,198)
  80.       IF (SYM .EQ. NAME) CALL PRNTID(SYN,1)
  81.       IF (SYM .EQ. NUM) WRITE(WTE,199) SYV
  82.   197 FORMAT(0X,A1)
  83.   198 FORMAT(0X,'EOL')
  84.   199 FORMAT(0X,G8.2)
  85.       RETURN
  86.       END
  87.       SUBROUTINE GETLIN
  88. C     GET A NEW LINE
  89.       INTEGER ALFA(52),ALFB(52),ALFL,CASE
  90.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  91.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  92.       CHARACTER CBUF(256),MYCHAR
  93.       COMMON /ALFS/ ALFA,ALFB,ALFL,CASE
  94.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  95.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  96.       INTEGER LRECL,EOL,SLASH,BSLASH,DOT,BLANK,RETU(4)
  97.       DATA EOL/99/,DOT/47/,BLANK/36/,RETU/27,14,29,30/
  98.       DATA SLASH/44/,BSLASH/45/,LRECL/80/
  99. C
  100.    10 L = LPT(1)
  101.    11 DO 12 J = 1, LRECL
  102.          CBUF(J)=' '
  103.    12 CONTINUE
  104.       READ(RIO,101,END=50,ERR=10) (CBUF(J),J=1,LRECL)
  105.       DO 1000 J=1,LRECL
  106.         BUF(J)=ICHAR(CBUF(J))
  107.  1000 CONTINUE
  108.   101 FORMAT(80A1)
  109.       N = LRECL+1
  110.    15 N = N-1
  111.       IF (BUF(N) .EQ. ALFA(BLANK+1)) GO TO 15
  112.       DO 1001 J=1,N
  113.         CBUF(J)=MYCHAR(BUF(J))
  114.  1001 CONTINUE
  115.       IF (MOD(LCT(4),2) .EQ. 1) WRITE(WTE,102) (CBUF(J),J=1,N)
  116.       IF (WIO .NE. 0) WRITE(WIO,102) (CBUF(J),J=1,N)
  117.   102 FORMAT(0X,80A1)
  118. C
  119.       DO 40 J = 1, N
  120.          DO 20 K = 1, ALFL
  121.            IF (BUF(J).EQ.ALFA(K) .OR. BUF(J).EQ.ALFB(K)) GO TO 30
  122.    20    CONTINUE
  123.          K = EOL+1
  124.          CALL XCHAR(BUF(J),K)
  125.          IF (K .GT. EOL) GO TO 10
  126.          IF (K .EQ. EOL) GO TO 45
  127.          IF (K .EQ. -1) L = L-1
  128.          IF (K .LE. 0) GO TO 40
  129. C
  130.    30    K = K-1
  131.          IF (K.EQ.SLASH .AND. BUF(J+1).EQ.BUF(J)) GO TO 45
  132.          IF (K.EQ.DOT .AND. BUF(J+1).EQ.BUF(J)) GO TO 11
  133.          IF (K.EQ.BSLASH .AND. N.EQ.1) GO TO 60
  134.          LIN(L) = K
  135.          IF (L.LT.1024) L = L+1
  136.          IF (L.EQ.1024) WRITE(WTE,33) L
  137.    33    FORMAT(0X,'INPUT BUFFER LIMIT IS ',I4,' CHARACTERS.')
  138.    40 CONTINUE
  139.    45 LIN(L) = EOL
  140.       LPT(6) = L
  141.       LPT(4) = LPT(1)
  142.       LPT(3) = 0
  143.       LPT(2) = 0
  144.       LCT(1) = 0
  145.       CALL GETCH
  146.       RETURN
  147. C
  148.    50 IF (RIO .EQ. RTE) GO TO 52
  149.       CALL PUTID(LIN(L),RETU)
  150.       L = L + 4 
  151.       GO TO 45 
  152.    52 CALL FILES(-RTE,BUF,BUF)
  153.       LIN(L) = EOL
  154.       RETURN
  155. C
  156.    60 N = LPT(6) - LPT(1)
  157.       DO 61 I = 1, N
  158.          J = L+I-1
  159.          K = LIN(J)
  160.          BUF(I) = ALFA(K+1)
  161.          IF (CASE.EQ.1 .AND. K.LT.36) BUF(I) = ALFB(K+1)
  162.    61 CONTINUE
  163.       CALL EDIT(BUF,N)
  164.       N = N + 1
  165.       GO TO 15
  166.       END
  167.       SUBROUTINE GETCH
  168. C     GET NEXT CHARACTER
  169.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  170.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  171.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  172.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  173.       INTEGER EOL
  174.       DATA EOL/99/
  175.       L = LPT(4)
  176.       CHAR = LIN(L)
  177.       IF (CHAR .NE. EOL) LPT(4) = L + 1
  178.       RETURN
  179.       END
  180.       SUBROUTINE GETVAL(S)
  181.       DOUBLE PRECISION S
  182. C     FORM NUMERICAL VALUE
  183.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  184.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  185.       S = 0.0D0
  186.    10 IF (CHAR .GT. 9) RETURN
  187.       S = 10.0D0*S + CHAR
  188.       CALL GETCH
  189.       GO TO 10
  190.       END
  191.       SUBROUTINE MATFN1
  192. C
  193. C     EVALUATE FUNCTIONS INVOLVING GAUSSIAN ELIMINATION
  194. C
  195.       DOUBLE PRECISION STKR(50005),STKI(50005)
  196.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
  197.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  198.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  199.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
  200.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  201.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  202.       DOUBLE PRECISION DTR(2),DTI(2),SR,SI,RCOND,T,T0,T1,FLOP,EPS,WASUM
  203. C
  204.       IF (DDT .EQ. 1) WRITE(WTE,100) FIN
  205.   100 FORMAT(0X,'MATFN1',I4)
  206. C
  207.       L = LSTK(TOP)
  208.       M = MSTK(TOP)
  209.       N = NSTK(TOP)
  210.       IF (FIN .EQ. -1) GO TO 10
  211.       IF (FIN .EQ. -2) GO TO 20
  212.       GO TO (30,40,50,60,70,80,85),FIN
  213. C
  214. C     MATRIX RIGHT DIVISION, A/A2
  215.    10 L2 = LSTK(TOP+1)
  216.       M2 = MSTK(TOP+1)
  217.       N2 = NSTK(TOP+1)
  218.       IF (M2 .NE. N2) CALL ERROR(20)
  219.       IF (ERR .GT. 0) RETURN
  220.       IF (M*N .EQ. 1) GO TO 16
  221.       IF (N .NE. N2) CALL ERROR(11)
  222.       IF (ERR .GT. 0) RETURN
  223.       L3 = L2 + M2*N2
  224.       ERR = L3+N2 - LSTK(BOT)
  225.       IF (ERR .GT. 0) CALL ERROR(17)
  226.       IF (ERR .GT. 0) RETURN
  227.       CALL WGECO(STKR(L2),STKI(L2),M2,N2,BUF,RCOND,STKR(L3),STKI(L3))
  228.       IF (RCOND .EQ. 0.0D0) CALL ERROR(19)
  229.       IF (ERR .GT. 0) RETURN
  230.       T = FLOP(1.0D0 + RCOND)
  231.       IF (T.EQ.1.0D0 .AND. FUN.NE.21) WRITE(WTE,11) RCOND
  232.       IF (T.EQ.1.0D0 .AND. FUN.NE.21 .AND. WIO.NE.0) WRITE(WIO,11) RCOND
  233.    11 FORMAT(0X,'WARNING.'
  234.      $  /0X,'MATRIX IS CLOSE TO SINGULAR OR BADLY SCALED.'
  235.      $  /0X,'RESULTS MAY BE INACCURATE.  RCOND =', 1PD13.4/)
  236.       IF (T.EQ.1.0D0 .AND. FUN.EQ.21) WRITE(WTE,12) RCOND
  237.       IF (T.EQ.1.0D0 .AND. FUN.EQ.21 .AND. WIO.NE.0) WRITE(WIO,12) RCOND
  238.    12 FORMAT(0X,'WARNING.'
  239.      $  /0X,'EIGENVECTORS ARE BADLY CONDITIONED.'
  240.      $  /0X,'RESULTS MAY BE INACCURATE.  RCOND =', 1PD13.4/)
  241.       DO 15 I = 1, M
  242.          DO 13 J = 1, N
  243.             LS = L+I-1+(J-1)*M
  244.             LL = L3+J-1
  245.             STKR(LL) = STKR(LS)
  246.             STKI(LL) = -STKI(LS)
  247.    13    CONTINUE
  248.          CALL WGESL(STKR(L2),STKI(L2),M2,N2,BUF,STKR(L3),STKI(L3),1)
  249.          DO 14 J = 1, N
  250.             LL = L+I-1+(J-1)*M
  251.             LS = L3+J-1
  252.             STKR(LL) = STKR(LS)
  253.             STKI(LL) = -STKI(LS)
  254.    14    CONTINUE
  255.    15 CONTINUE
  256.       IF (FUN .NE. 21) GO TO 99
  257. C
  258. C     CHECK FOR IMAGINARY ROUNDOFF IN MATRIX FUNCTIONS
  259.       SR = WASUM(N*N,STKR(L),STKR(L),1)
  260.       SI = WASUM(N*N,STKI(L),STKI(L),1)
  261.       EPS = STKR(VSIZE-4)
  262.       T = EPS*SR
  263.       IF (DDT .EQ. 18) WRITE(WTE,115) SR,SI,EPS,T
  264.   115 FORMAT(0X,'SR,SI,EPS,T',1P4D13.4)
  265.       IF (SI .LE. EPS*SR) CALL RSET(N*N,0.0D0,STKI(L),1)
  266.       GO TO 99
  267. C
  268.    16 SR = STKR(L)
  269.       SI = STKI(L)
  270.       N = N2
  271.       M = N
  272.       MSTK(TOP) = N
  273.       NSTK(TOP) = N
  274.       CALL WCOPY(N*N,STKR(L2),STKI(L2),1,STKR(L),STKI(L),1)
  275.       GO TO 30
  276. C
  277. C     MATRIX LEFT DIVISION A BACKSLASH A2
  278.    20 L2 = LSTK(TOP+1)
  279.       M2 = MSTK(TOP+1)
  280.       N2 = NSTK(TOP+1)
  281.       IF (M .NE. N) CALL ERROR(20)
  282.       IF (ERR .GT. 0) RETURN
  283.       IF (M2*N2 .EQ. 1) GO TO 26
  284.       L3 = L2 + M2*N2
  285.       ERR = L3+N - LSTK(BOT)
  286.       IF (ERR .GT. 0) CALL ERROR(17)
  287.       IF (ERR .GT. 0) RETURN
  288.       CALL WGECO(STKR(L),STKI(L),M,N,BUF,RCOND,STKR(L3),STKI(L3))
  289.       IF (RCOND .EQ. 0.0D0) CALL ERROR(19)
  290.       IF (ERR .GT. 0) RETURN
  291.       T = FLOP(1.0D0 + RCOND)
  292.       IF (T .EQ. 1.0D0) WRITE(WTE,11) RCOND
  293.       IF (T.EQ.1.0D0 .AND. WIO.NE.0) WRITE(WIO,11) RCOND
  294.       IF (M2 .NE. N) CALL ERROR(12)
  295.       IF (ERR .GT. 0) RETURN
  296.       DO 23 J = 1, N2
  297.          LJ = L2+(J-1)*M2
  298.          CALL WGESL(STKR(L),STKI(L),M,N,BUF,STKR(LJ),STKI(LJ),0)
  299.    23 CONTINUE
  300.       NSTK(TOP) = N2
  301.       CALL WCOPY(M2*N2,STKR(L2),STKI(L2),1,STKR(L),STKI(L),1)
  302.       GO TO 99
  303.    26 SR = STKR(L2)
  304.       SI = STKI(L2)
  305.       GO TO 30
  306. C
  307. C     INV
  308. C
  309.    30 IF (M .NE. N) CALL ERROR(20)
  310.       IF (ERR .GT. 0) RETURN
  311.       IF (DDT .EQ. 17) GO TO 32
  312.       DO 31 J = 1, N
  313.       DO 31 I = 1, N
  314.         LS = L+I-1+(J-1)*N
  315.         T0 = STKR(LS)
  316.         T1 = I+J-1
  317.         T1 = FLOP(1.0D0/T1)
  318.         IF (T0 .NE. T1) GO TO 32
  319.    31 CONTINUE
  320.       GO TO 72
  321.    32 L3 = L + N*N
  322.       ERR = L3+N - LSTK(BOT)
  323.       IF (ERR .GT. 0) CALL ERROR(17)
  324.       IF (ERR .GT. 0) RETURN
  325.       CALL WGECO(STKR(L),STKI(L),M,N,BUF,RCOND,STKR(L3),STKI(L3))
  326.       IF (RCOND .EQ. 0.0D0) CALL ERROR(19)
  327.       IF (ERR .GT. 0) RETURN
  328.       T = FLOP(1.0D0 + RCOND)
  329.       IF (T .EQ. 1.0D0) WRITE(WTE,11) RCOND
  330.       IF (T.EQ.1.0D0 .AND. WIO.NE.0) WRITE(WIO,11) RCOND
  331.       CALL WGEDI(STKR(L),STKI(L),M,N,BUF,DTR,DTI,STKR(L3),STKI(L3),1)
  332.       IF (FIN .LT. 0) CALL WSCAL(N*N,SR,SI,STKR(L),STKI(L),1)
  333.       GO TO 99
  334. C
  335. C     DET
  336. C
  337.    40 IF (M .NE. N) CALL ERROR(20)
  338.       IF (ERR .GT. 0) RETURN
  339.       CALL WGEFA(STKR(L),STKI(L),M,N,BUF,INFO)
  340.       CALL WGEDI(STKR(L),STKI(L),M,N,BUF,DTR,DTI,SR,SI,10)
  341.       K = IDINT(DTR(2))
  342.       KA = IABS(K)+2
  343.       T = 1.0D0
  344.       DO 41 I = 1, KA
  345.          T = T/10.0D0
  346.          IF (T .EQ. 0.0D0) GO TO 42
  347.    41 CONTINUE
  348.       STKR(L) = DTR(1)*10.D0**K
  349.       STKI(L) = DTI(1)*10.D0**K
  350.       MSTK(TOP) = 1
  351.       NSTK(TOP) = 1
  352.       GO TO 99
  353.    42 IF (DTI(1) .EQ. 0.0D0) WRITE(WTE,43) DTR(1),K
  354.       IF (DTI(1) .NE. 0.0D0) WRITE(WTE,44) DTR(1),DTI(1),K
  355.    43 FORMAT(0X,'DET =  ',F7.4,7H * 10**,I4)
  356.    44 FORMAT(0X,'DET =  ',F7.4,' + ',F7.4,' i ',7H * 10**,I4)
  357.       STKR(L) = DTR(1)
  358.       STKI(L) = DTI(1)
  359.       STKR(L+1) = DTR(2)
  360.       STKI(L+1) = 0.0D0
  361.       MSTK(TOP) = 1
  362.       NSTK(TOP) = 2
  363.       GO TO 99
  364. C
  365. C     RCOND
  366. C
  367.    50 IF (M .NE. N) CALL ERROR(20)
  368.       IF (ERR .GT. 0) RETURN
  369.       L3 = L + N*N
  370.       ERR = L3+N - LSTK(BOT)
  371.       IF (ERR .GT. 0) CALL ERROR(17)
  372.       IF (ERR .GT. 0) RETURN
  373.       CALL WGECO(STKR(L),STKI(L),M,N,BUF,RCOND,STKR(L3),STKI(L3))
  374.       STKR(L) = RCOND
  375.       STKI(L) = 0.0D0
  376.       MSTK(TOP) = 1
  377.       NSTK(TOP) = 1
  378.       IF (LHS .EQ. 1) GO TO 99
  379.       L = L + 1
  380.       CALL WCOPY(N,STKR(L3),STKI(L3),1,STKR(L),STKI(L),1)
  381.       TOP = TOP + 1
  382.       LSTK(TOP) = L
  383.       MSTK(TOP) = N
  384.       NSTK(TOP) = 1
  385.       GO TO 99
  386. C
  387. C     LU
  388. C
  389.    60 IF (M .NE. N) CALL ERROR(20)
  390.       IF (ERR .GT. 0) RETURN
  391.       CALL WGEFA(STKR(L),STKI(L),M,N,BUF,INFO)
  392.       IF (LHS .NE. 2) GO TO 99
  393.       NN = N*N
  394.       IF (TOP+1 .GE. BOT) CALL ERROR(18)
  395.       IF (ERR .GT. 0) RETURN
  396.       TOP = TOP+1
  397.       LSTK(TOP) = L + NN
  398.       MSTK(TOP) = N
  399.       NSTK(TOP) = N
  400.       ERR = L+NN+NN - LSTK(BOT)
  401.       IF (ERR .GT. 0) CALL ERROR(17)
  402.       IF (ERR .GT. 0) RETURN
  403.       DO 64 KB = 1, N
  404.         K = N+1-KB
  405.         DO 61 I = 1, N
  406.           LL = L+I-1+(K-1)*N
  407.           LU = LL + NN
  408.           IF (I .LE. K) STKR(LU) = STKR(LL)
  409.           IF (I .LE. K) STKI(LU) = STKI(LL)
  410.           IF (I .GT. K) STKR(LU) = 0.0D0
  411.           IF (I .GT. K) STKI(LU) = 0.0D0
  412.           IF (I .LT. K) STKR(LL) = 0.0D0
  413.           IF (I .LT. K) STKI(LL) = 0.0D0
  414.           IF (I .EQ. K) STKR(LL) = 1.0D0
  415.           IF (I .EQ. K) STKI(LL) = 0.0D0
  416.           IF (I .GT. K) STKR(LL) = -STKR(LL)
  417.           IF (I .GT. K) STKI(LL) = -STKI(LL)
  418.    61   CONTINUE
  419.         I = BUF(K)
  420.         IF (I .EQ. K) GO TO 64
  421.         LI = L+I-1+(K-1)*N
  422.         LK = L+K-1+(K-1)*N
  423.         CALL WSWAP(N-K+1,STKR(LI),STKI(LI),N,STKR(LK),STKI(LK),N)
  424.    64 CONTINUE
  425.       GO TO 99
  426. C
  427. C     HILBERT
  428.    70 N = IDINT(STKR(L))
  429.       MSTK(TOP) = N
  430.       NSTK(TOP) = N
  431.    72 CALL HILBER(STKR(L),N,N)
  432.       CALL RSET(N*N,0.0D0,STKI(L),1)
  433.       IF (FIN .LT. 0) CALL WSCAL(N*N,SR,SI,STKR(L),STKI(L),1)
  434.       GO TO 99
  435. C
  436. C     CHOLESKY
  437.    80 IF (M .NE. N) CALL ERROR(20)
  438.       IF (ERR .GT. 0) RETURN
  439.       CALL WPOFA(STKR(L),STKI(L),M,N,ERR)
  440.       IF (ERR .NE. 0) CALL ERROR(29)
  441.       IF (ERR .GT. 0) RETURN
  442.       DO 81 J = 1, N
  443.         LL = L+J+(J-1)*M
  444.         CALL WSET(M-J,0.0D0,0.0D0,STKR(LL),STKI(LL),1)
  445.    81 CONTINUE
  446.       GO TO 99
  447. C
  448. C     RREF
  449.    85 IF (RHS .LT. 2) GO TO 86
  450.         TOP = TOP-1
  451.         L = LSTK(TOP)
  452.         IF (MSTK(TOP) .NE. M) CALL ERROR(5)
  453.         IF (ERR .GT. 0) RETURN
  454.         N = N + NSTK(TOP)
  455.    86 CALL RREF(STKR(L),STKI(L),M,M,N,STKR(VSIZE-4))
  456.       NSTK(TOP) = N
  457.       GO TO 99
  458. C
  459.    99 RETURN
  460.       END
  461.       SUBROUTINE MATFN2
  462. C
  463. C     EVALUATE ELEMENTARY FUNCTIONS AND FUNCTIONS INVOLVING
  464. C     EIGENVALUES AND EIGENVECTORS
  465. C
  466.       DOUBLE PRECISION STKR(50005),STKI(50005)
  467.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
  468.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  469.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  470.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
  471.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  472.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  473.       DOUBLE PRECISION PYTHAG,ROUND,TR,TI,SR,SI,POWR,POWI,FLOP
  474.       LOGICAL HERM,SCHUR,VECT,HESS
  475. C
  476.       IF (DDT .EQ. 1) WRITE(WTE,100) FIN
  477.   100 FORMAT(0X,'MATFN2',I4)
  478. C
  479. C     FUNCTIONS/FIN
  480. C     **   SIN  COS ATAN  EXP  SQRT LOG
  481. C      0    1    2    3    4    5    6
  482. C    EIG  SCHU HESS POLY ROOT
  483. C     11   12   13   14   15
  484. C    ABS  ROUN REAL IMAG CONJ
  485. C     21   22   23   24   25
  486.       IF (FIN .NE. 0) GO TO 05
  487.          L = LSTK(TOP+1)
  488.          POWR = STKR(L)
  489.          POWI = STKI(L)
  490.    05 L = LSTK(TOP)
  491.       M = MSTK(TOP)
  492.       N = NSTK(TOP)
  493.       IF (FIN .GE. 11 .AND. FIN .LE. 13) GO TO 10
  494.       IF (FIN .EQ. 14 .AND. (M.EQ.1 .OR. N.EQ.1)) GO TO 50
  495.       IF (FIN .EQ. 14) GO TO 10
  496.       IF (FIN .EQ. 15) GO TO 60
  497.       IF (FIN .GT. 20) GO TO 40
  498.       IF (M .EQ. 1 .OR. N .EQ. 1) GO TO 40
  499. C
  500. C     EIGENVALUES AND VECTORS
  501.    10 IF (M .NE. N) CALL ERROR(20)
  502.       IF (ERR .GT. 0) RETURN
  503.       SCHUR = FIN .EQ. 12
  504.       HESS = FIN .EQ. 13
  505.       VECT = LHS.EQ.2 .OR. FIN.LT.10
  506.       NN = N*N
  507.       L2 = L + NN
  508.       LD = L2 + NN
  509.       LE = LD + N
  510.       LW = LE + N
  511.       ERR = LW+N - LSTK(BOT)
  512.       IF (ERR .GT. 0) CALL ERROR(17)
  513.       IF (ERR .GT. 0) RETURN
  514.       CALL WCOPY(NN,STKR(L),STKI(L),1,STKR(L2),STKI(L2),1)
  515. C
  516. C     CHECK IF HERMITIAN
  517.       DO 15 J = 1, N
  518.       DO 15 I = 1, J
  519.          LS = L+I-1+(J-1)*N
  520.          LL = L+(I-1)*N+J-1
  521.          HERM = STKR(LL).EQ.STKR(LS) .AND. STKI(LL).EQ.-STKI(LS)
  522.          IF (.NOT. HERM) GO TO 30
  523.    15 CONTINUE
  524. C
  525. C     HERMITIAN EIGENVALUE PROBLEM
  526.       CALL WSET(NN,0.0D0,0.0D0,STKR(L),STKI(L),1)
  527.       CALL WSET(N,1.0D0,0.0D0,STKR(L),STKI(L),N+1)
  528.       CALL WSET(N,0.0D0,0.0D0,STKI(LD),STKI(LE),1)
  529.       JOB = 0
  530.       IF (VECT) JOB = 1
  531.       CALL HTRIDI(N,N,STKR(L2),STKI(L2),STKR(LD),STKR(LE),
  532.      $            STKR(LE),STKR(LW))
  533.       IF (.NOT.HESS) CALL IMTQL2(N,N,STKR(LD),STKR(LE),STKR(L),ERR,JOB)
  534.       IF (ERR .GT. 0) CALL ERROR(24)
  535.       IF (ERR .GT. 0) RETURN
  536.       IF (JOB .NE. 0)
  537.      $  CALL HTRIBK(N,N,STKR(L2),STKI(L2),STKR(LW),N,STKR(L),STKI(L))
  538.       GO TO 31
  539. C
  540. C     NON-HERMITIAN EIGENVALUE PROBLEM
  541.    30 CALL CORTH(N,N,1,N,STKR(L2),STKI(L2),STKR(LW),STKI(LW))
  542.       IF (.NOT.VECT .AND. HESS) GO TO 31
  543.       JOB = 0
  544.       IF (VECT) JOB = 2
  545.       IF (VECT .AND. SCHUR) JOB = 1
  546.       IF (HESS) JOB = 3
  547.       CALL COMQR3(N,N,1,N,STKR(LW),STKI(LW),STKR(L2),STKI(L2),
  548.      $            STKR(LD),STKI(LD),STKR(L),STKI(L),ERR,JOB)
  549.       IF (ERR .GT. 0) CALL ERROR(24)
  550.       IF (ERR .GT. 0) RETURN
  551. C
  552. C     VECTORS
  553.    31 IF (.NOT.VECT) GO TO 34
  554.       IF (TOP+1 .GE. BOT) CALL ERROR(18)
  555.       IF (ERR .GT. 0) RETURN
  556.       TOP = TOP+1
  557.       LSTK(TOP) = L2
  558.       MSTK(TOP) = N
  559.       NSTK(TOP) = N
  560. C
  561. C     DIAGONAL OF VALUES OR CANONICAL FORMS
  562.    34 IF (.NOT.VECT .AND. .NOT.SCHUR .AND. .NOT.HESS) GO TO 37
  563.       DO 36 J = 1, N
  564.          LJ = L2+(J-1)*N
  565.          IF (SCHUR .AND. (.NOT.HERM)) LJ = LJ+J
  566.          IF (HESS .AND. (.NOT.HERM)) LJ = LJ+J+1
  567.          LL = L2+J*N-LJ
  568.          CALL WSET(LL,0.0D0,0.0D0,STKR(LJ),STKI(LJ),1)
  569.    36 CONTINUE
  570.       IF (.NOT.HESS .OR. HERM)
  571.      $   CALL WCOPY(N,STKR(LD),STKI(LD),1,STKR(L2),STKI(L2),N+1)
  572.       LL = L2+1
  573.       IF (HESS .AND. HERM)
  574.      $   CALL WCOPY(N-1,STKR(LE+1),STKI(LE+1),1,STKR(LL),STKI(LL),N+1)
  575.       LL = L2+N
  576.       IF (HESS .AND. HERM)
  577.      $   CALL WCOPY(N-1,STKR(LE+1),STKI(LE+1),1,STKR(LL),STKI(LL),N+1)
  578.       IF (FIN .LT. 10) GO TO 42
  579.       IF (VECT .OR. .NOT.(SCHUR.OR.HESS)) GO TO 99
  580.       CALL WCOPY(NN,STKR(L2),STKI(L2),1,STKR(L),STKI(L),1)
  581.       GO TO 99
  582. C
  583. C     VECTOR OF EIGENVALUES
  584.    37 IF (FIN .EQ. 14) GO TO 52
  585.       CALL WCOPY(N,STKR(LD),STKI(LD),1,STKR(L),STKI(L),1)
  586.       NSTK(TOP) = 1
  587.       GO TO 99
  588. C
  589. C     ELEMENTARY FUNCTIONS
  590. C     FOR MATRICES.. X,D = EIG(A), FUN(A) = X*FUN(D)/X
  591.    40 INC = 1
  592.       N = M*N
  593.       L2 = L
  594.       GO TO 44
  595.    42 INC = N+1
  596.    44 DO 46 J = 1, N
  597.         LS = L2+(J-1)*INC
  598.         SR = STKR(LS)
  599.         SI = STKI(LS)
  600.         TI = 0.0D0
  601.         IF (FIN .NE. 0) GO TO 45
  602.           CALL WLOG(SR,SI,SR,SI)
  603.           CALL WMUL(SR,SI,POWR,POWI,SR,SI)
  604.           TR = DEXP(SR)*DCOS(SI)
  605.           TI = DEXP(SR)*DSIN(SI)
  606.    45   IF (FIN .EQ. 1) TR = DSIN(SR)*DCOSH(SI)
  607.         IF (FIN .EQ. 1) TI = DCOS(SR)*DSINH(SI)
  608.         IF (FIN .EQ. 2) TR = DCOS(SR)*DCOSH(SI)
  609.         IF (FIN .EQ. 2) TI = -DSIN(SR)*DSINH(SI)
  610.         IF (FIN .EQ. 3) CALL WATAN(SR,SI,TR,TI)
  611.         IF (FIN .EQ. 4) TR = DEXP(SR)*DCOS(SI)
  612.         IF (FIN .EQ. 4) TI = DEXP(SR)*DSIN(SI)
  613.         IF (FIN .EQ. 5) CALL WSQRT(SR,SI,TR,TI)
  614.         IF (FIN .EQ. 6) CALL WLOG(SR,SI,TR,TI)
  615.         IF (FIN .EQ. 21) TR = PYTHAG(SR,SI)
  616.         IF (FIN .EQ. 22) TR = ROUND(SR)
  617.         IF (FIN .EQ. 23) TR = SR
  618.         IF (FIN .EQ. 24) TR = SI
  619.         IF (FIN .EQ. 25) TR = SR
  620.         IF (FIN .EQ. 25) TI = -SI
  621.         IF (ERR .GT. 0) RETURN
  622.         STKR(LS) = FLOP(TR)
  623.         STKI(LS) = 0.0D0
  624.         IF (TI .NE. 0.0D0) STKI(LS) = FLOP(TI)
  625.    46 CONTINUE
  626.       IF (INC .EQ. 1) GO TO 99
  627.       DO 48 J = 1, N
  628.         LS = L2+(J-1)*INC
  629.         SR = STKR(LS)
  630.         SI = STKI(LS)
  631.         LS = L+(J-1)*N
  632.         LL = L2+(J-1)*N
  633.         CALL WCOPY(N,STKR(LS),STKI(LS),1,STKR(LL),STKI(LL),1)
  634.         CALL WSCAL(N,SR,SI,STKR(LS),STKI(LS),1)
  635.    48 CONTINUE
  636. C     SIGNAL MATFN1 TO DIVIDE BY EIGENVECTORS
  637.       FUN = 21
  638.       FIN = -1
  639.       TOP = TOP-1
  640.       GO TO 99
  641. C
  642. C     POLY
  643. C     FORM POLYNOMIAL WITH GIVEN VECTOR AS ROOTS
  644.    50 N = MAX0(M,N)
  645.       LD = L+N+1
  646.       CALL WCOPY(N,STKR(L),STKI(L),1,STKR(LD),STKI(LD),1)
  647. C
  648. C     FORM CHARACTERISTIC POLYNOMIAL
  649.    52 CALL WSET(N+1,0.0D0,0.0D0,STKR(L),STKI(L),1)
  650.       STKR(L) = 1.0D0
  651.       DO 56 J = 1, N
  652.          CALL WAXPY(J,-STKR(LD),-STKI(LD),STKR(L),STKI(L),-1,
  653.      $              STKR(L+1),STKI(L+1),-1)
  654.          LD = LD+1
  655.    56 CONTINUE
  656.       MSTK(TOP) = N+1
  657.       NSTK(TOP) = 1
  658.       GO TO 99
  659. C
  660. C     ROOTS
  661.    60 LL = L+M*N
  662.       STKR(LL) = -1.0D0
  663.       STKI(LL) = 0.0D0
  664.       K = -1
  665.    61 K = K+1
  666.       L1 = L+K
  667.       IF (DABS(STKR(L1))+DABS(STKI(L1)) .EQ. 0.0D0) GO TO 61
  668.       N = MAX0(M*N - K-1, 0)
  669.       IF (N .LE. 0) GO TO 65
  670.       L2 = L1+N+1
  671.       LW = L2+N*N
  672.       ERR = LW+N - LSTK(BOT)
  673.       IF (ERR .GT. 0) CALL ERROR(17)
  674.       IF (ERR .GT. 0) RETURN
  675.       CALL WSET(N*N+N,0.0D0,0.0D0,STKR(L2),STKI(L2),1)
  676.       DO 64 J = 1, N
  677.          LL = L2+J+(J-1)*N
  678.          STKR(LL) = 1.0D0
  679.          LS = L1+J
  680.          LL = L2+(J-1)*N
  681.          CALL WDIV(-STKR(LS),-STKI(LS),STKR(L1),STKI(L1),
  682.      $             STKR(LL),STKI(LL))
  683.          IF (ERR .GT. 0) RETURN
  684.    64 CONTINUE
  685.       CALL COMQR3(N,N,1,N,STKR(LW),STKI(LW),STKR(L2),STKI(L2),
  686.      $            STKR(L),STKI(L),TR,TI,ERR,0)
  687.       IF (ERR .GT. 0) CALL ERROR(24)
  688.       IF (ERR .GT. 0) RETURN
  689.    65 MSTK(TOP) = N
  690.       NSTK(TOP) = 1
  691.       GO TO 99
  692.    99 RETURN
  693.       END
  694.       SUBROUTINE MATFN3
  695. C
  696. C     EVALUATE FUNCTIONS INVOLVING SINGULAR VALUE DECOMPOSITION
  697. C
  698.       DOUBLE PRECISION STKR(50005),STKI(50005)
  699.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
  700.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  701.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  702.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
  703.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  704.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  705.       LOGICAL FRO,INF
  706.       DOUBLE PRECISION P,S,T,TOL,EPS
  707.       DOUBLE PRECISION WDOTCR,WDOTCI,PYTHAG,WNRM2,WASUM,FLOP
  708. C
  709.       IF (DDT .EQ. 1) WRITE(WTE,100) FIN
  710.   100 FORMAT(0X,'MATFN3',I4)
  711. C
  712.       IF (FIN.EQ.1 .AND. RHS.EQ.2) TOP = TOP-1
  713.       L = LSTK(TOP)
  714.       M = MSTK(TOP)
  715.       N = NSTK(TOP)
  716.       MN = M*N
  717.       GO TO (50,70,10,30,70), FIN
  718. C
  719. C     COND
  720. C
  721.    10 LD = L + M*N
  722.       L1 = LD + MIN0(M+1,N)
  723.       L2 = L1 + N
  724.       ERR = L2+MIN0(M,N) - LSTK(BOT)
  725.       IF (ERR .GT. 0) CALL ERROR(17)
  726.       IF (ERR .GT. 0) RETURN
  727.       CALL WSVDC(STKR(L),STKI(L),M,M,N,STKR(LD),STKI(LD),
  728.      $           STKR(L1),STKI(L1),T,T,1,T,T,1,STKR(L2),STKI(L2),
  729.      $           0,ERR)
  730.       IF (ERR .NE. 0) CALL ERROR(24)
  731.       IF (ERR .GT. 0) RETURN
  732.       S = STKR(LD)
  733.       LD = LD + MIN0(M,N) - 1
  734.       T = STKR(LD)
  735.       IF (T .EQ. 0.0D0) GO TO 13
  736.       STKR(L) = FLOP(S/T)
  737.       STKI(L) = 0.0D0
  738.       MSTK(TOP) = 1
  739.       NSTK(TOP) = 1
  740.       GO TO 99
  741.    13 WRITE(WTE,14)
  742.       IF (WIO .NE. 0) WRITE(WIO,14)
  743.    14 FORMAT(0X,'CONDITION IS INFINITE')
  744.       MSTK(TOP) = 0
  745.       GO TO 99
  746. C
  747. C     NORM
  748. C
  749.    30 P = 2.0D0
  750.       INF = .FALSE.
  751.       IF (RHS .NE. 2) GO TO 31
  752.       FRO = IDINT(STKR(L)).EQ.15 .AND. MN.GT.1
  753.       INF = IDINT(STKR(L)).EQ.18 .AND. MN.GT.1
  754.       IF (.NOT. FRO) P = STKR(L)
  755.       TOP = TOP-1
  756.       L = LSTK(TOP)
  757.       M = MSTK(TOP)
  758.       N = NSTK(TOP)
  759.       MN = M*N
  760.       IF (FRO) M = MN
  761.       IF (FRO) N = 1
  762.    31 IF (M .GT. 1 .AND. N .GT. 1) GO TO 40
  763.       IF (P .EQ. 1.0D0) GO TO 36
  764.       IF (P .EQ. 2.0D0) GO TO 38
  765.       I = IWAMAX(MN,STKR(L),STKI(L),1) + L - 1
  766.       S = DABS(STKR(I)) + DABS(STKI(I))
  767.       IF (INF .OR. S .EQ. 0.0D0) GO TO 49
  768.       T = 0.0D0
  769.       DO 33 I = 1, MN
  770.          LS = L+I-1
  771.          T = FLOP(T + (PYTHAG(STKR(LS),STKI(LS))/S)**P)
  772.    33 CONTINUE
  773.       IF (P .NE. 0.0D0) P = 1.0D0/P
  774.       S = FLOP(S*T**P)
  775.       GO TO 49
  776.    36 S = WASUM(MN,STKR(L),STKI(L),1)
  777.       GO TO 49
  778.    38 S = WNRM2(MN,STKR(L),STKI(L),1)
  779.       GO TO 49
  780. C
  781. C     MATRIX NORM
  782. C
  783.    40 IF (INF) GO TO 43
  784.       IF (P .EQ. 1.0D0) GO TO 46
  785.       IF (P .NE. 2.0D0) CALL ERROR(23)
  786.       IF (ERR .GT. 0) RETURN
  787.       LD = L + M*N
  788.       L1 = LD + MIN0(M+1,N)
  789.       L2 = L1 + N
  790.       ERR = L2+MIN0(M,N) - LSTK(BOT)
  791.       IF (ERR .GT. 0) CALL ERROR(17)
  792.       IF (ERR .GT. 0) RETURN
  793.       CALL WSVDC(STKR(L),STKI(L),M,M,N,STKR(LD),STKI(LD),
  794.      $           STKR(L1),STKI(L1),T,T,1,T,T,1,STKR(L2),STKI(L2),
  795.      $           0,ERR)
  796.       IF (ERR .NE. 0) CALL ERROR(24)
  797.       IF (ERR .GT. 0) RETURN
  798.       S = STKR(LD)
  799.       GO TO 49
  800.    43 S = 0.0D0
  801.       DO 45 I = 1, M
  802.          LI = L+I-1
  803.          T = WASUM(N,STKR(LI),STKI(LI),M)
  804.          S = DMAX1(S,T)
  805.    45 CONTINUE
  806.       GO TO 49
  807.    46 S = 0.0D0
  808.       DO 48 J = 1, N
  809.          LJ = L+(J-1)*M
  810.          T = WASUM(M,STKR(LJ),STKI(LJ),1)
  811.          S = DMAX1(S,T)
  812.    48 CONTINUE
  813.       GO TO 49
  814.    49 STKR(L) = S
  815.       STKI(L) = 0.0D0
  816.       MSTK(TOP) = 1
  817.       NSTK(TOP) = 1
  818.       GO TO 99
  819. C
  820. C     SVD
  821. C
  822.    50 IF (LHS .NE. 3) GO TO 52
  823.       K = M
  824.       IF (RHS .EQ. 2) K = MIN0(M,N)
  825.       LU = L + M*N
  826.       LD = LU + M*K
  827.       LV = LD + K*N
  828.       L1 = LV + N*N
  829.       L2 = L1 + N
  830.       ERR = L2+MIN0(M,N) - LSTK(BOT)
  831.       IF (ERR .GT. 0) CALL ERROR(17)
  832.       IF (ERR .GT. 0) RETURN
  833.       JOB = 11
  834.       IF (RHS .EQ. 2) JOB = 21
  835.       CALL WSVDC(STKR(L),STKI(L),M,M,N,STKR(LD),STKI(LD),
  836.      $        STKR(L1),STKI(L1),STKR(LU),STKI(LU),M,STKR(LV),STKI(LV),
  837.      $        N,STKR(L2),STKI(L2),JOB,ERR)
  838.       DO 51 JB = 1, N
  839.       DO 51 I = 1, K
  840.         J = N+1-JB
  841.         LL = LD+I-1+(J-1)*K
  842.         IF (I.NE.J) STKR(LL) = 0.0D0
  843.         STKI(LL) = 0.0D0
  844.         LS = LD+I-1
  845.         IF (I.EQ.J) STKR(LL) = STKR(LS)
  846.         LS = L1+I-1
  847.         IF (ERR.NE.0 .AND. I.EQ.J-1) STKR(LL) = STKR(LS)
  848.    51 CONTINUE
  849.       IF (ERR .NE. 0) CALL ERROR(24)
  850.       ERR = 0
  851.       CALL WCOPY(M*K+K*N+N*N,STKR(LU),STKI(LU),1,STKR(L),STKI(L),1)
  852.       MSTK(TOP) = M
  853.       NSTK(TOP) = K
  854.       IF (TOP+1 .GE. BOT) CALL ERROR(18)
  855.       IF (ERR .GT. 0) RETURN
  856.       TOP = TOP+1
  857.       LSTK(TOP) = L + M*K
  858.       MSTK(TOP) = K
  859.       NSTK(TOP) = N
  860.       IF (TOP+1 .GE. BOT) CALL ERROR(18)
  861.       IF (ERR .GT. 0) RETURN
  862.       TOP = TOP+1
  863.       LSTK(TOP) = L + M*K + K*N
  864.       MSTK(TOP) = N
  865.       NSTK(TOP) = N
  866.       GO TO 99
  867. C
  868.    52 LD = L + M*N
  869.       L1 = LD + MIN0(M+1,N)
  870.       L2 = L1 + N
  871.       ERR = L2+MIN0(M,N) - LSTK(BOT)
  872.       IF (ERR .GT. 0) CALL ERROR(17)
  873.       IF (ERR .GT. 0) RETURN
  874.       CALL WSVDC(STKR(L),STKI(L),M,M,N,STKR(LD),STKI(LD),
  875.      $           STKR(L1),STKI(L1),T,T,1,T,T,1,STKR(L2),STKI(L2),
  876.      $           0,ERR)
  877.       IF (ERR .NE. 0) CALL ERROR(24)
  878.       IF (ERR .GT. 0) RETURN
  879.       K = MIN0(M,N)
  880.       CALL WCOPY(K,STKR(LD),STKI(LD),1,STKR(L),STKI(L),1)
  881.       MSTK(TOP) = K
  882.       NSTK(TOP) = 1
  883.       GO TO 99
  884. C
  885. C     PINV AND RANK
  886. C
  887.    70 TOL = -1.0D0
  888.       IF (RHS .NE. 2) GO TO 71
  889.       TOL = STKR(L)
  890.       TOP = TOP-1
  891.       L = LSTK(TOP)
  892.       M = MSTK(TOP)
  893.       N = NSTK(TOP)
  894.    71 LU = L + M*N
  895.       LD = LU + M*M
  896.       IF (FIN .EQ. 5) LD = L + M*N
  897.       LV = LD + M*N
  898.       L1 = LV + N*N
  899.       IF (FIN .EQ. 5) L1 = LD + N
  900.       L2 = L1 + N
  901.       ERR = L2+MIN0(M,N) - LSTK(BOT)
  902.       IF (ERR .GT. 0) CALL ERROR(17)
  903.       IF (ERR .GT. 0) RETURN
  904.       IF (FIN .EQ. 2) JOB = 11
  905.       IF (FIN .EQ. 5) JOB = 0
  906.       CALL WSVDC(STKR(L),STKI(L),M,M,N,STKR(LD),STKI(LD),
  907.      $        STKR(L1),STKI(L1),STKR(LU),STKI(LU),M,STKR(LV),STKI(LV),
  908.      $        N,STKR(L2),STKI(L2),JOB,ERR)
  909.       IF (ERR .NE. 0) CALL ERROR(24)
  910.       IF (ERR .GT. 0) RETURN
  911.       EPS = STKR(VSIZE-4)
  912.       S = MAX0(M,N)
  913.       IF (TOL .LT. 0.0D0) TOL = FLOP(S*EPS*STKR(LD))
  914.       MN = MIN0(M,N)
  915.       K = 0
  916.       DO 72 J = 1, MN
  917.         LS = LD+J-1
  918.         S = STKR(LS)
  919.         IF (S .LE. TOL) GO TO 73
  920.         K = J
  921.         LL = LV+(J-1)*N
  922.         IF (FIN .EQ. 2) CALL WRSCAL(N,1.0D0/S,STKR(LL),STKI(LL),1)
  923.    72 CONTINUE
  924.    73 IF (FIN .EQ. 5) GO TO 78
  925.       DO 76 J = 1, M
  926.       DO 76 I = 1, N
  927.         LL = L+I-1+(J-1)*N
  928.         L1 = LV+I-1
  929.         L2 = LU+J-1
  930.         STKR(LL) = WDOTCR(K,STKR(L2),STKI(L2),M,STKR(L1),STKI(L1),N)
  931.         STKI(LL) = WDOTCI(K,STKR(L2),STKI(L2),M,STKR(L1),STKI(L1),N)
  932.    76 CONTINUE
  933.       MSTK(TOP) = N
  934.       NSTK(TOP) = M
  935.       GO TO 99
  936.    78 STKR(L) = K
  937.       STKI(L) = 0.0D0
  938.       MSTK(TOP) = 1
  939.       NSTK(TOP) = 1
  940.       GO TO 99
  941. C
  942.    99 RETURN
  943.       END
  944.