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