home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / tests / rsttest.f < prev    next >
Text File  |  1996-09-28  |  15KB  |  444 lines
  1. C
  2. C     THIS DRIVER TESTS  EISPACK  FOR THE CLASS OF REAL SYMMETRIC TRI-
  3. C     DIAGONAL MATRICES SUMMARIZING THE FIGURES OF MERIT FOR ALL PATHS.
  4. C
  5. C     THIS DRIVER IS CATALOGUED AS  EISPDRV4(RSTSUMAR).
  6. C
  7. C     THE DIMENSION OF  ST  SHOULD BE  NM  BY  2.
  8. C     THE DIMENSION OF  Z  SHOULD BE  NM  BY  NM.
  9. C     THE DIMENSION OF  W,E,E2,IND,RV1,RV2,RV3,RV4,RV5,RV6,
  10. C     W1,  AND  W2  SHOULD BE  NM.
  11. C     THE DIMENSION OF  STHOLD  SHOULD BE  NM  BY  2.
  12. C     HERE NM = 20.
  13. C
  14.       DOUBLE PRECISION Z(20,20),ST(20,2),STHOLD(20,2),W(20),E(20),
  15.      X        E2(20),RV1(20),RV2(20),RV3(20),RV4(20),RV5(20),RV6(20),
  16.      X        W1( 20),W2( 20),TCRIT( 8),EPSLON,RESDUL,MAXEIG,MAXDIF,
  17.      X        U,LB,UB,EPS1,DFL
  18.       REAL  XUB,XLB
  19.       INTEGER IND( 20),IERR( 6),ERROR
  20.       DATA IREAD1/1/,IREADC/5/,IWRITE/6/
  21. C
  22.       OPEN(UNIT=IREAD1,FILE='FILE51')
  23.       OPEN(UNIT=IREADC,FILE='FILE52')
  24.       REWIND IREAD1
  25.       REWIND IREADC
  26. C
  27.       NM = 20
  28.       MMB = 2
  29.       LCOUNT = 0
  30.       WRITE(IWRITE,1)
  31.     1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S
  32.      XTATISTICS//1H ,95(1H-)/96H ORDER TQL2   TQL1   IMTQL2 IMTQL1    LB
  33.      X      UB    M  IMTQLV   TSTURM   BISECT  M1 NO   RATQR  /1H ,
  34.      X95(1H-)//48H UNDER 'ORDER' IS THE ORDER OF EACH TEST MATRIX. //
  35.      X95H UNDER 'TQL2   TQL1' ARE THREE NUMBERS.  THE FIRST NUMBER, AN I
  36.      XNTEGER, IS THE ABSOLUTE SUM OF  /
  37.      X95H THE ERROR FLAGS RETURNED SEPARATELY FROM  TQL2  AND  TQL1.  TH
  38.      XE SECOND NUMBER IS THE MEASURE  /
  39.      X62H OF PERFORMANCE BASED UPON THE RESIDUAL COMPUTED FOR THE  TQL2,
  40.      X25H  PATH.  THE THIRD NUMBER        /
  41.      X62H MEASURES THE AGREEMENT OF THE EIGENVALUES FROM THE  TQL2  AND,
  42.      X14H  TQL1  PATHS.    //
  43.      X95H UNDER 'IMTQL2 IMTQL1' ARE THREE NUMBERS WITH MEANING LIKE THOS
  44.      XE UNDER  'TQL2   TQL1'.         //
  45.      X95H UNDER 'LB' AND 'UB' ARE THE INPUT VARIABLES SPECIFYING THE INT
  46.      XERVAL TO  BISECT  AND  TSTURM.  //
  47.      X61H UNDER 'M' IS THE NUMBER OF EIGENVALUES DETERMINED BY  BISECT,
  48.      X30H  AND  TSTURM  THAT LIE IN THE    /18H INTERVAL (LB,UB).//
  49.      X95H UNDER EACH OF 'IMTQLV', 'TSTURM', 'BISECT', AND 'RATQR' ARE TW
  50.      XO NUMBERS.  THE FIRST NUMBER,       )
  51.       WRITE(IWRITE,2)
  52.     2 FORMAT(
  53.      X95H AN INTEGER, IS THE ABSOLUTE SUM OF THE ERROR FLAGS RETURNED FR
  54.      XOM THE RESPECTIVE PATH.         /
  55.      X95H THE SECOND NUMBER IS THE MEASURE OF PERFORMANCE BASED UPON THE
  56.      X RESIDUAL COMPUTED FOR THE PATH.//
  57.      X95H UNDER 'M1' AND 'NO' ARE THE VARIABLES SPECIFYING THE LOWER BOU
  58.      XNDARY INDEX AND THE NUMBER       /
  59.      X28H OF EIGENVALUES TO  RATQR.  //
  60.      X62H -1.0  AS THE MEASURE OF PERFORMANCE IS PRINTED IF AN ERROR IN,
  61.      X27H THE CORRESPONDING PATH HAS        /
  62.      X47H PREVENTED THE COMPUTATION OF THE EIGENVECTORS. //
  63.      X62H THE  TQL2    PATH USES THE EISPACK CODES  TQL2  .            /
  64.      X62H THE  TQL1    PATH USES THE EISPACK CODES  TQL1  .            /
  65.      X62H THE  IMTQL2  PATH USES THE EISPACK CODES  IMTQL2,            /
  66.      X39H AS CALLED FROM DRIVER SUBROUTINE  RST. /
  67.      X62H THE  IMTQL1  PATH USES THE EISPACK CODES  IMTQL1,            /
  68.      X39H AS CALLED FROM DRIVER SUBROUTINE  RST. /
  69.      X57H THE  IMTQLV  PATH USES THE EISPACK CODES  IMTQLV-TINVIT.
  70.      X            )
  71.       WRITE(IWRITE,3)
  72.     3 FORMAT(
  73.      X50H THE  TSTURM  PATH USES THE EISPACK CODES  TSTURM.
  74.      X   /
  75.      X57H THE  BISECT  PATH USES THE EISPACK CODES  BISECT-TINVIT.
  76.      X             /
  77.      X57H THE  RATQR   PATH USES THE EISPACK CODES  RATQR -TINVIT.
  78.      X             /)
  79.       WRITE(IWRITE,15)
  80.    15 FORMAT(1X,21HD.P. VERSION 04/15/83 )
  81.     5 FORMAT( 53H1       TABULATION OF THE ERROR FLAG  ERROR  AND THE ,
  82.      X    31HMEASURE OF PERFORMANCE  Y  FOR /5X,
  83.      X    56HTHE  EISPACK  CODES.  THIS RUN DISPLAYS THESE STATISTICS ,
  84.      X    41H FOR REAL SYMMETRIC TRIDIAGONAL MATRICES.  /
  85.      X    55H0ORDER TQL2   TQL1   IMTQL2 IMTQL1    LB      UB    M   ,
  86.      X    40HIMTQLV   TSTURM   BISECT  M1 NO   RATQR )
  87.    10 CALL RMATIN(NM,MMB,N,MB,ST,STHOLD,0)
  88.       READ(IREADC,50) MM,LB,UB,M11,NO
  89.    50 FORMAT(I4,2D24.16,2(4X,I4))
  90. C
  91. C     MM,LB,UB,M11,  AND  NO  ARE READ FROM SYSIN AFTER THE MATRIX IS
  92. C     GENERATED.  MM,LB,  AND  UB  SPECIFY TO  BISECT  THE MAXIMUM
  93. C     NUMBER OF EIGENVALUES AND BOUNDS FOR THE INTERVAL WHICH IS TO
  94. C     BE SEARCHED.  M11  AND  NO  SPECIFY TO  RATQR  THE LOWER BOUNDARY
  95. C     INDEX AND THE NUMBER OF DESIRED EIGENVALUES.
  96. C
  97.       DO  230  ICALL = 1,10
  98.          IF( ICALL .NE. 1 )  CALL  RMATIN(NM,MMB,N,MB,ST,STHOLD,1)
  99. C
  100. C
  101. C     IF  TQL1  PATH (LABEL 80) IS TAKEN THEN  TQL2  PATH (LABEL 70)
  102. C     MUST ALSO BE TAKEN IN ORDER THAT THE MEASURE OF PERFORMANCE BE
  103. C     MEANINGFUL.
  104. C     IF  IMTQL1  PATH (LABEL 85) IS TAKEN THEN  IMTQL2  PATH (LABEL 75)
  105. C     MUST ALSO BE TAKEN IN ORDER THAT THE MEASURE OF PERFORMANCE BE
  106. C     MEANINGFUL.
  107. C     IF  TQL2  (IMTQL2)  PATH FAILS, THEN  TQL1  (IMTQL1)  PATH IS
  108. C     OMITTED AND PRINTOUT FLAGGED WITH  -1.0.
  109. C
  110.          GO TO  (70,75,80,85,89,90,95,230,110,230),  ICALL
  111. C
  112. C     RSTWZ  USING  TQL2
  113. C
  114.    70    ICT = 1
  115.          DO  72  I = 1,N
  116.             DO  71  J = 1,N
  117.    71         Z(I,J) = 0.0D0
  118.             Z(I,I) = 1.0D0
  119.             E(I) = ST(I,1)
  120.    72       W(I) = ST(I,2)
  121.          CALL  TQL2(NM,N,W,E,Z,ERROR)
  122.          IERR(ICT) = ERROR
  123.          M = ERROR - 1
  124.          IF( ERROR .NE. 0 ) GO TO 190
  125.          DO 73 I = 1,N
  126.    73       W1(I) = W(I)
  127.          M = N
  128.          GO TO  190
  129. C
  130. C     RSTWZ  USING  IMTQL2
  131. C     INVOKED FROM DRIVER SUBROUTINE  RST.
  132. C
  133.    75    ICT = 2
  134.          DO  77  I = 1,N
  135.             E(I) = ST(I,1)
  136.    77       W(I) = ST(I,2)
  137.          CALL  RST(NM,N,W,E,1,Z,ERROR)
  138.          IERR(ICT) = ERROR
  139.          M = ERROR - 1
  140.          IF( ERROR .NE. 0 ) GO TO 190
  141.          DO 78 I = 1,N
  142.             W2(I) = W(I)
  143.    78    CONTINUE
  144.          M = N
  145.          GO TO  190
  146. C
  147. C     RSTW  USING  TQL1
  148. C
  149.    80    ICT = 7
  150.          IF( IERR(1) .NE. 0 ) GO TO 200
  151.          DO  81  I = 1,N
  152.    81       W(I) = ST(I,2)
  153.          CALL  TQL1(N,W,ST(1,1),ERROR)
  154.          IERR(1) = ERROR
  155.          IF( ERROR .NE. 0 ) GO TO 200
  156.          MAXEIG = 0.0D0
  157.          MAXDIF = 0.0D0
  158.          DO 82 I = 1,N
  159.             IF( DABS(W(I)) .GT. MAXEIG ) MAXEIG = DABS(W(I))
  160.             U = DABS(W1(I) - W(I))
  161.             IF( U .GT. MAXDIF ) MAXDIF = U
  162.    82    CONTINUE
  163.          IF( MAXEIG .EQ. 0.0D0 ) MAXEIG = 1.0D0
  164.          DFL = N * 10
  165.          TCRIT(7) = MAXDIF/EPSLON(MAXEIG*DFL)
  166.          GO TO  230
  167. C
  168. C     RSTW  USING  IMTQL1
  169. C     INVOKED FROM DRIVER SUBROUTINE  RST.
  170. C
  171.    85    ICT = 8
  172.          IF( IERR(2) .NE. 0 ) GO TO 200
  173.          DO  86  I = 1,N
  174.          E(I) = ST(I,1)
  175.    86       W(I) = ST(I,2)
  176.          CALL  RST(NM,N,W,E,0,Z,ERROR)
  177.          IERR(2) = ERROR
  178.          IF( ERROR .NE. 0 ) GO TO 200
  179.          MAXEIG = 0.0D0
  180.          MAXDIF = 0.0D0
  181.          DO 87 I = 1,N
  182.             IF( DABS(W(I)) .GT. MAXEIG ) MAXEIG = DABS(W(I))
  183.             U = DABS(W2(I) - W(I))
  184.             IF( U .GT. MAXDIF ) MAXDIF = U
  185.    87    CONTINUE
  186.          IF( MAXEIG .EQ. 0.0D0 ) MAXEIG = 1.0D0
  187.          DFL = N * 10
  188.          TCRIT(8) = MAXDIF/EPSLON(MAXEIG*DFL)
  189.          GO TO  230
  190. C
  191. C     RSTW1Z  ( USAGE HERE COMPUTES ALL EIGENVECTORS )
  192. C
  193.    89    ICT = 3
  194.          DO 891 I = 2,N
  195.              E2(I) = ST(I,1) ** 2
  196.   891    CONTINUE
  197.          CALL IMTQLV(N,ST(1,2),ST(1,1),E2,W,IND,ERROR,RV1)
  198.          IERR(ICT) = ERROR
  199.          M = N
  200.          IF( ERROR .NE. 0 ) M = ERROR - 1
  201.   893    CALL TINVIT(NM,N,ST(1,2),ST(1,1),E2,M,W,IND,Z,ERROR,RV1,RV2,
  202.      X               RV3,RV4,RV6)
  203.          IERR(ICT) = IERR(ICT) + IABS(ERROR)
  204.          GO TO 190
  205. C
  206. C     RST1W1Z  USING  TSTURM
  207. C
  208.    90    ICT = 4
  209.          EPS1 = 0.0D0
  210.          DO  92  I = 2,N
  211.    92       E2(I) = ST(I,1)**2
  212.          CALL  TSTURM(NM,N,EPS1,ST(1,2),ST(1,1),E2,LB,UB,MM,M,W,Z,
  213.      X                ERROR,RV1,RV2,RV3,RV4,RV5,RV6)
  214.          IERR(ICT) = ERROR
  215.          XLB = LB
  216.          XUB = UB
  217.          IF( ERROR .EQ. 3*N + 1) GO TO 200
  218.          IF( ERROR .GT. 4*N ) M = ERROR - 4*N - 1
  219.          GO TO  190
  220. C
  221. C     RST1W1Z  USING  BISECT  AND  TINVIT
  222. C
  223.    95    ICT = 5
  224.          EPS1 = 0.0D0
  225.          DO  97  I = 2,N
  226.    97       E2(I) = ST(I,1)**2
  227.          CALL  BISECT(N,EPS1,ST(1,2),ST(1,1),E2,LB,UB,MM,M,W,IND,ERROR,
  228.      X                RV4,RV5)
  229.          IERR(ICT) = ERROR
  230.          MBISCT = M
  231.          XLB = LB
  232.          XUB = UB
  233.          IF( ERROR .NE. 0 ) GO TO 200
  234.          CALL  TINVIT(NM,N,ST(1,2),ST(1,1),E2,M,W,IND,Z,
  235.      X                ERROR,RV1,RV2,RV3,RV4,RV6)
  236.          IERR(ICT) = IABS(ERROR)
  237.          GO TO  190
  238. C
  239. C     RST1W1Z  USING  RATQR  AND  TINVIT
  240. C
  241.   110    ICT = 6
  242.          EPS1 = 0.0D0
  243.          DO 112  I = 2,N
  244.   112       E2(I) = ST(I,1)**2
  245.          CALL  RATQR(N,EPS1,ST(1,2),ST(1,1),E2,NO+M11-1,W,IND,
  246.      X               RV1,.TRUE.,0,ERROR)
  247.          DO 115 I = 1,NO
  248.             M = I + M11 - 1
  249.             W(I) = W(M)
  250.             IND(I) = IND(M)
  251.   115    CONTINUE
  252.          IERR(ICT) = ERROR
  253.          IF( ERROR .NE. 0 )  GO TO  200
  254.          M = NO
  255.          CALL  TINVIT(NM,N,ST(1,2),ST(1,1),E2,M,W,IND,Z,ERROR,
  256.      X                   RV1,RV2,RV3,RV4,RV6)
  257.          IERR(ICT) = IABS(ERROR)
  258. C
  259.   190    IF( M .EQ. 0 .AND. ERROR .NE. 0 ) GO TO 200
  260.          CALL  RSBWZR(NM,N,M,2,ST,W,Z,RV1,RESDUL)
  261.          DFL = N * 10
  262.          TCRIT(ICT) = RESDUL/EPSLON(DFL)
  263.          GO TO 230
  264.   200    TCRIT(ICT) = -1.0D0
  265.   230 CONTINUE
  266. C
  267.       IF( MOD(LCOUNT,35) .EQ. 0 ) WRITE(IWRITE,5)
  268.       LCOUNT = LCOUNT + 1
  269.       WRITE(IWRITE,240) N,IERR(1),TCRIT(1),TCRIT(7),IERR(2),TCRIT(2),
  270.      X             TCRIT(8),XLB,XUB,MBISCT,(IERR(I),TCRIT(I),I=3,5),
  271.      X             M11,NO,IERR(6),TCRIT(6)
  272.   240 FORMAT(I4,2(I3,2F6.3),2(1PE8.0),I3,3(I3,0PF6.3),3I3,F6.3)
  273.       GO TO  10
  274.       END
  275.       SUBROUTINE RSBWZR(NM,N,M,MB,A,W,Z,NORM,RESDUL)
  276. C
  277.       INTEGER  MB,MB1
  278.       DOUBLE PRECISION NORM(M), W(M), A(NM,MB), Z(NM,M), NORMA, TNORM,
  279.      X       S, SUM, SUMA, SUMZ, RESDUL
  280. C
  281. C     THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
  282. C     A*Z-Z*DIAG(W)  WHERE  A  IS A REAL SYMMETRIC BAND MATRIX,
  283. C     W  IS A VECTOR WHICH CONTAINS  M  EIGENVALUES OF  A, AND  Z
  284. C     IS AN ARRAY WHICH CONTAINS THE  M  CORRESPONDING EIGENVECTORS OF
  285. C     A.  ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS.
  286. C
  287. C     THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RSBWZR).
  288. C
  289. C     INPUT.
  290. C
  291. C        NM IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
  292. C           AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT;
  293. C
  294. C        N IS THE ORDER OF THE MATRIX  A;
  295. C
  296. C        M IS THE NUMBERS OF EIGENVECTORS WHOSE RESIDUALS ARE DESIRED;
  297. C
  298. C        MB  IS THE BAND WIDTH OF THE INPUT MATRIX  A .  BAND WIDTH IS
  299. C            DEFINED AS THE NUMBER OF ADJACENT DIAGONALS, INCLUDING THE
  300. C            PRINCIPAL DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO
  301. C            PORTION OF THE LOWER TRIANGLE OF THE MATRIX;
  302. C
  303. C        A(N,MB) IS AN ARRAY WHICH CONTAINS IN ITS COLUMNS THE
  304. C           SUBDIAGONAL AND DIAGONAL OF THE SYMMETRIC BAND
  305. C           MATRIX;
  306. C
  307. C        W(M) IS A VECTOR WHOSE FIRST  M  COMPONENTS CONTAIN EIGENVALUES
  308. C           OF  A;
  309. C
  310. C        Z(NM,M) IS AN ARRAY WHOSE FIRST  M  COLUMNS CONTAIN THE
  311. C           EIGENVECTORS OF  A  CORRESPONDING TO THE EIGENVALUES IN  W.
  312. C
  313. C     OUTPUT.
  314. C
  315. C        Z(NM,M) IS AN ARRAY WHICH CONTAINS THE NORMALIZED
  316. C           APPROXIMATE EIGENVECTORS OF  A.  THE EIGENVECTORS
  317. C           ARE NORMALIZED USING THE 1-NORM IN SUCH A WAY
  318. C           THAT THE FIRST ELEMENT WHOSE MAGNITUDE IS LARGER
  319. C           THAN THE NORM OF THE EIGENVECTOR DIVIDED BY  N  IS
  320. C           POSITIVE;
  321. C
  322. C        NORM(M) IS AN ARRAY SUCH THAT FOR EACH  K
  323. C           NORM(K) = !!A*Z(K)-Z(K)*W(K)!!/(!!A!!*!!Z(K)!!)
  324. C           WHERE  Z(K)  IS THE K-TH EIGENVECTOR;
  325. C
  326. C        RESDUL IS THE REAL NUMBER
  327. C           !!A*Z-Z*DIAG(W)!!/(!!A!!*!!Z!!).
  328. C
  329. C     ----------------------------------------------------------------
  330. C
  331.       RESDUL = 0.0D0
  332.       IF( M .EQ. 0 ) RETURN
  333.       NORMA = 0.0D0
  334. C
  335.       DO 40 I=1,N
  336.          J=I
  337.          SUMA = 0.0D0
  338.          IF(I .EQ. 1) GO TO 20
  339.          MB1 = MB-1
  340.          LSTART = MAX0(1,MB+1-I)
  341. C
  342.          DO 10 L=LSTART,MB1
  343.    10       SUMA = SUMA + DABS(A(I,L))
  344. C
  345.    20    LSTOP = MIN0(MB,N+1-J)
  346. C
  347.          DO 30 L=1,LSTOP
  348.             L1 = MB + 1 - L
  349.             SUMA = SUMA + DABS(A(J,L1))
  350.    30       J = J+1
  351. C
  352.    40    NORMA = DMAX1(SUMA,NORMA)
  353. C
  354.       IF(NORMA .EQ. 0.0D0) NORMA = 1.0D0
  355. C
  356.       DO 120 I=1,M
  357.          S = 0.0D0
  358.          SUMZ = 0.0D0
  359. C
  360.          DO 80 L=1,N
  361.             SUM = -W(I)*Z(L,I)
  362.             SUMZ = SUMZ + DABS(Z(L,I))
  363.             J = MAX0(0,L-MB)
  364.             IF(L .EQ. 1) GO TO 60
  365.             MB1 = MB-1
  366.             KSTART = MAX0(1,MB+1-L)
  367. C
  368.             DO 50 K=KSTART,MB1
  369.                J = J+1
  370.    50          SUM = SUM + A(L,K)*Z(J,I)
  371. C
  372.    60       KSTOP = MIN0(MB,N+1-L)
  373. C
  374.             DO 70 K=1,KSTOP
  375.                J = J+1
  376.                K1 = MB + 1 - K
  377.    70          SUM = SUM + A(J,K1)*Z(J,I)
  378. C
  379.    80       S = S + DABS(SUM)
  380. C
  381.          NORM(I) = SUMZ
  382.          IF (SUMZ .EQ. 0.0D0) GO TO 120
  383. C        ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE
  384. C                  WILL ALWAYS EXIST AN ELEMENT IN THE VECTOR Z(I)
  385. C                  LARGER THAN !!Z(I)!!/N..........
  386.          DO 90 L=1,N
  387.             IF(DABS(Z(L,I)) .GE. NORM(I)/N) GO TO 100
  388.    90       CONTINUE
  389. C
  390.   100    TNORM = DSIGN(NORM(I),Z(L,I))
  391. C
  392.          DO 110 L=1,N
  393.   110       Z(L,I) = Z(L,I)/TNORM
  394. C
  395.          NORM(I) = S/(NORM(I)*NORMA)
  396.   120    RESDUL = DMAX1(NORM(I),RESDUL)
  397. C
  398.       RETURN
  399.       END
  400.       SUBROUTINE RMATIN(NM,MMB,N,MB,ST,STHOLD,INITIL)
  401. C
  402. C     THIS INPUT SUBROUTINE READS A REAL SYMMETRIC BAND MATRIX
  403. C     FROM SYSIN OF ORDER  N,  AND BAND WIDTH  MB .
  404. C     TO GENERATE THE MATRIX  ST  INITIALLY,  INITIL  IS TO BE 0.
  405. C     TO REGENERATE THE MATRIX  ST  FOR THE PURPOSE OF THE RESIDUAL
  406. C     CALCULATION,  INITIL  IS TO BE  1.
  407. C
  408. C     THIS ROUTINE IS CATALOGUED AS  EISPDRV4(RSBREADI).
  409. C
  410.       DOUBLE PRECISION ST(NM,MMB),STHOLD(NM,MMB)
  411.       INTEGER  IA( 5)
  412.       DATA IREADA/1/,IWRITE/6/
  413. C
  414.       IF( INITIL .EQ. 1 )  GO TO  30
  415.       READ(IREADA,5) N,MB
  416.     5 FORMAT(2I6)
  417.       IF( N .EQ. 0 )  GO TO  70
  418.       DO 8 I = 1,N
  419.          DO 7 J = 1,MB
  420.             ST(I,J) = 0.0D0
  421.     7    CONTINUE
  422.     8 CONTINUE
  423.       DO 15 I=1,N
  424.          MBB = MIN0(MB,N-I+1)
  425.          READ(IREADA,10) (IA(J),J=1,MBB)
  426.    10    FORMAT(6I12)
  427.          DO 11 J=1,MBB
  428.            M = MB+1-J
  429.            K = I+J-1
  430.    11      ST(K,M) = IA(J)
  431.    15 CONTINUE
  432.       DO  20  I = 1,N
  433.          DO  20  J = 1,MB
  434.    20      STHOLD(I,J) = ST(I,J)
  435.       RETURN
  436.    30 DO  40  I = 1,N
  437.          DO  40  J = 1,MB
  438.    40      ST(I,J) = STHOLD(I,J)
  439.       RETURN
  440.    70 WRITE(IWRITE,80)
  441.    80 FORMAT(46H0END OF DATA FOR SUBROUTINE  RMATIN(RSBREADI).  /1H1)
  442.       STOP
  443.       END
  444.