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