home *** CD-ROM | disk | FTP | other *** search
/ Crawly Crypt Collection 1 / crawlyvol1.bin / apps / math / matlab / src / mat2.f < prev    next >
Text File  |  1992-04-21  |  33KB  |  1,081 lines
  1.       SUBROUTINE FACTOR
  2.       DOUBLE PRECISION STKR(50005),STKI(50005)
  3.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
  4.       INTEGER IDS(4,32),PSTK(32),RSTK(32),PSIZE,PT,PTZ
  5.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  6.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  7.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
  8.       COMMON /RECU/ IDS,PSTK,RSTK,PSIZE,PT,PTZ
  9.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  10.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  11.       INTEGER SEMI,EOL,BLANK,R,ID(4),EXCNT,LPAREN,RPAREN
  12.       INTEGER STAR,DSTAR,COMMA,LESS,GREAT,QUOTE,NUM,NAME,ALFL
  13.       DATA DSTAR/54/,SEMI/39/,EOL/99/,BLANK/36/
  14.       DATA STAR/43/,COMMA/48/,LPAREN/37/,RPAREN/38/
  15.       DATA LESS/50/,GREAT/51/,QUOTE/49/,NUM/0/,NAME/1/,ALFL/52/
  16.       IF (DDT .EQ. 1) WRITE(WTE,100) PT,RSTK(PT),SYM
  17.   100 FORMAT(0X,'FACTOR',3I4)
  18.       R = RSTK(PT)
  19.       GO TO (99,99,99,99,99,99,99,01,01,25,45,65,99,99,99,55,75,32,37),R
  20.    01 IF (SYM.EQ.NUM .OR. SYM.EQ.QUOTE .OR.  SYM.EQ.LESS) GO TO 10
  21.       IF (SYM .EQ. GREAT) GO TO 30
  22.       EXCNT = 0
  23.       IF (SYM .EQ. NAME) GO TO 40
  24.       ID(1) = BLANK
  25.       IF (SYM .EQ. LPAREN) GO TO 42
  26.       CALL ERROR(2)
  27.       IF (ERR .GT. 0) RETURN
  28. C
  29. C     PUT SOMETHING ON THE STACK
  30.    10 L = 1
  31.       IF (TOP .GT. 0) L = LSTK(TOP) + MSTK(TOP)*NSTK(TOP)
  32.       IF (TOP+1 .GE. BOT) CALL ERROR(18)
  33.       IF (ERR .GT. 0) RETURN
  34.       TOP = TOP+1
  35.       LSTK(TOP) = L
  36.       IF (SYM .EQ. QUOTE) GO TO 15
  37.       IF (SYM .EQ. LESS) GO TO 20
  38. C
  39. C     SINGLE NUMBER, GETSYM STORED IT IN STKI
  40.       MSTK(TOP) = 1
  41.       NSTK(TOP) = 1
  42.       STKR(L) = STKI(VSIZE)
  43.       STKI(L) = 0.0D0
  44.       CALL GETSYM
  45.       GO TO 60
  46. C
  47. C     STRING
  48.    15 N = 0
  49.       LPT(4) = LPT(3)
  50.       CALL GETCH
  51.    16 IF (CHAR .EQ. QUOTE) GO TO 18
  52.    17 LN = L+N
  53.       IF (CHAR .EQ. EOL) CALL ERROR(31)
  54.       IF (ERR .GT. 0) RETURN
  55.       STKR(LN) = CHAR
  56.       STKI(LN) = 0.0D0
  57.       N = N+1
  58.       CALL GETCH
  59.       GO TO 16
  60.    18 CALL GETCH
  61.       IF (CHAR .EQ. QUOTE) GO TO 17
  62.       IF (N .LE. 0) CALL ERROR(31)
  63.       IF (ERR .GT. 0) RETURN
  64.       MSTK(TOP) = 1
  65.       NSTK(TOP) = N
  66.       CALL GETSYM
  67.       GO TO 60
  68. C
  69. C     EXPLICIT MATRIX
  70.    20 MSTK(TOP) = 0
  71.       NSTK(TOP) = 0
  72.    21 TOP = TOP + 1
  73.       LSTK(TOP) = LSTK(TOP-1) + MSTK(TOP-1)*NSTK(TOP-1)
  74.       MSTK(TOP) = 0
  75.       NSTK(TOP) = 0
  76.       CALL GETSYM
  77.    22 IF (SYM.EQ.SEMI .OR. SYM.EQ.GREAT .OR. SYM.EQ.EOL) GO TO 27
  78.       IF (SYM .EQ. COMMA) CALL GETSYM
  79.       PT = PT+1
  80.       RSTK(PT) = 10
  81. C     *CALL* EXPR
  82.       RETURN
  83.    25 PT = PT-1
  84.       TOP = TOP - 1
  85.       IF (MSTK(TOP) .EQ. 0) MSTK(TOP) = MSTK(TOP+1)
  86.       IF (MSTK(TOP) .NE. MSTK(TOP+1)) CALL ERROR(5)
  87.       IF (ERR .GT. 0) RETURN
  88.       NSTK(TOP) = NSTK(TOP) + NSTK(TOP+1)
  89.       GO TO 22
  90.    27 IF (SYM.EQ.SEMI .AND. CHAR.EQ.EOL) CALL GETSYM
  91.       CALL STACK1(QUOTE)
  92.       IF (ERR .GT. 0) RETURN
  93.       TOP = TOP - 1
  94.       IF (MSTK(TOP) .EQ. 0) MSTK(TOP) = MSTK(TOP+1)
  95.       IF (MSTK(TOP).NE.MSTK(TOP+1) .AND. MSTK(TOP+1).GT.0) CALL ERROR(6)
  96.       IF (ERR .GT. 0) RETURN 
  97.       NSTK(TOP) = NSTK(TOP) + NSTK(TOP+1) 
  98.       IF (SYM .EQ. EOL) CALL GETLIN
  99.       IF (SYM .NE. GREAT) GO TO 21
  100.       CALL STACK1(QUOTE)
  101.       IF (ERR .GT. 0) RETURN
  102.       CALL GETSYM
  103.       GO TO 60
  104. C
  105. C     MACRO STRING
  106.    30 CALL GETSYM
  107.       IF (SYM.EQ.LESS .AND. CHAR.EQ.EOL) CALL ERROR(28)
  108.       IF (ERR .GT. 0) RETURN
  109.       PT = PT+1
  110.       RSTK(PT) = 18
  111. C     *CALL* EXPR
  112.       RETURN
  113.    32 PT = PT-1
  114.       IF (SYM.NE.LESS .AND. SYM.NE.EOL) CALL ERROR(37)
  115.       IF (ERR .GT. 0) RETURN
  116.       IF (SYM .EQ. LESS) CALL GETSYM
  117.       K = LPT(6)
  118.       LIN(K+1) = LPT(1)
  119.       LIN(K+2) = LPT(2)
  120.       LIN(K+3) = LPT(6)
  121.       LPT(1) = K + 4
  122. C     TRANSFER STACK TO INPUT LINE
  123.       K = LPT(1)
  124.       L = LSTK(TOP)
  125.       N = MSTK(TOP)*NSTK(TOP)
  126.       DO 34 J = 1, N
  127.          LS = L + J-1
  128.          LIN(K) = IDINT(STKR(LS))
  129.          IF (LIN(K).LT.0 .OR. LIN(K).GE.ALFL) CALL ERROR(37)
  130.          IF (ERR .GT. 0) RETURN
  131.          IF (K.LT.1024) K = K+1
  132.          IF (K.EQ.1024) WRITE(WTE,33) K
  133.    33    FORMAT(0X,'INPUT BUFFER LIMIT IS ',I4,' CHARACTERS.')
  134.    34 CONTINUE
  135.       TOP = TOP-1
  136.       LIN(K) = EOL
  137.       LPT(6) = K
  138.       LPT(4) = LPT(1)
  139.       LPT(3) = 0
  140.       LPT(2) = 0
  141.       LCT(1) = 0
  142.       CHAR = BLANK
  143.       CALL GETSYM
  144.       PT = PT+1
  145.       RSTK(PT) = 19
  146. C     *CALL* EXPR
  147.       RETURN
  148.    37 PT = PT-1
  149.       K = LPT(1) - 4
  150.       LPT(1) = LIN(K+1)
  151.       LPT(4) = LIN(K+2)
  152.       LPT(6) = LIN(K+3)
  153.       CHAR = BLANK
  154.       CALL GETSYM
  155.       GO TO 60
  156. C
  157. C     FUNCTION OR MATRIX ELEMENT
  158.    40 CALL PUTID(ID,SYN)
  159.       CALL GETSYM
  160.       IF (SYM .EQ. LPAREN) GO TO 42
  161.       RHS = 0
  162.       CALL FUNS(ID)
  163.       IF (FIN .NE. 0) CALL ERROR(25)
  164.       IF (ERR .GT. 0) RETURN
  165.       CALL STACKG(ID)
  166.       IF (ERR .GT. 0) RETURN
  167.       IF (FIN .EQ. 7) GO TO 50
  168.       IF (FIN .EQ. 0) CALL PUTID(IDS(1,PT+1),ID)
  169.       IF (FIN .EQ. 0) CALL ERROR(4)
  170.       IF (ERR .GT. 0) RETURN
  171.       GO TO 60
  172. C
  173.    42 CALL GETSYM
  174.       EXCNT = EXCNT+1
  175.       PT = PT+1
  176.       PSTK(PT) = EXCNT
  177.       CALL PUTID(IDS(1,PT),ID)
  178.       RSTK(PT) = 11
  179. C     *CALL* EXPR
  180.       RETURN
  181.    45 CALL PUTID(ID,IDS(1,PT))
  182.       EXCNT = PSTK(PT)
  183.       PT = PT-1
  184.       IF (SYM .EQ. COMMA) GO TO 42
  185.       IF (SYM .NE. RPAREN) CALL ERROR(3)
  186.       IF (ERR .GT. 0) RETURN
  187.       IF (SYM .EQ. RPAREN) CALL GETSYM
  188.       IF (ID(1) .EQ. BLANK) GO TO 60
  189.       RHS = EXCNT
  190.       CALL STACKG(ID)
  191.       IF (ERR .GT. 0) RETURN
  192.       IF (FIN .EQ. 0) CALL FUNS(ID)
  193.       IF (FIN .EQ. 0) CALL ERROR(4)
  194.       IF (ERR .GT. 0) RETURN
  195. C
  196. C     EVALUATE MATRIX FUNCTION
  197.    50 PT = PT+1
  198.       RSTK(PT) = 16
  199. C     *CALL* MATFN
  200.       RETURN
  201.    55 PT = PT-1
  202.       GO TO 60
  203. C
  204. C     CHECK FOR QUOTE (TRANSPOSE) AND ** (POWER)
  205.    60 IF (SYM .NE. QUOTE) GO TO 62
  206.          I = LPT(3) - 2
  207.          IF (LIN(I) .EQ. BLANK) GO TO 90
  208.          CALL STACK1(QUOTE)
  209.          IF (ERR .GT. 0) RETURN
  210.          CALL GETSYM
  211.    62 IF (SYM.NE.STAR .OR. CHAR.NE.STAR) GO TO 90
  212.       CALL GETSYM
  213.       CALL GETSYM
  214.       PT = PT+1
  215.       RSTK(PT) = 12
  216. C     *CALL* FACTOR
  217.       GO TO 01
  218.    65 PT = PT-1
  219.       CALL STACK2(DSTAR)
  220.       IF (ERR .GT. 0) RETURN
  221.       IF (FUN .NE. 2) GO TO 90
  222. C     MATRIX POWER, USE EIGENVECTORS
  223.       PT = PT+1
  224.       RSTK(PT) = 17
  225. C     *CALL* MATFN
  226.       RETURN
  227.    75 PT = PT-1
  228.    90 RETURN
  229.    99 CALL ERROR(22)
  230.       IF (ERR .GT. 0) RETURN
  231.       RETURN
  232.       END
  233.       SUBROUTINE FUNS(ID)
  234.       INTEGER ID(4)
  235. C
  236. C     SCAN FUNCTION LIST
  237. C
  238.       DOUBLE PRECISION STKR(50005),STKI(50005)
  239.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
  240.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  241.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  242.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
  243.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  244.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  245.       LOGICAL EQID
  246.       INTEGER FUNL,FUNN(4,57),FUNP(57)
  247.       DATA FUNL/57/
  248. C
  249. C    1  ABS   ATAN  BASE  CHAR  
  250. C    2  CHOL  CHOP  COND  CONJ  
  251. C    3  COS   DET   DIAG  DIAR  
  252. C    4  DISP  EIG   EPS   EXEC  
  253. C    5  EXP   EYE   FLOP  HESS  
  254. C    6  HILB  IMAG  INV   KRON
  255. C    7  LINE  LOAD  LOG   LU 
  256. C    8  MAGIC NORM  ONES  ORTH
  257. C    9  PINV  PLOT  POLY  PRINT 
  258. C    $  PROD  QR    RAND  RANK  
  259. C    1  RAT   RCOND REAL  ROOT  
  260. C    2  ROUND RREF  SAVE  SCHUR 
  261. C    3  SIN   SIZE  SQRT  SUM
  262. C    4  SVD   TRIL  TRIU  USER  
  263. C    5  DEBUG
  264. C
  265.       DATA FUNN/
  266.      1  10,11,28,36, 10,29,10,23, 11,10,28,14, 12,17,10,27, 
  267.      2  12,17,24,21, 12,17,24,25, 12,24,23,13, 12,24,23,19, 
  268.      3  12,24,28,36, 13,14,29,36, 13,18,10,16, 13,18,10,27, 
  269.      4  13,18,28,25, 14,18,16,36, 14,25,28,36, 14,33,14,12,
  270.      5  14,33,25,36, 14,34,14,36, 15,21,24,25, 17,14,28,28, 
  271.      6  17,18,21,11, 18,22,10,16, 18,23,31,36, 20,27,24,23,
  272.      7  21,18,23,14, 21,24,10,13, 21,24,16,36, 21,30,36,36,
  273.      8  22,10,16,18, 23,24,27,22, 24,23,14,28, 24,27,29,17,
  274.      9  25,18,23,31, 25,21,24,29, 25,24,21,34, 25,27,18,23, 
  275.      $  25,27,24,13, 26,27,36,36, 27,10,23,13, 27,10,23,20, 
  276.      1  27,10,29,36, 27,12,24,23, 27,14,10,21, 27,24,24,29, 
  277.      2  27,24,30,23, 27,27,14,15, 28,10,31,14, 28,12,17,30, 
  278.      3  28,18,23,36, 28,18,35,14, 28,26,27,29, 28,30,22,36,
  279.      4  28,31,13,36, 29,27,18,21, 29,27,18,30, 30,28,14,27, 
  280.      5  13,14,11,30/
  281. C
  282.       DATA FUNP/
  283.      1  221,203,507,509, 106,609,303,225, 202,102,602,505,
  284.      4  506,211,000,501, 204,606,000,213, 105,224,101,611, 
  285.      7  508,503,206,104, 601,304,608,402, 302,510,214,504,
  286.      $  604,401,607,305, 511,103,223,215, 222,107,502,212,
  287.      3  201,610,205,603, 301,614,615,605, 512/
  288. C
  289.       IF (ID(1).EQ.0) CALL PRNTID(FUNN,FUNL-1)
  290.       IF (ID(1).EQ.0) RETURN
  291. C
  292.       DO 10 K = 1, FUNL
  293.          IF (EQID(ID,FUNN(1,K))) GO TO 20
  294.    10 CONTINUE
  295.       FIN = 0 
  296.       RETURN 
  297. C
  298.    20 FIN = MOD(FUNP(K),100)
  299.       FUN = FUNP(K)/100
  300.       IF (RHS.EQ.0 .AND. FUNP(K).EQ.606) FIN = 0
  301.       IF (RHS.EQ.0 .AND. FUNP(K).EQ.607) FIN = 0
  302.       RETURN
  303.       END
  304.       SUBROUTINE STACKP(ID)
  305.       INTEGER ID(4)
  306. C
  307. C     PUT VARIABLES INTO STORAGE
  308. C
  309.       DOUBLE PRECISION STKR(50005),STKI(50005)
  310.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
  311.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  312.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  313.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
  314.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  315.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  316.       LOGICAL EQID
  317.       INTEGER SEMI
  318.       DATA SEMI/39/
  319.       IF (DDT .EQ. 1) WRITE(WTE,100) ID
  320.   100 FORMAT(0X,'STACKP',4I4)
  321.       IF (TOP .LE. 0) CALL ERROR(1)
  322.       IF (ERR .GT. 0) RETURN
  323.       CALL FUNS(ID)
  324.       IF (FIN .NE. 0) CALL ERROR(25)
  325.       IF (ERR .GT. 0) RETURN
  326.       M = MSTK(TOP)
  327.       N = NSTK(TOP)
  328.       IF (M .GT. 0) L = LSTK(TOP)
  329.       IF (M .LT. 0) CALL ERROR(14)
  330.       IF (ERR .GT. 0) RETURN
  331.       IF (M .EQ. 0 .AND. N .NE. 0) GO TO 99
  332.       MN = M*N
  333.       LK = 0
  334.       MK = 1
  335.       NK = 0
  336.       LT = 0
  337.       MT = 0
  338.       NT = 0
  339. C
  340. C     DOES VARIABLE ALREADY EXIST
  341.       CALL PUTID(IDSTK(1,BOT-1),ID)
  342.       K = LSIZE+1
  343.    05 K = K-1
  344.       IF (.NOT.EQID(IDSTK(1,K),ID)) GO TO 05
  345.       IF (K .EQ. BOT-1) GO TO 30
  346.       LK = LSTK(K)
  347.       MK = MSTK(K)
  348.       NK = NSTK(K)
  349.       MNK = MK*NK
  350.       IF (RHS .EQ. 0) GO TO 20
  351.       IF (RHS .GT. 2) CALL ERROR(15)
  352.       IF (ERR .GT. 0) RETURN
  353.       MT = MK
  354.       NT = NK
  355.       LT = L + MN
  356.       ERR = LT + MNK - LSTK(BOT)
  357.       IF (ERR .GT. 0) CALL ERROR(17)
  358.       IF (ERR .GT. 0) RETURN
  359.       CALL WCOPY(MNK,STKR(LK),STKI(LK),1,STKR(LT),STKI(LT),1)
  360. C
  361. C     DOES IT FIT
  362.    20 IF (RHS.EQ.0 .AND. MN.EQ.MNK) GO TO 40
  363.       IF (K .GE. LSIZE-3) CALL ERROR(13)
  364.       IF (ERR .GT. 0) RETURN
  365. C
  366. C     SHIFT STORAGE
  367.       IF (K .EQ. BOT) GO TO 25
  368.       LS = LSTK(BOT)
  369.       LL = LS + MNK
  370.       CALL WCOPY(LK-LS,STKR(LS),STKI(LS),-1,STKR(LL),STKI(LL),-1)
  371.       KM1 = K-1
  372.       DO 24 IB = BOT, KM1
  373.         I = BOT+KM1-IB
  374.         CALL PUTID(IDSTK(1,I+1),IDSTK(1,I))
  375.         MSTK(I+1) = MSTK(I)
  376.         NSTK(I+1) = NSTK(I)
  377.         LSTK(I+1) = LSTK(I)+MNK
  378.    24 CONTINUE
  379. C
  380. C     DESTROY OLD VARIABLE
  381.    25 BOT = BOT+1
  382. C
  383. C     CREATE NEW VARIABLE
  384.    30 IF (MN .EQ. 0) GO TO 99
  385.       IF (BOT-2 .LE. TOP) CALL ERROR(18)
  386.       IF (ERR .GT. 0) RETURN
  387.       K = BOT-1
  388.       CALL PUTID(IDSTK(1,K), ID)
  389.       IF (RHS .EQ. 1) GO TO 50
  390.       IF (RHS .EQ. 2) GO TO 55
  391. C
  392. C     STORE
  393.    40 IF (K .LT. LSIZE) LSTK(K) = LSTK(K+1) - MN
  394.       MSTK(K) = M
  395.       NSTK(K) = N
  396.       LK = LSTK(K)
  397.       CALL WCOPY(MN,STKR(L),STKI(L),-1,STKR(LK),STKI(LK),-1)
  398.       GO TO 90
  399. C
  400. C     VECT(ARG)
  401.    50 IF (MSTK(TOP-1) .LT. 0) GO TO 59
  402.       MN1 = 1
  403.       MN2 = 1
  404.       L1 = 0
  405.       L2 = 0
  406.       IF (N.NE.1 .OR. NK.NE.1) GO TO 52
  407.       L1 = LSTK(TOP-1)
  408.       M1 = MSTK(TOP-1)
  409.       MN1 = M1*NSTK(TOP-1)
  410.       M2 = -1
  411.       GO TO 60
  412.    52 IF (M.NE.1 .OR. MK.NE.1) CALL ERROR(15)
  413.       IF (ERR .GT. 0) RETURN
  414.       L2 = LSTK(TOP-1)
  415.       M2 = MSTK(TOP-1)
  416.       MN2 = M2*NSTK(TOP-1)
  417.       M1 = -1
  418.       GO TO 60
  419. C
  420. C     MATRIX(ARG,ARG)
  421.    55 IF (MSTK(TOP-1).LT.0 .AND. MSTK(TOP-2).LT.0) GO TO 59
  422.       L2 = LSTK(TOP-1)
  423.       M2 = MSTK(TOP-1)
  424.       MN2 = M2*NSTK(TOP-1)
  425.       IF (M2 .LT. 0) MN2 = N
  426.       L1 = LSTK(TOP-2)
  427.       M1 = MSTK(TOP-2)
  428.       MN1 = M1*NSTK(TOP-2)
  429.       IF (M1 .LT. 0) MN1 = M
  430.       GO TO 60
  431. C
  432.    59 IF (MN .NE. MNK) CALL ERROR(15)
  433.       IF (ERR .GT. 0) RETURN
  434.       LK = LSTK(K)
  435.       CALL WCOPY(MN,STKR(L),STKI(L),-1,STKR(LK),STKI(LK),-1)
  436.       GO TO 90
  437. C
  438.    60 IF (MN1.NE.M .OR. MN2.NE.N) CALL ERROR(15)
  439.       IF (ERR .GT. 0) RETURN
  440.       LL = 1
  441.       IF (M1 .LT. 0) GO TO 62
  442.       DO 61 I = 1, MN1
  443.          LS = L1+I-1
  444.          MK = MAX0(MK,IDINT(STKR(LS)))
  445.          LL = MIN0(LL,IDINT(STKR(LS)))
  446.    61 CONTINUE
  447.    62 MK = MAX0(MK,M)
  448.       IF (M2 .LT. 0) GO TO 64
  449.       DO 63 I = 1, MN2
  450.          LS = L2+I-1
  451.          NK = MAX0(NK,IDINT(STKR(LS)))
  452.          LL = MIN0(LL,IDINT(STKR(LS)))
  453.    63 CONTINUE
  454.    64 NK = MAX0(NK,N)
  455.       IF (LL .LT. 1) CALL ERROR(21)
  456.       IF (ERR .GT. 0) RETURN
  457.       MNK = MK*NK
  458.       LK = LSTK(K+1) - MNK
  459.       ERR = LT + MT*NT - LK
  460.       IF (ERR .GT. 0) CALL ERROR(17)
  461.       IF (ERR .GT. 0) RETURN
  462.       LSTK(K) = LK
  463.       MSTK(K) = MK
  464.       NSTK(K) = NK
  465.       CALL WSET(MNK,0.0D0,0.0D0,STKR(LK),STKI(LK),1)
  466.       IF (NT .LT. 1) GO TO 67
  467.       DO 66 J = 1, NT
  468.          LS = LT+(J-1)*MT
  469.          LL = LK+(J-1)*MK
  470.          CALL WCOPY(MT,STKR(LS),STKI(LS),-1,STKR(LL),STKI(LL),-1)
  471.    66 CONTINUE
  472.    67 DO 68 J = 1, N
  473.       DO 68 I = 1, M
  474.         LI = L1+I-1
  475.         IF (M1 .GT. 0) LI = L1 + IDINT(STKR(LI)) - 1
  476.         LJ = L2+J-1
  477.         IF (M2 .GT. 0) LJ = L2 + IDINT(STKR(LJ)) - 1
  478.         LL = LK+LI-L1+(LJ-L2)*MK
  479.         LS = L+I-1+(J-1)*M
  480.         STKR(LL) = STKR(LS)
  481.         STKI(LL) = STKI(LS)
  482.    68 CONTINUE
  483.       GO TO 90
  484. C
  485. C     POP STACK AND PRINT IF DESIRED 
  486.    90 IF (K .EQ. BOT-1) BOT = BOT-1
  487.    99 IF (M .NE. 0) TOP = TOP - 1 - RHS
  488.       IF (M .EQ. 0) TOP = TOP - 1
  489.       IF (M .EQ. 0) RETURN
  490.       IF (SYM.NE.SEMI .AND. LCT(3).EQ.0) CALL PRINT(ID,K)
  491.       IF (SYM.EQ.SEMI .AND. LCT(3).EQ.1) CALL PRINT(ID,K)
  492.       RETURN
  493.       END
  494.       SUBROUTINE STACKG(ID)
  495.       INTEGER ID(4)
  496. C
  497. C     GET VARIABLES FROM STORAGE
  498. C
  499.       DOUBLE PRECISION STKR(50005),STKI(50005)
  500.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
  501.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  502.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  503.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
  504.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  505.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  506.       LOGICAL EQID
  507.       IF (DDT .EQ. 1) WRITE(WTE,100) ID
  508.   100 FORMAT(0X,'STACKG',4I4)
  509.       CALL PUTID(IDSTK(1,BOT-1), ID)
  510.       K = LSIZE+1
  511.    10 K = K-1
  512.       IF (.NOT.EQID(IDSTK(1,K), ID)) GO TO 10
  513.       IF (K .GE. LSIZE-1 .AND. RHS .GT. 0) GO TO 98
  514.       IF (K .EQ. BOT-1) GO TO 98
  515.       LK = LSTK(K)
  516.       IF (RHS .EQ. 1) GO TO 40
  517.       IF (RHS .EQ. 2) GO TO 60
  518.       IF (RHS .GT. 2) CALL ERROR(21)
  519.       IF (ERR .GT. 0) RETURN
  520.       L = 1
  521.       IF (TOP .GT. 0) L = LSTK(TOP) + MSTK(TOP)*NSTK(TOP)
  522.       IF (TOP+1 .GE. BOT) CALL ERROR(18)
  523.       IF (ERR .GT. 0) RETURN
  524.       TOP = TOP+1
  525. C
  526. C     LOAD VARIABLE TO TOP OF STACK
  527.       LSTK(TOP) = L
  528.       MSTK(TOP) = MSTK(K)
  529.       NSTK(TOP) = NSTK(K)
  530.       MN = MSTK(K)*NSTK(K)
  531.       ERR = L+MN - LSTK(BOT)
  532.       IF (ERR .GT. 0) CALL ERROR(17)
  533.       IF (ERR .GT. 0) RETURN
  534. C     IF RAND, MATFN6 GENERATES RANDOM NUMBER
  535.       IF (K .EQ. LSIZE) GO TO 97
  536.       CALL WCOPY(MN,STKR(LK),STKI(LK),1,STKR(L),STKI(L),1)
  537.       GO TO 99
  538. C
  539. C     VECT(ARG)
  540.    40 IF (MSTK(TOP) .EQ. 0) GO TO 99
  541.       L = LSTK(TOP)
  542.       MN = MSTK(TOP)*NSTK(TOP)
  543.       MNK = MSTK(K)*NSTK(K)
  544.       IF (MSTK(TOP) .LT. 0) MN = MNK
  545.       DO 50 I = 1, MN
  546.         LL = L+I-1
  547.         LS = LK+I-1
  548.         IF (MSTK(TOP) .GT. 0) LS = LK + IDINT(STKR(LL)) - 1
  549.         IF (LS .LT. LK .OR. LS .GE. LK+MNK) CALL ERROR(21)
  550.         IF (ERR .GT. 0) RETURN
  551.         STKR(LL) = STKR(LS)
  552.         STKI(LL) = STKI(LS)
  553.    50 CONTINUE
  554.       MSTK(TOP) = 1
  555.       NSTK(TOP) = 1
  556.       IF (MSTK(K) .GT. 1) MSTK(TOP) = MN
  557.       IF (MSTK(K) .EQ. 1) NSTK(TOP) = MN
  558.       GO TO 99
  559. C
  560. C     MATRIX(ARG,ARG)
  561.    60 TOP = TOP-1
  562.       L = LSTK(TOP)
  563.       IF (MSTK(TOP+1) .EQ. 0) MSTK(TOP) = 0
  564.       IF (MSTK(TOP) .EQ. 0) GO TO 99
  565.       L2 = LSTK(TOP+1)
  566.       M = MSTK(TOP)*NSTK(TOP)
  567.       IF (MSTK(TOP) .LT. 0) M = MSTK(K)
  568.       N = MSTK(TOP+1)*NSTK(TOP+1)
  569.       IF (MSTK(TOP+1) .LT. 0) N = NSTK(K)
  570.       L3 = L2 + N
  571.       MK = MSTK(K)
  572.       MNK = MSTK(K)*NSTK(K)
  573.       DO 70 J = 1, N
  574.       DO 70 I = 1, M
  575.         LI = L+I-1
  576.         IF (MSTK(TOP) .GT. 0) LI = L + IDINT(STKR(LI)) - 1
  577.         LJ = L2+J-1
  578.         IF (MSTK(TOP+1) .GT. 0) LJ = L2 + IDINT(STKR(LJ)) - 1
  579.         LS = LK + LI-L + (LJ-L2)*MK
  580.         IF (LS.LT.LK .OR. LS.GE.LK+MNK) CALL ERROR(21)
  581.         IF (ERR .GT. 0) RETURN
  582.         LL = L3 + I-1 + (J-1)*M
  583.         STKR(LL) = STKR(LS)
  584.         STKI(LL) = STKI(LS)
  585.    70 CONTINUE
  586.       MN = M*N
  587.       CALL WCOPY(MN,STKR(L3),STKI(L3),1,STKR(L),STKI(L),1)
  588.       MSTK(TOP) = M
  589.       NSTK(TOP) = N
  590.       GO TO 99
  591.    97 FIN = 7
  592.       FUN = 6
  593.       RETURN
  594.    98 FIN = 0
  595.       RETURN
  596.    99 FIN = -1
  597.       FUN = 0
  598.       RETURN
  599.       END
  600.       SUBROUTINE STACK1(OP)
  601.       INTEGER OP
  602. C
  603. C     UNARY OPERATIONS
  604. C
  605.       DOUBLE PRECISION STKR(50005),STKI(50005)
  606.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
  607.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  608.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  609.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
  610.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  611.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  612.       INTEGER QUOTE
  613.       DATA QUOTE/49/
  614.       IF (DDT .EQ. 1) WRITE(WTE,100) OP
  615.   100 FORMAT(0X,'STACK1',I4)
  616.       L = LSTK(TOP)
  617.       M = MSTK(TOP)
  618.       N = NSTK(TOP)
  619.       MN = M*N
  620.       IF (MN .EQ. 0) GO TO 99
  621.       IF (OP .EQ. QUOTE) GO TO 30
  622. C
  623. C     UNARY MINUS
  624.       CALL WRSCAL(MN,-1.0D0,STKR(L),STKI(L),1)
  625.       GO TO 99
  626. C
  627. C     TRANSPOSE
  628.    30 LL = L + MN
  629.       ERR = LL+MN - LSTK(BOT)
  630.       IF (ERR .GT. 0) CALL ERROR(17)
  631.       IF (ERR .GT. 0) RETURN
  632.       CALL WCOPY(MN,STKR(L),STKI(L),1,STKR(LL),STKI(LL),1)
  633.       M = NSTK(TOP)
  634.       N = MSTK(TOP)
  635.       MSTK(TOP) = M
  636.       NSTK(TOP) = N
  637.       DO 50 I = 1, M
  638.       DO 50 J = 1, N
  639.         LS = L+MN+(J-1)+(I-1)*N
  640.         LL = L+(I-1)+(J-1)*M
  641.         STKR(LL) = STKR(LS)
  642.         STKI(LL) = -STKI(LS)
  643.    50 CONTINUE
  644.       GO TO 99
  645.    99 RETURN
  646.       END
  647.       SUBROUTINE STACK2(OP)
  648.       INTEGER OP
  649. C
  650. C     BINARY AND TERNARY OPERATIONS
  651. C
  652.       DOUBLE PRECISION STKR(50005),STKI(50005)
  653.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
  654.       INTEGER IDS(4,32),PSTK(32),RSTK(32),PSIZE,PT,PTZ
  655.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  656.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  657.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
  658.       COMMON /RECU/ IDS,PSTK,RSTK,PSIZE,PT,PTZ
  659.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  660.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  661.       DOUBLE PRECISION WDOTUR,WDOTUI
  662.       DOUBLE PRECISION SR,SI,E1,ST,E2,T,FLOP
  663.       INTEGER PLUS,MINUS,STAR,DSTAR,SLASH,BSLASH,DOT,COLON
  664.       DATA PLUS/41/,MINUS/42/,STAR/43/,DSTAR/54/,SLASH/44/
  665.       DATA BSLASH/45/,DOT/47/,COLON/40/
  666. C
  667.       IF (DDT .EQ. 1) WRITE(WTE,100) OP
  668.   100 FORMAT(0X,'STACK2',I4)
  669.       L2 = LSTK(TOP)
  670.       M2 = MSTK(TOP)
  671.       N2 = NSTK(TOP)
  672.       TOP = TOP-1
  673.       L = LSTK(TOP)
  674.       M = MSTK(TOP)
  675.       N = NSTK(TOP)
  676.       FUN = 0
  677.       IF (OP .EQ. PLUS) GO TO 01
  678.       IF (OP .EQ. MINUS) GO TO 03
  679.       IF (OP .EQ. STAR) GO TO 05
  680.       IF (OP .EQ. DSTAR) GO TO 30
  681.       IF (OP .EQ. SLASH) GO TO 20
  682.       IF (OP .EQ. BSLASH) GO TO 25
  683.       IF (OP .EQ. COLON) GO TO 60
  684.       IF (OP .GT. 2*DOT) GO TO 80
  685.       IF (OP .GT. DOT) GO TO 70
  686. C
  687. C     ADDITION
  688.    01 IF (M .LT. 0) GO TO 50
  689.       IF (M2 .LT. 0) GO TO 52
  690.       IF (M .NE. M2) CALL ERROR(8)
  691.       IF (ERR .GT. 0) RETURN
  692.       IF (N .NE. N2) CALL ERROR(8)
  693.       IF (ERR .GT. 0) RETURN
  694.       CALL WAXPY(M*N,1.0D0,0.0D0,STKR(L2),STKI(L2),1,
  695.      $            STKR(L),STKI(L),1)
  696.       GO TO 99
  697. C
  698. C     SUBTRACTION
  699.    03 IF (M .LT. 0) GO TO 54
  700.       IF (M2 .LT. 0) GO TO 56
  701.       IF (M .NE. M2) CALL ERROR(9)
  702.       IF (ERR .GT. 0) RETURN
  703.       IF (N .NE. N2) CALL ERROR(9)
  704.       IF (ERR .GT. 0) RETURN
  705.       CALL WAXPY(M*N,-1.0D0,0.0D0,STKR(L2),STKI(L2),1,
  706.      $            STKR(L),STKI(L),1)
  707.       GO TO 99
  708. C
  709. C     MULTIPLICATION
  710.    05 IF (M2*M2*N2 .EQ. 1) GO TO 10
  711.       IF (M*N .EQ. 1) GO TO 11
  712.       IF (M2*N2 .EQ. 1) GO TO 10
  713.       IF (N .NE. M2) CALL ERROR(10)
  714.       IF (ERR .GT. 0) RETURN
  715.       MN = M*N2
  716.       LL = L + MN
  717.       ERR = LL+M*N+M2*N2 - LSTK(BOT)
  718.       IF (ERR .GT. 0) CALL ERROR(17)
  719.       IF (ERR .GT. 0) RETURN
  720.       CALL WCOPY(M*N+M2*N2,STKR(L),STKI(L),-1,STKR(LL),STKI(LL),-1)
  721.       DO 08 J = 1, N2
  722.       DO 08 I = 1, M
  723.         K1 = L + MN + (I-1)
  724.         K2 = L2 + MN + (J-1)*M2
  725.         K = L + (I-1) + (J-1)*M
  726.         STKR(K) = WDOTUR(N,STKR(K1),STKI(K1),M,STKR(K2),STKI(K2),1)
  727.         STKI(K) = WDOTUI(N,STKR(K1),STKI(K1),M,STKR(K2),STKI(K2),1)
  728.    08 CONTINUE
  729.       NSTK(TOP) = N2
  730.       GO TO 99
  731. C
  732. C     MULTIPLICATION BY SCALAR
  733.    10 SR = STKR(L2)
  734.       SI = STKI(L2)
  735.       L1 = L
  736.       GO TO 13
  737.    11 SR = STKR(L)
  738.       SI = STKI(L)
  739.       L1 = L+1
  740.       MSTK(TOP) = M2
  741.       NSTK(TOP) = N2
  742.    13 MN = MSTK(TOP)*NSTK(TOP)
  743.       CALL WSCAL(MN,SR,SI,STKR(L1),STKI(L1),1)
  744.       IF (L1.NE.L)
  745.      $   CALL WCOPY(MN,STKR(L1),STKI(L1),1,STKR(L),STKI(L),1)
  746.       GO TO 99
  747. C
  748. C     RIGHT DIVISION
  749.    20 IF (M2*N2 .EQ. 1) GO TO 21
  750.       IF (M2 .EQ. N2) FUN = 1
  751.       IF (M2 .NE. N2) FUN = 4
  752.       FIN = -1
  753.       RHS = 2
  754.       GO TO 99
  755.    21 SR = STKR(L2)
  756.       SI = STKI(L2)
  757.       MN = M*N
  758.       DO 22 I = 1, MN
  759.          LL = L+I-1
  760.          CALL WDIV(STKR(LL),STKI(LL),SR,SI,STKR(LL),STKI(LL))
  761.          IF (ERR .GT. 0) RETURN
  762.    22 CONTINUE
  763.       GO TO 99
  764. C
  765. C     LEFT DIVISION
  766.    25 IF (M*N .EQ. 1) GO TO 26
  767.       IF (M .EQ. N) FUN = 1
  768.       IF (M .NE. N) FUN = 4
  769.       FIN = -2
  770.       RHS = 2
  771.       GO TO 99
  772.    26 SR = STKR(L)
  773.       SI = STKI(L)
  774.       MSTK(TOP) = M2
  775.       NSTK(TOP) = N2
  776.       MN = M2*N2
  777.       DO 27 I = 1, MN
  778.          LL = L+I-1
  779.          CALL WDIV(STKR(LL+1),STKI(LL+1),SR,SI,STKR(LL),STKI(LL))
  780.          IF (ERR .GT. 0) RETURN
  781.    27 CONTINUE
  782.       GO TO 99
  783. C
  784. C     POWER
  785.    30 IF (M2*N2 .NE. 1) CALL ERROR(30)
  786.       IF (ERR .GT. 0) RETURN
  787.       IF (M .NE. N) CALL ERROR(20)
  788.       IF (ERR .GT. 0) RETURN
  789.       NEXP = IDINT(STKR(L2))
  790.       T = NEXP
  791.       IF (STKR(L2).NE.T .OR. STKI(L2).NE.0.0D0 .OR. NEXP.LT.2) GO TO 39
  792.       MN = M*N
  793.       ERR = L2+MN+N - LSTK(BOT)
  794.       IF (ERR .GT. 0) CALL ERROR(17)
  795.       IF (ERR .GT. 0) RETURN
  796.       CALL WCOPY(MN,STKR(L),STKI(L),1,STKR(L2),STKI(L2),1)
  797.       L3 = L2+MN
  798.       DO 36 KEXP = 2, NEXP
  799.         DO 35 J = 1, N
  800.           LS = L+(J-1)*N
  801.           CALL WCOPY(N,STKR(LS),STKI(LS),1,STKR(L3),STKI(L3),1)
  802.           DO 34 I = 1, N
  803.             LS = L2+I-1
  804.             LL = L+I-1+(J-1)*N
  805.             STKR(LL) = WDOTUR(N,STKR(LS),STKI(LS),N,STKR(L3),STKI(L3),1)
  806.             STKI(LL) = WDOTUI(N,STKR(LS),STKI(LS),N,STKR(L3),STKI(L3),1)
  807.    34     CONTINUE
  808.    35   CONTINUE
  809.    36 CONTINUE
  810.       GO TO 99
  811. C
  812. C     NONINTEGER OR NONPOSITIVE POWER, USE EIGENVECTORS
  813.    39 FUN = 2
  814.       FIN = 0
  815.       GO TO 99
  816. C
  817. C     ADD OR SUBTRACT SCALAR
  818.    50 IF (M2 .NE. N2) CALL ERROR(8)
  819.       IF (ERR .GT. 0) RETURN
  820.       M = M2
  821.       N = N2
  822.       MSTK(TOP) = M
  823.       NSTK(TOP) = N
  824.       SR = STKR(L)
  825.       SI = STKI(L)
  826.       CALL WCOPY(M*N,STKR(L+1),STKI(L+1),1,STKR(L),STKI(L),1)
  827.       GO TO 58
  828.    52 IF (M .NE. N) CALL ERROR(8)
  829.       IF (ERR .GT. 0) RETURN
  830.       SR = STKR(L2)
  831.       SI = STKI(L2)
  832.       GO TO 58
  833.    54 IF (M2 .NE. N2) CALL ERROR(9)
  834.       IF (ERR .GT. 0) RETURN
  835.       M = M2
  836.       N = N2
  837.       MSTK(TOP) = M
  838.       NSTK(TOP) = N
  839.       SR = STKR(L)
  840.       SI = STKI(L)
  841.       CALL WCOPY(M*N,STKR(L+1),STKI(L+1),1,STKR(L),STKI(L),1)
  842.       CALL WRSCAL(M*N,-1.0D0,STKR(L),STKI(L),1)
  843.       GO TO 58
  844.    56 IF (M .NE. N) CALL ERROR(9)
  845.       IF (ERR .GT. 0) RETURN
  846.       SR = -STKR(L2)
  847.       SI = -STKI(L2)
  848.       GO TO 58
  849.    58 DO 59 I = 1, N
  850.          LL = L + (I-1)*(N+1)
  851.          STKR(LL) = FLOP(STKR(LL)+SR)
  852.          STKI(LL) = FLOP(STKI(LL)+SI)
  853.    59 CONTINUE
  854.       GO TO 99
  855. C
  856. C     COLON
  857.    60 E2 = STKR(L2)
  858.       ST = 1.0D0
  859.       N = 0
  860.       IF (RHS .LT. 3) GO TO 61
  861.       ST = STKR(L)
  862.       TOP = TOP-1
  863.       L = LSTK(TOP)
  864.       IF (ST .EQ. 0.0D0) GO TO 63
  865.    61 E1 = STKR(L)
  866. C     CHECK FOR CLAUSE
  867.       IF (RSTK(PT) .EQ. 3) GO TO 64
  868.       ERR = L + MAX0(3,IDINT((E2-E1)/ST)) - LSTK(BOT)
  869.       IF (ERR .GT. 0) CALL ERROR(17)
  870.       IF (ERR .GT. 0) RETURN
  871.    62 IF (ST .GT. 0.0D0 .AND. STKR(L) .GT. E2) GO TO 63
  872.       IF (ST .LT. 0.0D0 .AND. STKR(L) .LT. E2) GO TO 63
  873.         N = N+1
  874.         L = L+1
  875.         T = N
  876.         STKR(L) = E1 + T*ST
  877.         STKI(L) = 0.0D0
  878.         GO TO 62
  879.    63 NSTK(TOP) = N
  880.       MSTK(TOP) = 1
  881.       IF (N .EQ. 0) MSTK(TOP) = 0
  882.       GO TO 99
  883. C
  884. C     FOR CLAUSE 
  885.    64 STKR(L) = E1
  886.       STKR(L+1) = ST
  887.       STKR(L+2) = E2
  888.       MSTK(TOP) = -3
  889.       NSTK(TOP) = -1
  890.       GO TO 99
  891. C
  892. C     ELEMENTWISE OPERATIONS
  893.    70 OP = OP - DOT
  894.       IF (M.NE.M2 .OR. N.NE.N2) CALL ERROR(10)
  895.       IF (ERR .GT. 0) RETURN
  896.       MN = M*N
  897.       DO 72 I = 1, MN
  898.          J = L+I-1
  899.          K = L2+I-1
  900.          IF (OP .EQ. STAR)
  901.      $      CALL WMUL(STKR(J),STKI(J),STKR(K),STKI(K),STKR(J),STKI(J))
  902.          IF (OP .EQ. SLASH)
  903.      $      CALL WDIV(STKR(J),STKI(J),STKR(K),STKI(K),STKR(J),STKI(J))
  904.          IF (OP .EQ. BSLASH)
  905.      $      CALL WDIV(STKR(K),STKI(K),STKR(J),STKI(J),STKR(J),STKI(J))
  906.          IF (ERR .GT. 0) RETURN
  907.    72 CONTINUE
  908.       GO TO 99
  909. C
  910. C     KRONECKER
  911.    80 FIN = OP - 2*DOT - STAR + 11
  912.       FUN = 6
  913.       TOP = TOP + 1
  914.       RHS = 2
  915.       GO TO 99
  916. C
  917.    99 RETURN
  918.       END
  919.       SUBROUTINE PRINT(ID,K)
  920. C     PRIMARY OUTPUT ROUTINE
  921.       INTEGER ID(4),K
  922.       DOUBLE PRECISION STKR(50005),STKI(50005)
  923.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
  924.       INTEGER ALFA(52),ALFB(52),ALFL,CASE
  925.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  926.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  927.       CHARACTER CLS
  928.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
  929.       COMMON /ALFS/ ALFA,ALFB,ALFL,CASE
  930.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  931.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  932.       DOUBLE PRECISION S,TR,TI,PR(12),PI(12)
  933.       INTEGER FNO(11),FNL(11),SIG(12),PLUS,MINUS,BLANK,TYP,F
  934.       CHARACTER MYCHAR,CSIG(12)
  935.       DATA PLUS/41/,MINUS/42/,BLANK/36/
  936. C     FORMAT NUMBERS AND LENGTHS
  937.       DATA FNO /11,12,21,22,23,24,31,32,33,34,-1/
  938.       DATA FNL /12, 6, 8, 4, 6, 3, 4, 2, 3, 1, 1/
  939. C     FMT   1       2       3       4       5
  940. C         SHORT   LONG   SHORT E  LONG E    Z
  941. C     TYP   1       2       3
  942. C         INTEGER  REAL   COMPLEX
  943.       IF (DDT .EQ. 1) WRITE(WTE,100) ID
  944.   100 FORMAT(0X,'PRINT ',4I4)
  945.       IF (LCT(1) .LT. 0) GO TO 99
  946.       L = LSTK(K)
  947.       M = MSTK(K)
  948.       N = NSTK(K)
  949.       MN = M*N
  950.       TYP = 1
  951.       S = 0.0D0
  952.       DO 10 I = 1, MN
  953.         LS = L+I-1
  954.         TR = STKR(LS)
  955.         TI = STKI(LS)
  956.         KTR = DMIN1(DABS(TR),2147483647D0)
  957.         S = DMAX1(S,DABS(TR),DABS(TI))
  958.         IF (KTR .NE. DABS(TR)) TYP = MAX0(2,TYP)
  959.         IF (TI .NE. 0.0D0) TYP = 3
  960.    10 CONTINUE
  961.       IF (S .NE. 0.0D0) S = DLOG10(S)
  962.       KS = IDINT(S)
  963.       IF (-2 .LE. KS .AND. KS .LE. 1) KS = 0
  964.       IF (KS .EQ. 2 .AND. FMT .EQ. 1 .AND. TYP .EQ. 2) KS = 0
  965.       IF (TYP .EQ. 1 .AND. KS .LE. 2) F = 1
  966.       IF (TYP .EQ. 1 .AND. KS .GT. 2) F = 2
  967.       IF (TYP .EQ. 1 .AND. KS .GT. 9) TYP = 2
  968.       IF (TYP .EQ. 2) F = FMT + 2
  969.       IF (TYP .EQ. 3) F = FMT + 6
  970.       IF (MN.EQ.1 .AND. KS.NE.0 .AND. FMT.LT.3 .AND. TYP.NE.1) F = F+2
  971.       IF (FMT .EQ. 5) F = 11
  972.       JINC = FNL(F)
  973.       F = FNO(F)
  974.       S = 1.0D0
  975.       IF (F.EQ.21 .OR. F.EQ.22 .OR. F.EQ.31 .OR. F.EQ.32) S = 10.0D0**KS
  976.       LS = ((N-1)/JINC+1)*M + 2
  977.       IF (LCT(1) + LS .LE. LCT(2)) GO TO 20
  978.          LCT(1) = 0
  979.          WRITE(WTE,43) LS
  980.          READ(RTE,44,END=19) CLS
  981.          LS=ICHAR(CLS)
  982.          IF (LS .EQ. ALFA(BLANK+1)) GO TO 20
  983.          LCT(1) = -1
  984.          GO TO 99
  985.    19    CALL FILES(-RTE,BUF,BUF)
  986.    20 CONTINUE
  987.       WRITE(WTE,44)
  988.       IF (WIO .NE. 0) WRITE(WIO,44)
  989.       CALL PRNTID(ID,-1)
  990.       LCT(1) = LCT(1)+2
  991.       LUNIT = WTE
  992.    50 IF (S .NE. 1.0D0) WRITE(LUNIT,41) S
  993.       DO 80 J1 = 1, N, JINC
  994.         J2 = MIN0(N, J1+JINC-1)
  995.         WRITE(LUNIT,44)
  996.         IF (N .GT. JINC) WRITE(LUNIT,42) J1,J2
  997.         DO 70 I = 1, M
  998.           JM = J2-J1+1
  999.           DO 60 J = 1, JM
  1000.              LS = L+I-1+(J+J1-2)*M
  1001.              IF (F .LT. 20) BUF(J) = STKR(LS)
  1002.              IF (F .LT. 20) GO TO 60
  1003.              PR(J) = STKR(LS)/S
  1004.              PI(J) = DABS(STKI(LS)/S)
  1005.              SIG(J) = ALFA(PLUS+1)
  1006.              IF (STKI(LS) .LT. 0.0D0) SIG(J) = ALFA(MINUS+1)
  1007.    60     CONTINUE
  1008.           DO 1000 J=1,JM
  1009.              CSIG(J)=MYCHAR(SIG(J))
  1010.  1000     CONTINUE
  1011.           IF (F .EQ. 11) WRITE(LUNIT,11)(BUF(J),J=1,JM)
  1012.           IF (F .EQ. 12) WRITE(LUNIT,12)(BUF(J),J=1,JM)
  1013.           IF (F .EQ. 21) WRITE(LUNIT,21)(PR(J),J=1,JM)
  1014.           IF (F .EQ. 22) WRITE(LUNIT,22)(PR(J),J=1,JM)
  1015.           IF (F .EQ. 23) WRITE(LUNIT,23)(PR(J),J=1,JM)
  1016.           IF (F .EQ. 24) WRITE(LUNIT,24)(PR(J),J=1,JM)
  1017.           IF (F .EQ. 31) WRITE(LUNIT,31)(PR(J),CSIG(J),PI(J),J=1,JM)
  1018.           IF (F .EQ. 32) WRITE(LUNIT,32)(PR(J),CSIG(J),PI(J),J=1,JM)
  1019.           IF (F .EQ. 33) WRITE(LUNIT,33)(PR(J),CSIG(J),PI(J),J=1,JM)
  1020.           IF (F .EQ. 34) WRITE(LUNIT,34)(PR(J),CSIG(J),PI(J),J=1,JM)
  1021.           IF (F .EQ. -1) CALL FORMZ(LUNIT,STKR(LS),STKI(LS))
  1022.           LCT(1) = LCT(1)+1
  1023.    70   CONTINUE
  1024.    80 CONTINUE
  1025.       IF (LUNIT.EQ.WIO .OR. WIO.EQ.0) GO TO 99
  1026.       LUNIT = WIO
  1027.       GO TO 50
  1028.    99 RETURN
  1029. C
  1030.    11 FORMAT(0X,12I6)
  1031.    12 FORMAT(0X,6I12)
  1032.    21 FORMAT(0X,F9.4,7F10.4)
  1033.    22 FORMAT(0X,F19.15,3F20.15)
  1034.    23 FORMAT(0X,1P6D13.4)
  1035.    24 FORMAT(0X,1P3D24.15)
  1036.    31 FORMAT(0X,4(F9.4,' ',A1,F7.4,'i'))
  1037.    32 FORMAT(0X,F19.15,A1,F18.15,'i',F20.15,A1,F18.15,'i')
  1038.    33 FORMAT(0X,3(1PD13.4,' ',A1,1PD10.4,'i'))
  1039.    34 FORMAT(0X,1PD24.15,' ',A1,1PD21.15,'i')
  1040.    41 FORMAT(/0X,' ',1PD9.1,2H *)
  1041.    42 FORMAT(0X,'    COLUMNS',I3,' THRU',I3)
  1042.    43 FORMAT(/0X,'AT LEAST ',I5,' MORE LINES.',
  1043.      $       '  ENTER BLANK LINE TO CONTINUE OUTPUT.')
  1044.    44 FORMAT(A1)
  1045. C
  1046.       END
  1047.       SUBROUTINE PRNTID(ID,ARGCNT)
  1048. C     PRINT VARIABLE NAMES
  1049.       INTEGER ID(4,1),ARGCNT
  1050.       INTEGER ALFA(52),ALFB(52),ALFL,CASE
  1051.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  1052.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  1053.       COMMON /ALFS/ ALFA,ALFB,ALFL,CASE
  1054.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  1055.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  1056.       CHARACTER CBUF(256)
  1057.       CHARACTER MYCHAR
  1058.       INTEGER EQUAL
  1059.       DATA EQUAL/46/
  1060.       J1 = 1
  1061.    10 J2 = MIN0(J1+7,IABS(ARGCNT))
  1062.       L = 0
  1063.       DO 15 J = J1,J2
  1064.       DO 15 I = 1, 4
  1065.       K = ID(I,J)+1
  1066.       L = L+1
  1067.       BUF(L) = ALFA(K)
  1068.    15 CONTINUE
  1069.       IF (ARGCNT .EQ. -1) L=L+1
  1070.       IF (ARGCNT .EQ. -1) BUF(L) = ALFA(EQUAL+1)
  1071.       DO 1000 I=1,L
  1072.         CBUF(I)=MYCHAR(BUF(I))
  1073.  1000 CONTINUE
  1074.       WRITE(WTE,20) (CBUF(I),I=1,L)
  1075.       IF (WIO .NE. 0) WRITE(WIO,20) (CBUF(I),I=1,L)
  1076.    20 FORMAT(0X,8(4A1,2H  ))
  1077.       J1 = J1+8
  1078.       IF (J1 .LE. IABS(ARGCNT)) GO TO 10
  1079.       RETURN
  1080.       END
  1081.