home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / fortran / linpkdrv.arc / ST.FOR < prev    next >
Text File  |  1984-01-12  |  16KB  |  546 lines
  1. C     MAIN PROGRAM
  2.       INTEGER LUNIT
  3. C     ALLOW 5000 UNDERFLOWS.
  4. C     CALL TRAPS(0,0,5001,0,0)
  5. C
  6. C     OUTPUT UNIT NUMBER
  7. C
  8.       LUNIT = 6
  9. C
  10.       CALL STRTS(LUNIT)
  11.       STOP
  12.       END
  13.       SUBROUTINE STRTS(LUNIT)
  14.       INTEGER LUNIT
  15. C     LUNIT IS THE OUTPUT UNIT NUMBER
  16. C
  17. C     TESTS
  18. C        STRCO, STRSL
  19. C
  20. C     LINPACK. THIS VERSION DATED 08/14/78 .
  21. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  22. C
  23. C     SUBROUTINES AND FUNCTIONS
  24. C
  25. C     LINPACK STRCO,STRSL
  26. C     EXTERNAL STRXX,SMACH
  27. C     BLAS SAXPY,SDOT,SSCAL,SASUM
  28. C     FORTRAN ABS,AMAX1,FLOAT,MAX0
  29. C
  30. C     INTERNAL VARIABLES
  31. C
  32.       REAL A(15,15),B(15),BT(15),X(15),XEXACT(15),XT(15),Z(15)
  33.       REAL AINV(15,15),DET(2)
  34.       REAL SDOT,STUFF,T
  35.       REAL ANORM,AINORM,RCOND,COND,COND1,SMACH,EPS
  36.       REAL ENORM,ETNORM,RNORM,RTNORM,XNORM,XTNORM,EN,SASUM
  37.       REAL FNORM,ONEPX,Q(7),QS(7)
  38.       INTEGER I,INFO,J,JOB,KASE,KOUNT,KSING,LDA,ML,MU,N,NPRINT
  39.       INTEGER KSUSP(7),IQ(7)
  40. C
  41.       LDA = 15
  42. C
  43. C     WRITE MATRIX AND SOLUTIONS IF  N .LE. NPRINT
  44. C
  45.       NPRINT = 3
  46. C
  47.       WRITE (LUNIT,380)
  48.       WRITE (LUNIT,730)
  49. C
  50.       DO 10 I = 1, 7
  51.          KSUSP(I) = 0
  52.    10 CONTINUE
  53.       KSING = 0
  54. C
  55. C     SET EPS TO ROUNDING UNIT
  56. C
  57.       EPS = SMACH(1)
  58.       WRITE (LUNIT,390) EPS
  59.       WRITE (LUNIT,370)
  60. C
  61. C     START MAIN LOOP
  62. C
  63.       KASE = 1
  64.    20 CONTINUE
  65. C
  66. C        GENERATE TEST MATRIX
  67. C
  68.          CALL STRXX(A,LDA,N,KASE,LUNIT)
  69. C
  70. C        N = 0 SIGNALS NO MORE TEST MATRICES
  71. C
  72. C     ...EXIT
  73.          IF (N .LE. 0) GO TO 360
  74.          ANORM = 0.0E0
  75.          DO 30 J = 1, N
  76.             ANORM = AMAX1(ANORM,SASUM(N,A(1,J),1))
  77.    30    CONTINUE
  78.          WRITE (LUNIT,570) ANORM
  79. C
  80.          IF (N .GT. NPRINT) GO TO 50
  81.             WRITE (LUNIT,370)
  82.             DO 40 I = 1, N
  83.                WRITE (LUNIT,600) (A(I,J), J = 1, N)
  84.    40       CONTINUE
  85.             WRITE (LUNIT,370)
  86.    50    CONTINUE
  87. C
  88. C        GENERATE EXACT SOLUTION
  89. C
  90.          XEXACT(1) = 1.0E0
  91.          IF (N .GE. 2) XEXACT(2) = 0.0E0
  92.          IF (N .LE. 2) GO TO 70
  93.             DO 60 I = 3, N
  94.                XEXACT(I) = -XEXACT(I-2)
  95.    60       CONTINUE
  96.    70    CONTINUE
  97. C
  98. C        GENERATE R.H.S.
  99. C
  100.          DO 90 I = 1, N
  101.             B(I) = 0.0E0
  102.             BT(I) = 0.0E0
  103.             DO 80 J = 1, N
  104.                B(I) = B(I) + A(I,J)*XEXACT(J)
  105.                BT(I) = BT(I) + A(J,I)*XEXACT(J)
  106.    80       CONTINUE
  107.             X(I) = B(I)
  108.             XT(I) = BT(I)
  109.    90    CONTINUE
  110. C
  111. C        UPPER OR LOWER TRIANGULAR
  112. C
  113.          ML = 0
  114.          MU = 0
  115.          DO 120 J = 1, N
  116.             DO 110 I = 1, N
  117.                IF (A(I,J) .EQ. 0.0E0) GO TO 100
  118.                   IF (I .LT. J) MU = MAX0(MU,J-I)
  119.                   IF (I .GT. J) ML = MAX0(ML,I-J)
  120.   100          CONTINUE
  121.   110       CONTINUE
  122.   120    CONTINUE
  123.          WRITE (LUNIT,670) ML,MU
  124.          IF (ML .NE. 0 .AND. MU .NE. 0) GO TO 350
  125.             IF (MU .EQ. 0) JOB = 0
  126.             IF (ML .EQ. 0) JOB = 1
  127.             IF (JOB .EQ. 0) WRITE (LUNIT,710)
  128.             IF (JOB .EQ. 1) WRITE (LUNIT,720)
  129.             STUFF = 4095.0E0
  130.             DO 140 J = 1, N
  131.                DO 130 I = 1, N
  132.                   IF (I .LT. J .AND. JOB .EQ. 0) A(I,J) = STUFF
  133.                   IF (I .GT. J .AND. JOB .EQ. 1) A(I,J) = STUFF
  134.   130          CONTINUE
  135.   140       CONTINUE
  136. C
  137. C           ESTIMATE CONDITION
  138. C
  139.             CALL STRCO(A,LDA,N,RCOND,Z,JOB)
  140. C
  141. C           OUTPUT NULL VECTOR IF N .LE. NPRINT
  142. C
  143.             IF (N .GT. NPRINT) GO TO 160
  144.                WRITE (LUNIT,620)
  145.                DO 150 I = 1, N
  146.                   WRITE (LUNIT,630) Z(I)
  147.   150          CONTINUE
  148.                WRITE (LUNIT,370)
  149.   160       CONTINUE
  150. C
  151. C
  152. C           TEST FOR SINGULARITY
  153. C
  154.             IF (RCOND .GT. 0.0E0) GO TO 170
  155.                WRITE (LUNIT,610) RCOND
  156.                WRITE (LUNIT,400)
  157.                KSING = KSING + 1
  158.             GO TO 340
  159.   170       CONTINUE
  160.                COND = 1.0E0/RCOND
  161.                WRITE (LUNIT,420) COND
  162.                ONEPX = 1.0E0 + RCOND
  163.                IF (ONEPX .EQ. 1.0E0) WRITE (LUNIT,410)
  164. C
  165. C              COMPUTE INVERSE, DETERMINANT AND COND1 = TRUE CONDITION
  166. C
  167.                DO 190 J = 1, N
  168.                   DO 180 I = 1, N
  169.                      AINV(I,J) = A(I,J)
  170.   180             CONTINUE
  171.   190          CONTINUE
  172.                CALL STRDI(AINV,LDA,N,DET,110+JOB,INFO)
  173.                AINORM = 0.0E0
  174.                DO 200 J = 1, N
  175.                   IF (JOB .EQ. 0)
  176.      *               AINORM = AMAX1(AINORM,SASUM(N-J+1,AINV(J,J),1))
  177.                   IF (JOB .EQ. 1)
  178.      *               AINORM = AMAX1(AINORM,SASUM(J,AINV(1,J),1))
  179.   200          CONTINUE
  180.                COND1 = ANORM*AINORM
  181.                WRITE (LUNIT,430) COND1
  182.                WRITE (LUNIT,650) DET(1)
  183.                WRITE (LUNIT,660) DET(2)
  184. C
  185. C              SOLVE  A*X = B  AND  TRANS(A)*XT = BT
  186. C
  187.                CALL STRSL(A,LDA,N,X,JOB,INFO)
  188.                CALL STRSL(A,LDA,N,XT,JOB+10,INFO)
  189. C
  190.                IF (N .GT. NPRINT) GO TO 230
  191.                   WRITE (LUNIT,440)
  192.                   DO 210 I = 1, N
  193.                      WRITE (LUNIT,640) X(I)
  194.   210             CONTINUE
  195.                   WRITE (LUNIT,450)
  196.                   DO 220 I = 1, N
  197.                      WRITE (LUNIT,640) XT(I)
  198.   220             CONTINUE
  199.                   WRITE (LUNIT,370)
  200.   230          CONTINUE
  201. C
  202. C              RESTORE ZEROS IN OTHER TRIANGLE
  203. C
  204.                DO 260 J = 1, N
  205.                   DO 250 I = 1, N
  206.                      IF (A(I,J) .NE. STUFF) GO TO 240
  207.                         A(I,J) = 0.0E0
  208.                         AINV(I,J) = 0.0E0
  209.   240                CONTINUE
  210.   250             CONTINUE
  211.   260          CONTINUE
  212. C
  213. C              COMPUTE ERRORS AND RESIDUALS
  214. C                 E  =  X - XEXACT
  215. C                 ET =  XT - XEXACT
  216. C                 R  =  B - A*X
  217. C                 RT =  BT - A*XT
  218. C                 AI = A*INV(A) - I
  219. C
  220.                XNORM = SASUM(N,X,1)
  221.                XTNORM = SASUM(N,XT,1)
  222.                ENORM = 0.0E0
  223.                ETNORM = 0.0E0
  224.                FNORM = 0.0E0
  225.                DO 270 J = 1, N
  226.                   ENORM = ENORM + ABS(X(J)-XEXACT(J))
  227.                   ETNORM = ETNORM + ABS(XT(J)-XEXACT(J))
  228.                   T = -X(J)
  229.                   CALL SAXPY(N,T,A(1,J),1,B,1)
  230.                   BT(J) = BT(J) - SDOT(N,A(1,J),1,XT,1)
  231.   270          CONTINUE
  232.                RNORM = SASUM(N,B,1)
  233.                RTNORM = SASUM(N,BT,1)
  234. C
  235. C
  236. C              A*INV(A) - I
  237. C
  238.                AINORM = 0.0E0
  239.                DO 300 J = 1, N
  240.                   DO 280 I = 1, N
  241.                      B(I) = 0.0E0
  242.   280             CONTINUE
  243.                   DO 290 K = 1, N
  244.                      T = AINV(K,J)
  245.                      CALL SAXPY(N,T,A(1,K),1,B,1)
  246.   290             CONTINUE
  247.                   B(J) = B(J) - 1.0E0
  248.                   AINORM = AMAX1(AINORM,SASUM(N,B,1))
  249.   300          CONTINUE
  250. C
  251.                WRITE (LUNIT,460) ENORM,ETNORM
  252.                WRITE (LUNIT,470) RNORM,RTNORM
  253.                WRITE (LUNIT,580) AINORM
  254. C
  255. C              COMPUTE TEST RATIOS
  256. C
  257.                Q(1) = COND/COND1
  258.                Q(2) = COND1/COND
  259.                Q(3) = ENORM/(EPS*COND*XNORM)
  260.                Q(4) = ETNORM/(EPS*COND*XTNORM)
  261.                Q(5) = RNORM/(EPS*ANORM*XNORM)
  262.                Q(6) = RTNORM/(EPS*ANORM*XTNORM)
  263.                Q(7) = AINORM/(EPS*COND)
  264.                WRITE (LUNIT,370)
  265.                WRITE (LUNIT,480)
  266.                WRITE (LUNIT,370)
  267.                WRITE (LUNIT,540)
  268.                WRITE (LUNIT,550)
  269.                WRITE (LUNIT,560)
  270.                WRITE (LUNIT,370)
  271.                WRITE (LUNIT,590) (Q(I), I = 1, 7)
  272.                WRITE (LUNIT,370)
  273. C
  274. C              LOOK FOR SUSPICIOUS RATIOS
  275. C
  276.                QS(1) = 1.0E0 + 4.0E0*EPS
  277.                QS(2) = 10.0E0
  278.                EN = FLOAT(N)
  279.                IF (N .EQ. 1) EN = 2.0E0
  280.                DO 310 I = 3, 7
  281.                   QS(I) = EN
  282.   310          CONTINUE
  283.                KOUNT = 0
  284.                DO 330 I = 1, 7
  285.                   IQ(I) = 0
  286.                   IF (Q(I) .LE. QS(I)) GO TO 320
  287.                      IQ(I) = 1
  288.                      KSUSP(I) = KSUSP(I) + 1
  289.                      KOUNT = KOUNT + 1
  290.   320             CONTINUE
  291.   330          CONTINUE
  292.                IF (KOUNT .EQ. 0) WRITE (LUNIT,690)
  293.                IF (KOUNT .NE. 0) WRITE (LUNIT,700) (IQ(I), I = 1, 7)
  294.                WRITE (LUNIT,370)
  295.   340       CONTINUE
  296.   350    CONTINUE
  297. C
  298.          WRITE (LUNIT,490)
  299.          KASE = KASE + 1
  300.       GO TO 20
  301.   360 CONTINUE
  302. C
  303. C     FINISH MAIN LOOP
  304. C
  305. C     SUMMARY
  306. C
  307.       WRITE (LUNIT,500)
  308.       KASE = KASE - 1
  309.       WRITE (LUNIT,510) KASE
  310.       WRITE (LUNIT,520) KSING
  311.       WRITE (LUNIT,530) KSUSP
  312.       WRITE (LUNIT,680)
  313.       RETURN
  314. C
  315. C     MOST FORMATS, ALSO SOME IN STRXX
  316. C
  317.   370 FORMAT (1H )
  318.   380 FORMAT (22H1LINPACK TESTER, STR**)
  319.   390 FORMAT ( / 18H MACHINE EPSILON =, 1PE13.5)
  320.   400 FORMAT ( / 19H EXACT SINGULARITY. /)
  321.   410 FORMAT ( / 16H MAYBE SINGULAR. /)
  322.   420 FORMAT (14H COND        =, 1PE13.5)
  323.   430 FORMAT (14H ACTUAL COND =, 1PE13.5)
  324.   440 FORMAT ( / 4H X =)
  325.   450 FORMAT ( / 5H XT =)
  326.   460 FORMAT (14H ERROR NORMS =, 1P2E13.5)
  327.   470 FORMAT (14H RESID NORMS =, 1P2E13.5)
  328.   480 FORMAT (26H TEST RATIOS.. E = MACHEPS)
  329.   490 FORMAT ( / 14H ************* /)
  330.   500 FORMAT (8H1SUMMARY)
  331.   510 FORMAT (18H NUMBER OF TESTS =, I4)
  332.   520 FORMAT (30H NUMBER OF SINGULAR MATRICES =, I4)
  333.   530 FORMAT (30H NUMBER OF SUSPICIOUS RATIOS =, 7I4)
  334.   540 FORMAT (30H     COND     ACTUAL    ERROR ,
  335.      *        40H   ERROR-T    RESID    RESID-T   A*AI-I )
  336.   550 FORMAT (7(10H   -------))
  337.   560 FORMAT (30H    ACTUAL     COND   E*COND*X,
  338.      *        40H  E*COND*X    E*A*X     E*A*X    E*COND )
  339.   570 FORMAT (14H NORM(A)     =, 1PE13.5)
  340.   580 FORMAT (14H NORM(A*AI-I)=, 1PE13.5)
  341.   590 FORMAT (7F10.4)
  342.   600 FORMAT (1H , 6G11.4)
  343.   610 FORMAT (14H 1/COND      =, 1PE13.5)
  344.   620 FORMAT ( / 7H NULL =)
  345.   630 FORMAT (2G14.6)
  346.   640 FORMAT (2G14.6)
  347.   650 FORMAT (14H DET FRACT   =, 2F9.5)
  348.   660 FORMAT (14H DET EXPON   =, 2F9.0)
  349.   670 FORMAT (5H ML =, I2, 6H  MU =, I2)
  350.   680 FORMAT ( / 12H END OF TEST)
  351.   690 FORMAT (21H NO SUSPICIOUS RATIOS)
  352.   700 FORMAT (I8, 5I10 / 7X, 28H1 INDICATES SUSPICIOUS RATIO)
  353.   710 FORMAT (26H LOWER TRIANGULAR, JOB = 0)
  354.   720 FORMAT (26H UPPER TRIANGULAR, JOB = 1)
  355.   730 FORMAT (29H THIS VERSION DATED 08/14/78.)
  356.       END
  357.       SUBROUTINE STRXX(A,LDA,N,KASE,LUNIT)
  358.       INTEGER LDA,N,KASE,LUNIT
  359.       REAL A(LDA,1)
  360. C
  361. C     GENERATES REAL TRIANGULAR TEST MATRICES
  362. C
  363. C     EXTERNAL SMACH
  364. C     FORTRAN FLOAT
  365.       REAL T1,T2
  366.       REAL SMACH,HUGE,TINY
  367.       INTEGER I,J
  368. C
  369.       GO TO (10,10,10,50,50,70,70,70,110,150,200,240,280,320,360,410,
  370.      *       460), KASE
  371. C
  372. C     KASE 1, 2 AND 3
  373. C
  374.    10 CONTINUE
  375.          N = 3*KASE
  376.          WRITE (LUNIT,20) KASE,N
  377.    20    FORMAT (5H KASE, I3, 3X, 16HHILBERT-HALF     / 4H N =, I4)
  378.          DO 40 J = 1, N
  379.             DO 30 I = 1, N
  380.                A(I,J) = 0.0E0
  381.                IF (I .GE. J - 3 .AND. I .LE. J)
  382.      *            A(I,J) = 1.0E0/FLOAT(I+J-1)
  383.    30       CONTINUE
  384.    40    CONTINUE
  385.       GO TO 470
  386. C
  387. C     KASE 4 AND 5
  388. C
  389.    50 CONTINUE
  390.          N = 1
  391.          WRITE (LUNIT,60) KASE,N
  392.    60    FORMAT (5H KASE, I3, 3X, 16HMONOELEMENTAL    / 4H N =, I4)
  393.          IF (KASE .EQ. 4) A(1,1) = 3.0E0
  394.          IF (KASE .EQ. 5) A(1,1) = 0.0E0
  395.       GO TO 470
  396. C
  397. C     KASE 6, 7 AND 8
  398. C
  399.    70 CONTINUE
  400.          N = 15
  401.          WRITE (LUNIT,80) KASE,N
  402.    80    FORMAT (5H KASE, I3, 3X, 16HBIDIAGONAL       / 4H N =, I4)
  403.          T1 = 0.0E0
  404.          T2 = 0.0E0
  405.          IF (KASE .EQ. 7) T1 = 100.0E0
  406.          IF (KASE .EQ. 8) T2 = 100.0E0
  407.          DO 100 I = 1, N
  408.             DO 90 J = 1, N
  409.                A(I,J) = 0.0E0
  410.                IF (I .EQ. J) A(I,I) = 4.0E0
  411.                IF (I .EQ. J - 1) A(I,J) = T1
  412.                IF (I .EQ. J + 1) A(I,J) = T2
  413.    90       CONTINUE
  414.   100    CONTINUE
  415.       GO TO 470
  416. C
  417. C     KASE 9
  418. C
  419.   110 CONTINUE
  420.          N = 5
  421.          WRITE (LUNIT,120) KASE,N
  422.   120    FORMAT (5H KASE, I3, 3X, 16HHALF OF RANK ONE / 4H N =, I4)
  423.          DO 140 I = 1, N
  424.             DO 130 J = 1, N
  425.                A(I,J) = 0.0E0
  426.                IF (I .GE. J) A(I,J) = 10.0E0**(I - J)
  427.   130       CONTINUE
  428.   140    CONTINUE
  429.       GO TO 470
  430. C
  431. C     KASE 10
  432. C
  433.   150 CONTINUE
  434.          N = 4
  435.          WRITE (LUNIT,160) KASE,N
  436.   160    FORMAT (5H KASE, I3, 3X, 16HZERO COLUMN      / 4H N =, I4)
  437.          DO 190 I = 1, N
  438.             DO 180 J = 1, N
  439.                A(I,J) = 0.0E0
  440.                IF (I .LT. J) GO TO 170
  441.                   T1 = FLOAT(J-3)
  442.                   T2 = FLOAT(I)
  443.                   A(I,J) = T1/T2
  444.   170          CONTINUE
  445.   180       CONTINUE
  446.   190    CONTINUE
  447.       GO TO 470
  448. C
  449. C     KASE 11
  450. C
  451.   200 CONTINUE
  452.          N = 5
  453.          WRITE (LUNIT,210) KASE,N
  454.   210    FORMAT (5H KASE, I3, 3X, 16HTEST COND        / 4H N =, I4)
  455.          DO 230 I = 1, N
  456.             DO 220 J = 1, N
  457.                IF (I .EQ. J) A(I,J) = 1.0E0
  458.                IF (I .GT. J) A(I,J) = 0.0E0
  459.                IF (I .LT. J) A(I,J) = -1.0E0
  460.   220       CONTINUE
  461.   230    CONTINUE
  462.       GO TO 470
  463. C
  464. C     KASE 12
  465. C
  466.   240 CONTINUE
  467.          N = 3
  468.          WRITE (LUNIT,250) KASE,N
  469.   250    FORMAT (5H KASE, I3, 3X, 16HIDENTITY         / 4H N =, I4)
  470.          DO 270 I = 1, N
  471.             DO 260 J = 1, N
  472.                IF (I .EQ. J) A(I,I) = 1.0E0
  473.                IF (I .NE. J) A(I,J) = 0.0E0
  474.   260       CONTINUE
  475.   270    CONTINUE
  476.       GO TO 470
  477. C
  478. C     KASE 13
  479. C
  480.   280 CONTINUE
  481.          N = 6
  482.          WRITE (LUNIT,290) KASE,N
  483.   290    FORMAT (5H KASE, I3, 3X, 16HUPPER TRIANGULAR / 4H N =, I4)
  484.          DO 310 I = 1, N
  485.             DO 300 J = 1, N
  486.                IF (I .GT. J) A(I,J) = 0.0E0
  487.                IF (I .LE. J) A(I,J) = FLOAT(I-J-1)*1.0E0
  488.   300       CONTINUE
  489.   310    CONTINUE
  490.       GO TO 470
  491. C
  492. C     KASE 14
  493. C
  494.   320 CONTINUE
  495.          N = 6
  496.          WRITE (LUNIT,330) KASE,N
  497.   330    FORMAT (5H KASE, I3, 3X, 16HLOWER TRIANGULAR / 4H N =, I4)
  498.          DO 350 I = 1, N
  499.             DO 340 J = 1, N
  500.                IF (I .LT. J) A(I,J) = 0.0E0
  501.                IF (I .GE. J) A(I,J) = FLOAT(I-J+1)*1.0E0
  502.   340       CONTINUE
  503.   350    CONTINUE
  504.       GO TO 470
  505. C
  506. C     KASE 15
  507. C
  508.   360 CONTINUE
  509.          N = 5
  510.          WRITE (LUNIT,370) KASE,N
  511.   370    FORMAT (5H KASE, I3, 3X, 16HNEAR UNDERFLOW   / 4H N =, I4)
  512.          TINY = SMACH(2)
  513.          WRITE (LUNIT,380) TINY
  514.   380    FORMAT (14H TINY        =, 1PE13.5)
  515.          DO 400 I = 1, N
  516.             DO 390 J = 1, N
  517.                A(I,J) = 0.0E0
  518.                IF (I .LE. J) A(I,J) = TINY*(FLOAT(J)/FLOAT(I))
  519.   390       CONTINUE
  520.   400    CONTINUE
  521.       GO TO 470
  522. C
  523. C     KASE 16
  524. C
  525.   410 CONTINUE
  526.          N = 5
  527.          WRITE (LUNIT,420) KASE,N
  528.   420    FORMAT (5H KASE, I3, 3X, 16HNEAR OVERFLOW    / 4H N =, I4)
  529.          HUGE = SMACH(3)
  530.          WRITE (LUNIT,430) HUGE
  531.   430    FORMAT (14H HUGE        =, 1PE16.5)
  532.          DO 450 I = 1, N
  533.             DO 440 J = 1, N
  534.                A(I,J) = 0.0E0
  535.                IF (I .GE. J) A(I,J) = HUGE*(FLOAT(J)/FLOAT(I))
  536.   440       CONTINUE
  537.   450    CONTINUE
  538.       GO TO 470
  539. C
  540.   460 CONTINUE
  541.          N = 0
  542.   470 CONTINUE
  543.       RETURN
  544. C
  545.       END
  546.