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