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