home *** CD-ROM | disk | FTP | other *** search
/ Crawly Crypt Collection 1 / crawlyvol1.bin / apps / math / matlab / src / mat4.f < prev    next >
Text File  |  1992-04-21  |  29KB  |  996 lines
  1.       SUBROUTINE MATFN4
  2. C
  3. C     EVALUATE FUNCTIONS INVOLVING QR DECOMPOSITION (LEAST SQUARES)
  4. C
  5.       DOUBLE PRECISION STKR(50005),STKI(50005)
  6.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
  7.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  8.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  9.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
  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 S,T,TOL,EPS,FLOP
  13.       INTEGER QUOTE
  14.       DATA QUOTE/49/
  15. C
  16.       IF (DDT .EQ. 1) WRITE(WTE,100) FIN
  17.   100 FORMAT(0X,'MATFN4',I4)
  18. C
  19.       L = LSTK(TOP)
  20.       M = MSTK(TOP)
  21.       N = NSTK(TOP)
  22.       IF (FIN .EQ. -1) GO TO 10
  23.       IF (FIN .EQ. -2) GO TO 20
  24.       GO TO 40
  25. C
  26. C     RECTANGULAR MATRIX RIGHT DIVISION, A/A2
  27.    10 L2 = LSTK(TOP+1)
  28.       M2 = MSTK(TOP+1)
  29.       N2 = NSTK(TOP+1)
  30.       TOP = TOP + 1
  31.       IF (N.GT.1 .AND. N.NE.N2) CALL ERROR(11)
  32.       IF (ERR .GT. 0) RETURN
  33.       CALL STACK1(QUOTE)
  34.       IF (ERR .GT. 0) RETURN
  35.       LL = L2+M2*N2
  36.       CALL WCOPY(M*N,STKR(L),STKI(L),1,STKR(LL),STKI(LL),1)
  37.       CALL WCOPY(M*N+M2*N2,STKR(L2),STKI(L2),1,STKR(L),STKI(L),1)
  38.       LSTK(TOP) = L+M2*N2
  39.       MSTK(TOP) = M
  40.       NSTK(TOP) = N
  41.       CALL STACK1(QUOTE)
  42.       IF (ERR .GT. 0) RETURN
  43.       TOP = TOP - 1
  44.       M = N2
  45.       N = M2
  46.       GO TO 20
  47. C
  48. C     RECTANGULAR MATRIX LEFT DIVISION A BACKSLASH A2
  49. C
  50.    20 L2 = LSTK(TOP+1)
  51.       M2 = MSTK(TOP+1)
  52.       N2 = NSTK(TOP+1)
  53.       IF (M2*N2 .GT. 1) GO TO 21
  54.         M2 = M
  55.         N2 = M
  56.         ERR = L2+M*M - LSTK(BOT)
  57.         IF (ERR .GT. 0) CALL ERROR(17)
  58.         IF (ERR .GT. 0) RETURN
  59.         CALL WSET(M*M-1,0.0D0,0.0D0,STKR(L2+1),STKI(L2+1),1)
  60.         CALL WCOPY(M,STKR(L2),STKI(L2),0,STKR(L2),STKI(L2),M+1)
  61.    21 IF (M2 .NE. M) CALL ERROR(12)
  62.       IF (ERR .GT. 0) RETURN
  63.       L3 = L2 + MAX0(M,N)*N2
  64.       L4 = L3 + N
  65.       ERR = L4 + N - LSTK(BOT)
  66.       IF (ERR .GT. 0) CALL ERROR(17)
  67.       IF (ERR .GT. 0) RETURN
  68.       IF (M .GT. N) GO TO 23
  69.       DO 22 JB = 1, N2
  70.         J = N+1-JB
  71.         LS = L2 + (J-1)*M
  72.         LL = L2 + (J-1)*N
  73.         CALL WCOPY(M,STKR(LS),STKI(LS),-1,STKR(LL),STKI(LL),-1)
  74.    22 CONTINUE
  75.    23 DO 24 J = 1, N
  76.         BUF(J) = 0
  77.    24 CONTINUE
  78.       CALL WQRDC(STKR(L),STKI(L),M,M,N,STKR(L4),STKI(L4),
  79.      $           BUF,STKR(L3),STKI(L3),1)
  80.       K = 0
  81.       EPS = STKR(VSIZE-4) 
  82.       S = MAX0(M,N)
  83.       T = DABS(STKR(L))+DABS(STKI(L))
  84.       TOL = FLOP(EPS*S*T)
  85.       MN = MIN0(M,N)
  86.       DO 27 J = 1, MN
  87.         LS = L+J-1+(J-1)*M
  88.         T = DABS(STKR(LS)) + DABS(STKI(LS))
  89.         IF (T .GT. TOL) K = J
  90.    27 CONTINUE
  91.       IF (K .LT. MN) WRITE(WTE,28) K,TOL
  92.       IF (K.LT.MN .AND. WIO.NE.0) WRITE(WIO,28) K,TOL
  93.    28 FORMAT(0X,'RANK DEFICIENT,  RANK =',I4,',  TOL =',1PD13.4)
  94.       MN = MAX0(M,N)
  95.       DO 29 J = 1, N2
  96.         LS = L2+(J-1)*MN
  97.         CALL WQRSL(STKR(L),STKI(L),M,M,K,STKR(L4),STKI(L4),
  98.      $             STKR(LS),STKI(LS),T,T,STKR(LS),STKI(LS),
  99.      $             STKR(LS),STKI(LS),T,T,T,T,100,INFO)
  100.         LL = LS+K
  101.         CALL WSET(N-K,0.0D0,0.0D0,STKR(LL),STKI(LL),1)
  102.    29 CONTINUE
  103.       DO 31 J = 1, N
  104.         BUF(J) = -BUF(J)
  105.    31 CONTINUE
  106.       DO 35 J = 1, N
  107.         IF (BUF(J) .GT. 0) GO TO 35
  108.         K = -BUF(J)
  109.         BUF(J) = K
  110.    33   CONTINUE
  111.           IF (K .EQ. J) GO TO 34
  112.           LS = L2+J-1
  113.           LL = L2+K-1
  114.           CALL WSWAP(N2,STKR(LS),STKI(LS),MN,STKR(LL),STKI(LL),MN)
  115.           BUF(K) = -BUF(K)
  116.           K = BUF(K)
  117.           GO TO 33
  118.    34   CONTINUE
  119.    35 CONTINUE
  120.       DO 36 J = 1, N2
  121.         LS = L2+(J-1)*MN
  122.         LL = L+(J-1)*N
  123.         CALL WCOPY(N,STKR(LS),STKI(LS),1,STKR(LL),STKI(LL),1)
  124.    36 CONTINUE
  125.       MSTK(TOP) = N
  126.       NSTK(TOP) = N2
  127.       IF (FIN .EQ. -1) CALL STACK1(QUOTE)
  128.       IF (ERR .GT. 0) RETURN
  129.       GO TO 99
  130. C
  131. C     QR
  132. C
  133.    40 MM = MAX0(M,N)
  134.       LS = L + MM*MM
  135.       IF (LHS.EQ.1 .AND. FIN.EQ.1) LS = L
  136.       LE = LS + M*N
  137.       L4 = LE + MM
  138.       ERR = L4+MM - LSTK(BOT)
  139.       IF (ERR .GT. 0) CALL ERROR(17)
  140.       IF (ERR .GT. 0) RETURN
  141.       IF (LS.NE.L) CALL WCOPY(M*N,STKR(L),STKI(L),1,STKR(LS),STKI(LS),1)
  142.       JOB = 1
  143.       IF (LHS.LT.3) JOB = 0
  144.       DO 42 J = 1, N
  145.         BUF(J) = 0
  146.    42 CONTINUE
  147.       CALL WQRDC(STKR(LS),STKI(LS),M,M,N,STKR(L4),STKI(L4),
  148.      $            BUF,STKR(LE),STKI(LE),JOB)
  149.       IF (LHS.EQ.1 .AND. FIN.EQ.1) GO TO 99
  150.       CALL WSET(M*M,0.0D0,0.0D0,STKR(L),STKI(L),1)
  151.       CALL WSET(M,1.0D0,0.0D0,STKR(L),STKI(L),M+1)
  152.       DO 43 J = 1, M
  153.         LL = L+(J-1)*M
  154.         CALL WQRSL(STKR(LS),STKI(LS),M,M,N,STKR(L4),STKI(L4),
  155.      $             STKR(LL),STKI(LL),STKR(LL),STKI(LL),T,T,
  156.      $             T,T,T,T,T,T,10000,INFO)
  157.    43 CONTINUE
  158.       IF (FIN .EQ. 2) GO TO 99
  159.       NSTK(TOP) = M
  160.       DO 45 J = 1, N
  161.         LL = LS+J+(J-1)*M
  162.         CALL WSET(M-J,0.0D0,0.0D0,STKR(LL),STKI(LL),1)
  163.    45 CONTINUE
  164.       IF (TOP+1 .GE. BOT) CALL ERROR(18)
  165.       IF (ERR .GT. 0) RETURN
  166.       TOP = TOP+1
  167.       LSTK(TOP) = LS
  168.       MSTK(TOP) = M
  169.       NSTK(TOP) = N
  170.       IF (LHS .EQ. 2) GO TO 99
  171.       CALL WSET(N*N,0.0D0,0.0D0,STKR(LE),STKI(LE),1)
  172.       DO 47 J = 1, N
  173.         LL = LE+BUF(J)-1+(J-1)*N
  174.         STKR(LL) = 1.0D0
  175.    47 CONTINUE
  176.       IF (TOP+1 .GE. BOT) CALL ERROR(18)
  177.       IF (ERR .GT. 0) RETURN
  178.       TOP = TOP+1
  179.       LSTK(TOP) = LE
  180.       MSTK(TOP) = N
  181.       NSTK(TOP) = N
  182.       GO TO 99
  183. C
  184.    99 RETURN
  185.       END
  186.       SUBROUTINE MATFN5
  187. C
  188. C     FILE HANDLING AND OTHER I/O
  189. C
  190.       DOUBLE PRECISION STKR(50005),STKI(50005)
  191.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
  192.       INTEGER ALFA(52),ALFB(52),ALFL,CASE
  193.       INTEGER IDS(4,32),PSTK(32),RSTK(32),PSIZE,PT,PTZ
  194.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  195.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  196.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
  197.       COMMON /ALFS/ ALFA,ALFB,ALFL,CASE
  198.       COMMON /RECU/ IDS,PSTK,RSTK,PSIZE,PT,PTZ
  199.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  200.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  201.       INTEGER EOL,CH,BLANK,FLAG,TOP2,PLUS,MINUS,QUOTE,SEMI,LRAT,MRAT
  202.       CHARACTER CCH,MYCHAR,CBUF(256)
  203.       INTEGER ID(4)
  204.       DOUBLE PRECISION EPS,B,S,T,FLOP,WASUM
  205.       LOGICAL TEXT
  206.       DATA EOL/99/,BLANK/36/,PLUS/41/,MINUS/42/,QUOTE/49/,SEMI/39/
  207.       DATA LRAT/5/,MRAT/100/
  208. C
  209.       IF (DDT .EQ. 1) WRITE(WTE,100) FIN
  210.   100 FORMAT(0X,'MATFN5',I4)
  211. C     FUNCTIONS/FIN
  212. C     EXEC SAVE LOAD PRIN DIAR DISP BASE LINE CHAR PLOT RAT  DEBU
  213. C      1    2    3    4    5    6    7    8    9   10   11   12
  214.       L = LSTK(TOP)
  215.       M = MSTK(TOP)
  216.       N = NSTK(TOP)
  217.       IF (FIN .GT. 5) GO TO 15
  218. C
  219. C     CONVERT FILE NAME
  220.       MN = M*N
  221.       FLAG = 3
  222.       IF (SYM .EQ. SEMI) FLAG = 0
  223.       IF (RHS .LT. 2) GO TO 12
  224.          FLAG = IDINT(STKR(L))
  225.          TOP2 = TOP
  226.          TOP = TOP-1
  227.          L = LSTK(TOP)
  228.          MN = MSTK(TOP)*NSTK(TOP)
  229.    12 LUN = -1
  230.       IF (MN.EQ.1 .AND. STKR(L).LT.10.0D0) LUN = IDINT(STKR(L))
  231.       IF (LUN .EQ. RTE .AND. FIN .NE. 1) CALL ERROR(7)
  232.       IF (LUN .EQ. WTE .OR. LUN .EQ. HIO) CALL ERROR(7)
  233.       IF (ERR .GT. 0) RETURN
  234.       IF (LUN .GE. 0) GO TO 15
  235.       DO 14 J = 1, 64
  236.          LS = L+J-1
  237.          IF (J .LE. MN) CH = IDINT(STKR(LS))
  238.          IF (J .GT. MN) CH = BLANK
  239.          IF (CH.LT.0 .OR. CH.GE.ALFL) CALL ERROR(38)
  240.          IF (ERR .GT. 0) RETURN
  241.          IF (CASE .EQ. 0) BUF(J) = ALFA(CH+1)
  242.          IF (CASE .EQ. 1) BUF(J) = ALFB(CH+1)
  243.    14 CONTINUE
  244.       IOSTAT = 0
  245.       GO TO 15
  246. C
  247.    15 GO TO (20,30,35,25,27,60,65,70,50,80,40,95),FIN
  248. C
  249. C     EXEC
  250.    20 IF (LUN .EQ. 0) GO TO 23
  251.       K = LPT(6)
  252.       LIN(K+1) = LPT(1)
  253.       LIN(K+2) = LPT(3)
  254.       LIN(K+3) = LPT(6)
  255.       LIN(K+4) = PTZ
  256.       LIN(K+5) = RIO
  257.       LIN(K+6) = LCT(4)
  258.       LPT(1) = K + 7
  259.       LCT(4) = FLAG
  260.       PTZ = PT - 4
  261.       IF (RIO .EQ. RTE) RIO = 10
  262.       RIO = RIO + 1
  263.       IF (LUN .GT. 0) RIO = LUN
  264.       IF (LUN .LT. 0) CALL FILES(RIO,BUF,IOSTAT)
  265.       IF (IOSTAT .NE. 0) CALL ERROR(7)
  266.       IF (ERR .GT. 0) GO TO 23
  267.       IF (FLAG .GE. 4) WRITE(WTE,22)
  268.    22 FORMAT(0X,'PAUSE MODE. ENTER BLANK LINES.')
  269.       SYM = EOL
  270.       MSTK(TOP) = 0
  271.       GO TO 99
  272. C
  273. C     EXEC(0)
  274.    23 RIO = RTE
  275.       ERR = 99
  276.       GO TO 99
  277. C
  278. C     PRINT
  279.    25 IF (LUN .LT. 0) CALL FILES(7,BUF,IOSTAT)
  280.       IF (IOSTAT .NE. 0) CALL ERROR(7)
  281.       IF (ERR .GT. 0) RETURN
  282.       K = WTE
  283.       IF (LUN .LT. 0) WTE = 7
  284.       IF (LUN .GT. 0) WTE = LUN
  285.       L = LCT(2)
  286.       LCT(2) = 9999
  287.       IF (RHS .GT. 1) CALL PRINT(SYN,TOP2)
  288.       IF (LUN .LT. 0) CALL FILES(-WTE,BUF,IOSTAT)
  289.       WTE = K
  290.       LCT(2) = L
  291.       MSTK(TOP) = 0
  292.       GO TO 99
  293. C
  294. C     DIARY     
  295.    27 IF (WIO .NE. 0 .AND. WIO .NE. LUN) CALL FILES(-WIO,BUF,IOSTAT)
  296.       IF (WIO .NE. LUN) WIO = 0
  297.       IF (LUN .LT. 0) CALL FILES(8,BUF,IOSTAT)
  298.       IF (IOSTAT .NE. 0) CALL ERROR(7)
  299.       IF (ERR .GT. 0) RETURN
  300.       IF (LUN .LT. 0) WIO = 8
  301.       IF (LUN .GT. 0) WIO = LUN
  302.       MSTK(TOP) = 0
  303.       GO TO 99
  304. C
  305. C     SAVE
  306.    30 IF (LUN .LT. 0) LUNIT = 1
  307.       IF (LUN .LT. 0) CALL FILES(LUNIT,BUF,IOSTAT)
  308.       IF (IOSTAT .NE. 0) CALL ERROR(7)
  309.       IF (ERR .GT. 0) RETURN
  310.       IF (LUN .GT. 0) LUNIT = LUN
  311.       K = LSIZE-4
  312.       IF (K .LT. BOT) K = LSIZE
  313.       IF (RHS .EQ. 2) K = TOP2
  314.       IF (RHS .EQ. 2) CALL PUTID(IDSTK(1,K),SYN)
  315.    32 L = LSTK(K)
  316.       M = MSTK(K)
  317.       N = NSTK(K)
  318.       DO 34 I = 1, 4
  319.          J = IDSTK(I,K)+1
  320.          BUF(I) = ALFA(J)
  321.    34 CONTINUE
  322.       IMG = 0
  323.       IF (WASUM(M*N,STKI(L),STKI(L),1) .NE. 0.0D0) IMG = 1
  324.       CALL SAVLOD(LUNIT,BUF,M,N,IMG,0,STKR(L),STKI(L))
  325.       K = K-1
  326.       IF (K .GE. BOT) GO TO 32
  327.       CALL FILES(-LUNIT,BUF,IOSTAT)
  328.       MSTK(TOP) = 0
  329.       GO TO 99
  330. C
  331. C     LOAD
  332.    35 IF (LUN .LT. 0) LUNIT = 2
  333.       IF (LUN .LT. 0) CALL FILES(LUNIT,BUF,IOSTAT)
  334.       IF (IOSTAT .NE. 0) CALL ERROR(7)
  335.       IF (ERR .GT. 0) RETURN
  336.       IF (LUN .GT. 0) LUNIT = LUN
  337.    36 JOB = LSTK(BOT) - L
  338.       CALL SAVLOD(LUNIT,ID,MSTK(TOP),NSTK(TOP),IMG,JOB,STKR(L),STKI(L))
  339.       MN = MSTK(TOP)*NSTK(TOP)
  340.       IF (MN .EQ. 0) GO TO 39
  341.       IF (IMG .EQ. 0) CALL RSET(MN,0.0D0,STKI(L),1)
  342.       DO 38 I = 1, 4
  343.          J = 0
  344.    37    J = J+1
  345.          IF (ID(I).NE.ALFA(J) .AND. J.LE.BLANK) GO TO 37
  346.          ID(I) = J-1
  347.    38 CONTINUE
  348.       SYM = SEMI
  349.       RHS = 0
  350.       CALL STACKP(ID)
  351.       TOP = TOP + 1
  352.       GO TO 36
  353.    39 CALL FILES(-LUNIT,BUF,IOSTAT)
  354.       MSTK(TOP) = 0
  355.       GO TO 99
  356. C
  357. C     RAT
  358.    40 IF (RHS .EQ. 2) GO TO 44
  359.       MN = M*N
  360.       L2 = L
  361.       IF (LHS .EQ. 2) L2 = L + MN
  362.       LW = L2 + MN
  363.       ERR = LW + LRAT - LSTK(BOT)
  364.       IF (ERR .GT. 0) CALL ERROR(17)
  365.       IF (ERR .GT. 0) RETURN
  366.       IF (LHS .EQ. 2) TOP = TOP + 1
  367.       LSTK(TOP) = L2
  368.       MSTK(TOP) = M
  369.       NSTK(TOP) = N
  370.       CALL RSET(LHS*MN,0.0D0,STKI(L),1)
  371.       DO 42 I = 1, MN
  372.          CALL RAT(STKR(L),LRAT,MRAT,S,T,STKR(LW))
  373.          STKR(L) = S
  374.          STKR(L2) = T 
  375.          IF (LHS .EQ. 1) STKR(L) = FLOP(S/T)
  376.          L = L + 1
  377.          L2 = L2 + 1
  378.    42 CONTINUE
  379.       GO TO 99
  380.    44 MRAT = IDINT(STKR(L))
  381.       LRAT = IDINT(STKR(L-1))
  382.       TOP = TOP - 1
  383.       MSTK(TOP) = 0
  384.       GO TO 99
  385. C
  386. C     CHAR
  387.    50 K = IABS(IDINT(STKR(L)))
  388.       IF (M*N.NE.1 .OR. K.GE.ALFL) CALL ERROR(36)
  389.       IF (ERR .GT. 0) RETURN
  390.       CH = ALFA(K+1)
  391.       IF (STKR(L) .LT. 0.0D0) CH = ALFB(K+1)
  392.       CCH=MYCHAR(CH)
  393.       WRITE(WTE,51) CCH
  394.    51 FORMAT(0X,'REPLACE CHARACTER ',A1)
  395.       READ(RTE,52) CCH
  396.       CH=ICHAR(CCH)
  397.    52 FORMAT(A1)
  398.       IF (STKR(L) .GE. 0.0D0) ALFA(K+1) = CH
  399.       IF (STKR(L) .LT. 0.0D0) ALFB(K+1) = CH
  400.       MSTK(TOP) = 0
  401.       GO TO 99
  402. C
  403. C     DISP
  404.    60 WRITE(WTE,61)
  405.       IF (WIO .NE. 0) WRITE(WIO,61)
  406.    61 FORMAT(0X,80A1)
  407.       IF (RHS .EQ. 2) GO TO 65
  408.       MN = M*N
  409.       TEXT = .TRUE.
  410.       DO 62 I = 1, MN
  411.         LS = L+I-1
  412.         CH = IDINT(STKR(LS))
  413.         T = CH
  414.         TEXT = TEXT .AND. (CH.GE.0) .AND. (CH.LT.ALFL)
  415.         TEXT = TEXT .AND. (T.EQ.STKR(LS))
  416.    62 CONTINUE
  417.       DO 64 I = 1, M
  418.       DO 63 J = 1, N
  419.         LS = L+I-1+(J-1)*M
  420.         IF (STKR(LS) .EQ. 0.0D0) CH = BLANK
  421.         IF (STKR(LS) .GT. 0.0D0) CH = PLUS
  422.         IF (STKR(LS) .LT. 0.0D0) CH = MINUS
  423.         IF (TEXT) CH = IDINT(STKR(LS))
  424.         BUF(J) = ALFA(CH+1)
  425.    63 CONTINUE
  426.       DO 1010 J=1,N
  427.         CBUF(J)=MYCHAR(BUF(J))
  428.  1010 CONTINUE
  429.       WRITE(WTE,61) (CBUF(J),J=1,N)
  430.       IF (WIO .NE. 0) WRITE(WIO,61) (CBUF(J),J=1,N)
  431.    64 CONTINUE
  432.       MSTK(TOP) = 0
  433.       GO TO 99
  434. C
  435. C     BASE
  436.    65 IF (RHS .NE. 2) CALL ERROR(39)
  437.       IF (STKR(L) .LE. 1.0D0) CALL ERROR(36)
  438.       IF (ERR .GT. 0) RETURN
  439.       B = STKR(L)
  440.       L2 = L
  441.       TOP = TOP-1
  442.       RHS = 1
  443.       L = LSTK(TOP)
  444.       M = MSTK(TOP)*NSTK(TOP)
  445.       EPS = STKR(VSIZE-4)
  446.       DO 66 I = 1, M
  447.          LS = L2+(I-1)*N
  448.          LL = L+I-1
  449.          CALL BASE(STKR(LL),B,EPS,STKR(LS),N)
  450.    66 CONTINUE
  451.       CALL RSET(M*N,0.0D0,STKI(L2),1)
  452.       CALL WCOPY(M*N,STKR(L2),STKI(L2),1,STKR(L),STKI(L),1)
  453.       MSTK(TOP) = N
  454.       NSTK(TOP) = M
  455.       CALL STACK1(QUOTE)
  456.       IF (FIN .EQ. 6) GO TO 60
  457.       GO TO 99
  458. C
  459. C     LINES
  460.    70 LCT(2) = IDINT(STKR(L))
  461.       MSTK(TOP) = 0
  462.       GO TO 99
  463. C
  464. C     PLOT
  465.    80 IF (RHS .GE. 2) GO TO 82
  466.       N = M*N
  467.       DO 81 I = 1, N
  468.          LL = L+I-1
  469.          STKI(LL) = I
  470.    81 CONTINUE
  471.       CALL PLOT(WTE,STKI(L),STKR(L),N,T,0,BUF)
  472.       IF (WIO .NE. 0) CALL PLOT(WIO,STKI(L),STKR(L),N,T,0,BUF)
  473.       MSTK(TOP) = 0
  474.       GO TO 99
  475.    82 IF (RHS .EQ. 2) K = 0
  476.       IF (RHS .EQ. 3) K = M*N
  477.       IF (RHS .GT. 3) K = RHS - 2
  478.       TOP = TOP - (RHS - 1)
  479.       N = MSTK(TOP)*NSTK(TOP)
  480.       IF (MSTK(TOP+1)*NSTK(TOP+1) .NE. N) CALL ERROR(5)
  481.       IF (ERR .GT. 0) RETURN
  482.       LX = LSTK(TOP)
  483.       LY = LSTK(TOP+1)
  484.       IF (RHS .GT. 3) L = LSTK(TOP+2)
  485.       CALL PLOT(WTE,STKR(LX),STKR(LY),N,STKR(L),K,BUF)
  486.       IF (WIO .NE. 0) CALL PLOT(WIO,STKR(LX),STKR(LY),N,STKR(L),K,BUF)
  487.       MSTK(TOP) = 0
  488.       GO TO 99
  489. C
  490. C     DEBUG
  491.    95 DDT = IDINT(STKR(L))
  492.       WRITE(WTE,96) DDT
  493.    96 FORMAT(0X,'DEBUG ',I4)
  494.       MSTK(TOP) = 0
  495.       GO TO 99
  496. C
  497.    99 RETURN
  498.       END
  499.       SUBROUTINE MATFN6
  500. C
  501. C     EVALUATE UTILITY FUNCTIONS
  502. C
  503.       DOUBLE PRECISION STKR(50005),STKI(50005)
  504.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
  505.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  506.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  507.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
  508.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  509.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  510.       INTEGER SEMI,ID(4),UNIFOR(4),NORMAL(4),SEED(4)
  511.       DOUBLE PRECISION EPS0,EPS,S,SR,SI,T
  512.       DOUBLE PRECISION FLOP,URAND
  513.       LOGICAL EQID
  514.       DATA SEMI/39/
  515.       DATA UNIFOR/30,23,18,15/,NORMAL/23,24,27,22/,SEED/28,14,14,13/
  516. C
  517.       IF (DDT .EQ. 1) WRITE(WTE,100) FIN
  518.   100 FORMAT(0X,'MATFN6',I4)
  519. C     FUNCTIONS/FIN
  520. C     MAGI DIAG SUM  PROD USER EYE  RAND ONES CHOP SIZE KRON  TRIL TRIU
  521. C       1    2    3    4    5    6    7    8    9   10  11-13  14   15
  522.       L = LSTK(TOP)
  523.       M = MSTK(TOP)
  524.       N = NSTK(TOP)
  525.       GO TO (75,80,65,67,70,90,90,90,60,77,50,50,50,80,80),FIN
  526. C
  527. C     KRONECKER PRODUCT
  528.    50 IF (RHS .NE. 2) CALL ERROR(39)
  529.       IF (ERR .GT. 0) RETURN
  530.       TOP = TOP - 1
  531.       L = LSTK(TOP)
  532.       MA = MSTK(TOP)
  533.       NA = NSTK(TOP)
  534.       LA = L + MAX0(M*N*MA*NA,M*N+MA*NA)
  535.       LB = LA + MA*NA
  536.       ERR = LB + M*N - LSTK(BOT)
  537.       IF (ERR .GT. 0) CALL ERROR(17)
  538.       IF (ERR .GT. 0) RETURN
  539. C     MOVE A AND B ABOVE RESULT
  540.       CALL WCOPY(MA*NA+M*N,STKR(L),STKI(L),1,STKR(LA),STKI(LA),1)
  541.       DO 54 JA = 1, NA
  542.         DO 53 J = 1, N
  543.           LJ = LB + (J-1)*M
  544.           DO 52 IA = 1, MA
  545. C           GET J-TH COLUMN OF B
  546.             CALL WCOPY(M,STKR(LJ),STKI(LJ),1,STKR(L),STKI(L),1)
  547. C           ADDRESS OF A(IA,JA)
  548.             LS = LA + IA-1 + (JA-1)*MA
  549.             DO 51 I = 1, M
  550. C             A(IA,JA) OP B(I,J)
  551.               IF (FIN .EQ. 11) CALL WMUL(STKR(LS),STKI(LS),
  552.      $           STKR(L),STKI(L),STKR(L),STKI(L))
  553.               IF (FIN .EQ. 12) CALL WDIV(STKR(LS),STKI(LS),
  554.      $           STKR(L),STKI(L),STKR(L),STKI(L))
  555.               IF (FIN .EQ. 13) CALL WDIV(STKR(L),STKI(L),
  556.      $           STKR(LS),STKI(LS),STKR(L),STKI(L))
  557.               IF (ERR .GT. 0) RETURN
  558.               L = L + 1
  559.    51       CONTINUE
  560.    52     CONTINUE
  561.    53   CONTINUE
  562.    54 CONTINUE
  563.       MSTK(TOP) = M*MA
  564.       NSTK(TOP) = N*NA
  565.       GO TO 99
  566. C
  567. C     CHOP
  568.    60 EPS0 = 1.0D0
  569.    61 EPS0 = EPS0/2.0D0
  570.       T = FLOP(1.0D0 + EPS0)
  571.       IF (T .GT. 1.0D0) GO TO 61
  572.       EPS0 = 2.0D0*EPS0
  573.       FLP(2) = IDINT(STKR(L))
  574.       IF (SYM .NE. SEMI) WRITE(WTE,62) FLP(2)
  575.    62 FORMAT(/0X,'CHOP ',I2,' PLACES.')
  576.       EPS = 1.0D0
  577.    63 EPS = EPS/2.0D0
  578.       T = FLOP(1.0D0 + EPS)
  579.       IF (T .GT. 1.0D0) GO TO 63
  580.       EPS = 2.0D0*EPS
  581.       T = STKR(VSIZE-4)
  582.       IF (T.LT.EPS .OR. T.EQ.EPS0) STKR(VSIZE-4) = EPS
  583.       MSTK(TOP) = 0
  584.       GO TO 99
  585. C
  586. C     SUM
  587.    65 SR = 0.0D0
  588.       SI = 0.0D0
  589.       MN = M*N
  590.       DO 66 I = 1, MN
  591.          LS = L+I-1
  592.          SR = FLOP(SR+STKR(LS))
  593.          SI = FLOP(SI+STKI(LS))
  594.    66 CONTINUE
  595.       GO TO 69
  596. C
  597. C     PROD
  598.    67 SR = 1.0D0
  599.       SI = 0.0D0
  600.       MN = M*N
  601.       DO 68 I = 1, MN
  602.          LS = L+I-1
  603.          CALL WMUL(STKR(LS),STKI(LS),SR,SI,SR,SI)
  604.    68 CONTINUE
  605.    69 STKR(L) = SR
  606.       STKI(L) = SI
  607.       MSTK(TOP) = 1
  608.       NSTK(TOP) = 1
  609.       GO TO 99
  610. C
  611. C     USER
  612.    70 S = 0.0D0
  613.       T = 0.0D0
  614.       IF (RHS .LT. 2) GO TO 72
  615.       IF (RHS .LT. 3) GO TO 71
  616.       T = STKR(L)
  617.       TOP = TOP-1
  618.       L = LSTK(TOP)
  619.       M = MSTK(TOP)
  620.       N = NSTK(TOP)
  621.    71 S = STKR(L)
  622.       TOP = TOP-1
  623.       L = LSTK(TOP)
  624.       M = MSTK(TOP)
  625.       N = NSTK(TOP)
  626.    72 CALL USER(STKR(L),M,N,S,T)
  627.       CALL RSET(M*N,0.0D0,STKI(L),1)
  628.       MSTK(TOP) = M
  629.       NSTK(TOP) = N
  630.       GO TO 99
  631. C
  632. C     MAGIC
  633.    75 N = MAX0(IDINT(STKR(L)),0)
  634.       IF (N .EQ. 2) N = 0
  635.       IF (N .GT. 0) CALL MAGIC(STKR(L),N,N)
  636.       CALL RSET(N*N,0.0D0,STKI(L),1)
  637.       MSTK(TOP) = N
  638.       NSTK(TOP) = N
  639.       GO TO 99
  640. C
  641. C     SIZE
  642.    77 STKR(L) = M
  643.       STKR(L+1) = N
  644.       STKI(L) = 0.0D0
  645.       STKI(L+1) = 0.0D0
  646.       MSTK(TOP) = 1
  647.       NSTK(TOP) = 2
  648.       IF (LHS .EQ. 1) GO TO 99
  649.       NSTK(TOP) = 1
  650.       TOP = TOP + 1
  651.       LSTK(TOP) = L+1
  652.       MSTK(TOP) = 1
  653.       NSTK(TOP) = 1
  654.       GO TO 99
  655. C
  656. C     DIAG, TRIU, TRIL
  657.    80 K = 0
  658.       IF (RHS .NE. 2) GO TO 81
  659.          K = IDINT(STKR(L))
  660.          TOP = TOP-1
  661.          L = LSTK(TOP)
  662.          M = MSTK(TOP)
  663.          N = NSTK(TOP)
  664.    81 IF (FIN .GE. 14) GO TO 85
  665.       IF (M .EQ. 1 .OR. N .EQ. 1) GO TO 83
  666.       IF (K.GE.0) MN=MIN0(M,N-K)
  667.       IF (K.LT.0) MN=MIN0(M+K,N)
  668.       MSTK(TOP) = MAX0(MN,0)
  669.       NSTK(TOP) = 1
  670.       IF (MN .LE. 0) GO TO 99
  671.       DO 82 I = 1, MN
  672.          IF (K.GE.0) LS = L+(I-1)+(I+K-1)*M
  673.          IF (K.LT.0) LS = L+(I-K-1)+(I-1)*M
  674.          LL = L+I-1
  675.          STKR(LL) = STKR(LS)
  676.          STKI(LL) = STKI(LS)
  677.    82 CONTINUE
  678.       GO TO 99
  679.    83 N = MAX0(M,N)+IABS(K)
  680.       ERR = L+N*N - LSTK(BOT)
  681.       IF (ERR .GT. 0) CALL ERROR(17)
  682.       IF (ERR .GT. 0) RETURN
  683.       MSTK(TOP) = N
  684.       NSTK(TOP) = N
  685.       DO 84 JB = 1, N
  686.       DO 84 IB = 1, N
  687.          J = N+1-JB
  688.          I = N+1-IB
  689.          SR = 0.0D0
  690.          SI = 0.0D0
  691.          IF (K.GE.0) LS = L+I-1
  692.          IF (K.LT.0) LS = L+J-1
  693.          LL = L+I-1+(J-1)*N
  694.          IF (J-I .EQ. K) SR = STKR(LS)
  695.          IF (J-I .EQ. K) SI = STKI(LS)
  696.          STKR(LL) = SR
  697.          STKI(LL) = SI
  698.    84 CONTINUE
  699.       GO TO 99
  700. C
  701. C     TRIL, TRIU
  702.    85 DO 87 J = 1, N
  703.          LD = L + J - K - 1 + (J-1)*M
  704.          IF (FIN .EQ. 14) LL = J - K - 1
  705.          IF (FIN .EQ. 14) LS = LD - LL
  706.          IF (FIN .EQ. 15) LL = M - J + K
  707.          IF (FIN .EQ. 15) LS = LD + 1
  708.          IF (LL .GT. 0) CALL WSET(LL,0.0D0,0.0D0,STKR(LS),STKI(LS),1)
  709.    87 CONTINUE
  710.       GO TO 99
  711. C
  712. C     EYE, RAND, ONES
  713.    90 IF (M.GT.1 .OR. RHS.EQ.0) GO TO 94
  714.       IF (RHS .NE. 2) GO TO 91
  715.         NN = IDINT(STKR(L))
  716.         TOP = TOP-1
  717.         L = LSTK(TOP)
  718.         N = NSTK(TOP)
  719.    91 IF (FIN.NE.7 .OR. N.LT.4) GO TO 93
  720.       DO 92 I = 1, 4
  721.         LS = L+I-1
  722.         ID(I) = IDINT(STKR(LS))
  723.    92 CONTINUE
  724.       IF (EQID(ID,UNIFOR).OR.EQID(ID,NORMAL)) GO TO 97
  725.       IF (EQID(ID,SEED)) GO TO 98
  726.    93 IF (N .GT. 1) GO TO 94
  727.       M = MAX0(IDINT(STKR(L)),0)
  728.       IF (RHS .EQ. 2) N = MAX0(NN,0)
  729.       IF (RHS .NE. 2) N = M
  730.       ERR = L+M*N - LSTK(BOT)
  731.       IF (ERR .GT. 0) CALL ERROR(17)
  732.       IF (ERR .GT. 0) RETURN
  733.       MSTK(TOP) = M
  734.       NSTK(TOP) = N
  735.       IF (M*N .EQ. 0) GO TO 99
  736.    94 DO 96 J = 1, N
  737.       DO 96 I = 1, M
  738.         LL = L+I-1+(J-1)*M
  739.         STKR(LL) = 0.0D0
  740.         STKI(LL) = 0.0D0
  741.         IF (I.EQ.J .OR. FIN.EQ.8) STKR(LL) = 1.0D0
  742.         IF (FIN.EQ.7 .AND. RAN(2).EQ.0) STKR(LL) = FLOP(URAND(RAN(1)))
  743.         IF (FIN.NE.7 .OR. RAN(2).EQ.0) GO TO 96
  744.    95      SR = 2.0D0*URAND(RAN(1))-1.0D0
  745.            SI = 2.0D0*URAND(RAN(1))-1.0D0
  746.            T = SR*SR + SI*SI
  747.            IF (T .GT. 1.0D0) GO TO 95
  748.         STKR(LL) = FLOP(SR*DSQRT(-2.0D0*DLOG(T)/T))
  749.    96 CONTINUE
  750.       GO TO 99
  751. C
  752. C     SWITCH UNIFORM AND NORMAL
  753.    97 RAN(2) = ID(1) - UNIFOR(1)
  754.       MSTK(TOP) = 0
  755.       GO TO 99
  756. C     
  757. C     SEED
  758.    98 IF (RHS .EQ. 2) RAN(1) = NN
  759.       STKR(L) = RAN(1)
  760.       MSTK(TOP) = 1
  761.       IF (RHS .EQ. 2) MSTK(TOP) = 0
  762.       NSTK(TOP) = 1
  763.       GO TO 99
  764. C
  765.    99 RETURN
  766.       END
  767.       SUBROUTINE ERROR(N)
  768.       INTEGER N
  769.       DOUBLE PRECISION STKR(50005),STKI(50005)
  770.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP
  771.       INTEGER IDS(4,32),PSTK(32),RSTK(32),PSIZE,PT,PTZ
  772.       INTEGER ALFA(52),ALFB(52),ALFL,CASE
  773.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  774.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)
  775.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP
  776.       COMMON /RECU/ IDS,PSTK,RSTK,PSIZE,PT,PTZ
  777.       COMMON /ALFS/ ALFA,ALFB,ALFL,CASE
  778.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  779.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN
  780.       CHARACTER CBUF(256),MYCHAR
  781.       CHARACTER ERRMSG(8),BLH,BEL
  782.       DATA ERRMSG /'/','-','-','E','R','R','O','R'/,BLH/' '/,BEL/'◆'/
  783. C     SET BEL TO CTRL-G IF POSSIBLE
  784. C
  785.       K = LPT(2) - LPT(1)
  786.       IF (K .LT. 1) K = 1
  787.       LUNIT = WTE
  788.    98 WRITE(LUNIT,100) (BLH,I=1,K),(ERRMSG(I),I=1,8),BEL
  789.   100 FORMAT(0X,80A1)
  790.       GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,
  791.      $      23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40),N
  792. C
  793.     1 WRITE(LUNIT,101)
  794.   101 FORMAT(0X,'IMPROPER MULTIPLE ASSIGNMENT')
  795.       GO TO 99
  796.     2 WRITE(LUNIT,102)
  797.   102 FORMAT(0X,'IMPROPER FACTOR')
  798.       GO TO 99
  799.     3 WRITE(LUNIT,103)
  800.   103 FORMAT(0X,'EXPECT RIGHT PARENTHESIS')
  801.       GO TO 99
  802.     4 DO 94 I = 1, 4
  803.          K = IDS(I,PT+1)
  804.          BUF(I) = ALFA(K+1)
  805.      CBUF(I)=MYCHAR(BUF(I))
  806.    94 CONTINUE
  807.       WRITE(LUNIT,104) (CBUF(I),I=1,4)
  808.   104 FORMAT(0X,'UNDEFINED VARIABLE: ',4A1)
  809.       GO TO 99
  810.     5 WRITE(LUNIT,105)
  811.   105 FORMAT(0X,'COLUMN LENGTHS DO NOT MATCH')
  812.       GO TO 99
  813.     6 WRITE(LUNIT,106)
  814.   106 FORMAT(0X,'ROW LENGTHS DO NOT MATCH')
  815.       GO TO 99
  816.     7 K = 65
  817.    77 K = K-1
  818.       IF (BUF(K) .EQ. ICHAR(BLH)) GO TO 77
  819.       DO 1075 I=1,K
  820.         CBUF(I)=MYCHAR(BUF(I))
  821.  1075 CONTINUE
  822.       WRITE(LUNIT,107) (CBUF(I),I=1,K)
  823.   107 FORMAT(0X,'CAN NOT OPEN FILE: ',64A1)
  824.       GO TO 99
  825.     8 WRITE(LUNIT,108)
  826.   108 FORMAT(0X,'INCOMPATIBLE FOR ADDITION')
  827.       GO TO 99
  828.     9 WRITE(LUNIT,109)
  829.   109 FORMAT(0X,'INCOMPATIBLE FOR SUBTRACTION')
  830.       GO TO 99
  831.    10 WRITE(LUNIT,110)
  832.   110 FORMAT(0X,'INCOMPATIBLE FOR MULTIPLICATION')
  833.       GO TO 99
  834.    11 WRITE(LUNIT,111)
  835.   111 FORMAT(0X,'INCOMPATIBLE FOR RIGHT DIVISION')
  836.       GO TO 99
  837.    12 WRITE(LUNIT,112)
  838.   112 FORMAT(0X,'INCOMPATIBLE FOR LEFT DIVISION')
  839.       GO TO 99
  840.    13 WRITE(LUNIT,113)
  841.   113 FORMAT(0X,'IMPROPER ASSIGNMENT TO PERMANENT VARIABLE')
  842.       GO TO 99
  843.    14 WRITE(LUNIT,114)
  844.   114 FORMAT(0X,'EYE-DENTITY UNDEFINED BY CONTEXT')
  845.       GO TO 99
  846.    15 WRITE(LUNIT,115)
  847.   115 FORMAT(0X,'IMPROPER ASSIGNMENT TO SUBMATRIX')
  848.       GO TO 99
  849.    16 WRITE(LUNIT,116)
  850.   116 FORMAT(0X,'IMPROPER COMMAND')
  851.       GO TO 99
  852.    17 LB = VSIZE - LSTK(BOT) + 1
  853.       LT = ERR + LSTK(BOT)
  854.       WRITE(LUNIT,117) LB,LT,VSIZE
  855.   117 FORMAT(0X,'TOO MUCH MEMORY REQUIRED'
  856.      $  /0X,' ',I7,' VARIABLES,',I7,' TEMPORARIES,',I7,' AVAILABLE.')
  857.       GO TO 99
  858.    18 WRITE(LUNIT,118)
  859.   118 FORMAT(0X,'TOO MANY NAMES')
  860.       GO TO 99
  861.    19 WRITE(LUNIT,119)
  862.   119 FORMAT(0X,'MATRIX IS SINGULAR TO WORKING PRECISION')
  863.       GO TO 99
  864.    20 WRITE(LUNIT,120)
  865.   120 FORMAT(0X,'MATRIX MUST BE SQUARE')
  866.       GO TO 99
  867.    21 WRITE(LUNIT,121)
  868.   121 FORMAT(0X,'SUBSCRIPT OUT OF RANGE')
  869.       GO TO 99
  870.    22 WRITE(LUNIT,122) (RSTK(I),I=1,PT)
  871.   122 FORMAT(0X,'RECURSION DIFFICULTIES',10I4)
  872.       GO TO 99
  873.    23 WRITE(LUNIT,123)
  874.   123 FORMAT(0X,'ONLY 1, 2 OR INF NORM OF MATRIX')
  875.       GO TO 99
  876.    24 WRITE(LUNIT,124)
  877.   124 FORMAT(0X,'NO CONVERGENCE')
  878.       GO TO 99
  879.    25 WRITE(LUNIT,125)
  880.   125 FORMAT(0X,'CAN NOT USE FUNCTION NAME AS VARIABLE')
  881.       GO TO 99
  882.    26 WRITE(LUNIT,126)
  883.   126 FORMAT(0X,'TOO COMPLICATED (STACK OVERFLOW)')
  884.       GO TO 99
  885.    27 WRITE(LUNIT,127)
  886.   127 FORMAT(0X,'DIVISION BY ZERO IS A NO-NO')
  887.       GO TO 99
  888.    28 WRITE(LUNIT,128)
  889.   128 FORMAT(0X,'EMPTY MACRO')
  890.       GO TO 99
  891.    29 WRITE(LUNIT,129)
  892.   129 FORMAT(0X,'NOT POSITIVE DEFINITE')
  893.       GO TO 99
  894.    30 WRITE(LUNIT,130)
  895.   130 FORMAT(0X,'IMPROPER EXPONENT')
  896.       GO TO 99
  897.    31 WRITE(LUNIT,131)
  898.   131 FORMAT(0X,'IMPROPER STRING')
  899.       GO TO 99
  900.    32 WRITE(LUNIT,132)
  901.   132 FORMAT(0X,'SINGULARITY OF LOG OR ATAN')
  902.       GO TO 99
  903.    33 WRITE(LUNIT,133)
  904.   133 FORMAT(0X,'TOO MANY COLONS')
  905.       GO TO 99
  906.    34 WRITE(LUNIT,134)
  907.   134 FORMAT(0X,'IMPROPER FOR CLAUSE')
  908.       GO TO 99
  909.    35 WRITE(LUNIT,135)
  910.   135 FORMAT(0X,'IMPROPER WHILE OR IF CLAUSE')
  911.       GO TO 99
  912.    36 WRITE(LUNIT,136)
  913.   136 FORMAT(0X,'ARGUMENT OUT OF RANGE')
  914.       GO TO 99
  915.    37 WRITE(LUNIT,137)
  916.   137 FORMAT(0X,'IMPROPER MACRO')
  917.       GO TO 99
  918.    38 WRITE(LUNIT,138)
  919.   138 FORMAT(0X,'IMPROPER FILE NAME')
  920.       GO TO 99
  921.    39 WRITE(LUNIT,139)
  922.   139 FORMAT(0X,'INCORRECT NUMBER OF ARGUMENTS')
  923.       GO TO 99
  924.    40 WRITE(LUNIT,140)
  925.   140 FORMAT(0X,'EXPECT STATEMENT TERMINATOR')
  926.       GO TO 99
  927. C
  928.    99 ERR = N
  929.       IF (LUNIT.EQ.WIO .OR. WIO.EQ.0) RETURN
  930.       LUNIT = WIO
  931.       GO TO 98
  932.       END
  933.       DOUBLE PRECISION FUNCTION PYTHAG(A,B)
  934.       DOUBLE PRECISION A,B
  935. C
  936. C     FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW
  937. C
  938.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  939.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  940.       DOUBLE PRECISION P,R,S,T,U
  941.       P = DMAX1(DABS(A),DABS(B))
  942.       IF (P .EQ. 0.0D0) GO TO 20
  943.       R = (DMIN1(DABS(A),DABS(B))/P)**2
  944.       IF (DDT .EQ. 25) WRITE(WTE,1)
  945.       IF (DDT .EQ. 25) WRITE(WTE,2) P,R
  946.     1 FORMAT(0X,'PYTHAG',1P2D23.15)
  947.     2 FORMAT(0X,1P2D23.15)
  948.    10 CONTINUE
  949.          T = 4.0D0 + R
  950.          IF (T .EQ. 4.0D0) GO TO 20
  951.          S = R/T
  952.          U = 1.0D0 + 2.0D0*S
  953.          P = U*P
  954.          R = (S/U)**2 * R
  955.          IF (DDT .EQ. 25) WRITE(WTE,2) P,R
  956.       GO TO 10
  957.    20 PYTHAG = P
  958.       RETURN
  959.       END
  960.       SUBROUTINE RAT(X,LEN,MAXD,A,B,D)
  961.       INTEGER LEN,MAXD
  962.       DOUBLE PRECISION X,A,B,D(LEN)
  963. C
  964. C     A/B = CONTINUED FRACTION APPROXIMATION TO X
  965. C           USING  LEN  TERMS EACH LESS THAN MAXD
  966. C 
  967.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO
  968.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO
  969.       DOUBLE PRECISION S,T,Z,ROUND 
  970.       Z = X
  971.       S = MAXD
  972.       DO 10 I = 1, LEN
  973.          K = I
  974.          D(K) = ROUND(Z)  
  975.          Z = Z - D(K)
  976.          IF (S*DABS(Z) .LE. 1.0D0) GO TO 20
  977.          Z = 1.0D0/Z
  978.    10 CONTINUE
  979.    20 T = D(K)
  980.       S = 1.0D0
  981.       IF (K .LT. 2) GO TO 40
  982.       DO 30 IB = 2, K
  983.          I = K+1-IB
  984.          Z = T
  985.          T = D(I)*T + S
  986.          S = Z
  987.    30 CONTINUE
  988.    40 IF (S .LT. 0.0D0) T = -T
  989.       IF (S .LT. 0.0D0) S = -S
  990.       IF (DDT .EQ. 27) WRITE(WTE,50) X,T,S,(D(I),I=1,K)
  991.    50 FORMAT(/0X,1PD23.15,0PF8.0,' /',F8.0,4X,6F5.0/(0X,45X,6F5.0))
  992.       A = T
  993.       B = S
  994.       RETURN
  995.       END
  996.