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