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 / cgtest.f < prev    next >
Text File  |  1996-09-28  |  21KB  |  587 lines
  1. C
  2. C     THIS DRIVER TESTS  EISPACK  FOR THE CLASS OF COMPLEX GENERAL
  3. C     MATRICES SUMMARIZING THE FIGURES OF MERIT FOR ALL PATHS.
  4. C
  5. C     THIS DRIVER IS CATALOGUED AS  EISPDRV4(CGSUMARY).
  6. C
  7. C     THE DIMENSION OF  AR,AI,ZR,ZI,ASAVER,ASAVEI,RM1,  AND  RM2 SHOULD
  8. C     BE  NM  BY  NM.
  9. C     THE DIMENSION OF  WR,WI,WR1,WI1,SELECT,SLHOLD,INT,SCALE,ORTR,ORTI,
  10. C     RV1  AND  RV2  SHOULD BE  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        ASAVER( 20, 20),ASAVEI( 20, 20),RM1( 20, 20),RM2( 20, 20),
  16.      X        ORTR( 20),ORTI( 20),WR( 20),WI( 20),
  17.      X        WR1( 20),WI1( 20),SCALE( 20),RV1( 20),RV2( 20),
  18.      X        TCRIT( 8),EPSLON,RESDUL,DFL
  19.       DOUBLE PRECISION ARHOLD( 20, 20),AIHOLD( 20, 20)
  20.       INTEGER  INT( 20),IERR( 8),MATCH( 4),LOW,UPP,ERROR,SAME,DIFF
  21.       LOGICAL  SELECT( 20),SLHOLD( 20)
  22.       DATA IREAD1/1/,IREADC/5/,IWRITE/6/
  23.       DATA  SAME,DIFF,LR,QR/4HSAME,4HDIFF,1HL,1HQ/
  24. C
  25.       OPEN(UNIT=IREAD1,FILE='FILE41')
  26.       OPEN(UNIT=IREADC,FILE='FILE42')
  27.       REWIND IREAD1
  28.       REWIND IREADC
  29. C
  30.       NM = 20
  31.       LCOUNT = 0
  32.       WRITE(IWRITE,1)
  33.     1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S
  34.      XTATISTICS//1H ,95(1H-)/ 26X,9HBALANCING,15X,12HNO BALANCING,10X,
  35.      X16HVECTORS SELECTED/ 18X,3(24H------------------------,1X)/
  36.      X16H ORDER LOW UPP T,2X,
  37.      X2(25H  COM-R2 COM-R   CINVIT     )/1H ,95(1H-)//
  38.      X95H 'BALANCING' IS THE OPTION THAT EMPLOYS SUBROUTINE  CBAL  TO  B
  39.      XALANCE THE MATRIX BEFORE THE    /
  40.      X95H EIGENVALUES ARE COMPUTED AND  CBABK2  TO BACK TRANSFORM THE EI
  41.      XGENVECTORS.                     //
  42.      X95H 'VECTORS COMPUTED' IDENTIFIES BY POSITION, IN THE SET RETURNED
  43.      X BY  COM-R, THOSE EIGENVALUES   /
  44.      X95H FOR WHICH  CINVIT  COMPUTED THE ASSOCIATED EIGENVECTORS.  T  I
  45.      XNDEXES AN EIGENVALUE FOR WHICH  /
  46.      X66H THE EIGENVECTOR WAS COMPUTED AND  F  INDEXES A PASSED EIGENVAL
  47.      XUE.    // 51H UNDER 'ORDER' IS THE ORDER OF EACH TEST MATRIX.   //
  48.      X95H UNDER 'LOW' AND 'UPP' ARE INTEGERS INDICATING THE BOUNDARY IND
  49.      XICES FOR THE BALANCED MATRIX.   //
  50.      X95H UNDER 'T' IS THE LETTER L OR Q INDICATING THE USE OF ELEMENTAR
  51.      XY OR UNITARY TRANSFORMATIONS.   /)
  52.       WRITE(IWRITE,2)
  53.     2 FORMAT(
  54.      X61H PERFORMANCE REPORTED UNDER 'COM-R2 COM-R' IS THAT OF  COMLR2 ,
  55.      X31H  AND  COMLR  ON LINE HEADED BY   /
  56.      X60H 'L' AND THAT OF  COMQR2  AND  COMQR  ON LINE HEADED BY 'Q'.//
  57.      X95H UNDER 'COM-R2 COM-R' ARE TWO NUMBERS AND A KEYWORD.  THE FIRST
  58.      X NUMBER, AN INTEGER, IS THE     /
  59.      X95H ABSOLUTE SUM OF THE ERROR FLAGS RETURNED SEPARATELY FROM  COM-
  60.      XR2  AND  COM-R.  THE SECOND     /
  61.      X95H NUMBER IS THE MEASURE OF PERFORMANCE BASED UPON THE RESIDUAL C
  62.      XOMPUTED FOR THE  COM-R2  PATH.  /
  63.      X95H THE KEYWORD REPORTS THE DUPLICATION OF THE COMPUTED EIGENVALUE
  64.      XS FROM  COM-R2  AND  COM-R.     /
  65.      X95H 'SAME' MEANS THAT THE EIGENVALUES ARE EXACT DUPLICATES.  'DIFF
  66.      X' MEANS THAT FOR AT LEAST ONE   /
  67.      X95H PAIR OF CORRESPONDING EIGENVALUES, THE MEMBERS OF THE PAIR ARE
  68.      X NOT IDENTICAL.                 //
  69.      X95H UNDER 'CINVIT' ARE TWO NUMBERS.  THE FIRST NUMBER, AN INTEGER,
  70.      X IS THE ABSOLUTE SUM OF THE     /
  71.      X62H ERROR FLAGS RETURNED FROM THE PATH.  THE SECOND NUMBER IS THE,
  72.      X29H MEASURE OF PERFORMANCE BASED    )
  73.       WRITE(IWRITE,3)
  74.     3 FORMAT(
  75.      X41H UPON THE RESIDUAL COMPUTED FOR THE PATH.//
  76.      X62H -1.0  AS THE MEASURE OF PERFORMANCE IS PRINTED IF AN ERROR IN,
  77.      X27H THE CORRESPONDING PATH HAS      /
  78.      X47H PREVENTED THE COMPUTATION OF THE EIGENVECTORS.//
  79.      X45H0FOR REDUCTIONS BY ELEMENTARY TRANSFORMATIONS/
  80.      X95H THE  COMLR2   PATH WITH    BALANCING USES THE EISPACK CODES  C
  81.      XBAL  -COMHES-COMLR2-CBABK2.     /
  82.      X95H THE   COMLR   PATH WITH    BALANCING USES THE EISPACK CODES  C
  83.      XBAL  -COMHES-COMLR .            /
  84.      X60H THE  CINVIT   PATH WITH    BALANCING USES THE EISPACK CODES ,
  85.      X44H  CBAL  -COMHES-COMLR -CINVIT-COMBAK-CBABK2. /
  86.      X95H THE  COMLR2   PATH WITH NO BALANCING USES THE EISPACK CODES  C
  87.      XOMHES-COMLR2.                   /
  88.      X95H THE   COMLR   PATH WITH NO BALANCING USES THE EISPACK CODES  C
  89.      XOMHES-COMLR .                   /
  90.      X95H THE  CINVIT   PATH WITH NO BALANCING USES THE EISPACK CODES  C
  91.      XOMHES-COMLR -CINVIT-COMBAK.     )
  92.       WRITE(IWRITE,4)
  93.     4 FORMAT(42H0FOR REDUCTIONS BY UNITARY TRANSFORMATIONS/
  94.      X95H THE  COMQR2   PATH WITH    BALANCING USES THE EISPACK CODES  C
  95.      XBAL  -CORTH -COMQR2-CBABK2,           /
  96.      X38H AS CALLED FROM DRIVER SUBROUTINE  CG. /
  97.      X95H THE   COMQR   PATH WITH    BALANCING USES THE EISPACK CODES  C
  98.      XBAL  -CORTH -COMQR ,            /
  99.      X38H AS CALLED FROM DRIVER SUBROUTINE  CG. /
  100.      X60H THE  CINVIT   PATH WITH    BALANCING USES THE EISPACK CODES ,
  101.      X44H  CBAL  -CORTH -COMQR -CINVIT-CORTB -CBABK2. /
  102.      X95H THE  COMQR2   PATH WITH NO BALANCING USES THE EISPACK CODES  C
  103.      XORTH -COMQR2.                   /
  104.      X95H THE   COMQR   PATH WITH NO BALANCING USES THE EISPACK CODES  C
  105.      XORTH -COMQR .                   /
  106.      X95H THE  CINVIT   PATH WITH NO BALANCING USES THE EISPACK CODES  C
  107.      XORTH -COMQR -CINVIT-CORTB .     )
  108.       WRITE(IWRITE,15)
  109.    15 FORMAT(1X,21HD.P. VERSION 04/15/83 )
  110.     5 FORMAT( 53H1       TABULATION OF THE ERROR FLAG  ERROR  AND THE ,
  111.      X    31HMEASURE OF PERFORMANCE  Y  FOR /5X,
  112.      X    56HTHE  EISPACK  CODES.  THIS RUN DISPLAYS THESE STATISTICS ,
  113.      X    33H FOR COMPLEX GENERAL MATRICES.      //
  114.      X    26X,9HBALANCING,15X,12HNO BALANCING,10X,16HVECTORS SELECTED/
  115.      X    18X,3(24H------------------------,1X)/
  116.      X    16H ORDER LOW UPP T,2X,2(25H  COM-R2 COM-R   CINVIT       ))
  117.    10 CALL CMATIN(NM,N,AR,AI,ARHOLD,AIHOLD,0)
  118.       READ(IREADC,50) MM,(SLHOLD(I),I=1,N)  
  119.    50 FORMAT(I4/(72L1))
  120. C
  121. C     MM  AND  SELECT  ARE READ FROM SYSIN AFTER THE MATRIX IS
  122. C     GENERATED.  MM  SPECIFIES TO  INVIT  THE MAXIMUM -UMBER OF
  123. C     EIGENVECTORS THAT WILL BE "OMPUTED.  SELECT  CONTAINS INFORMATION
  124. C     ABOUT WHICH EIGENVECTORS ARE DESIRED.
  125. C
  126.       DO  300  ICALL = 1,12
  127.          ICALL2 = MOD(ICALL,3)
  128.          IF( ICALL2 .NE. 0 ) GO TO 80
  129.          DO  70  I = 1,N
  130.    70    SELECT(I) = SLHOLD(I)
  131.    80    IF( ICALL .NE. 1 )  CALL  CMATIN(NM,N,AR,AI,ARHOLD,AIHOLD,1)
  132.          GO TO  (90,100,110,120,130,140,160,165,170,175,180,185), ICALL
  133. C
  134. C     IF  COM-R  PATH IS TAKEN THEN  COM-R2  PATH MUST ALSO BE TAKEN
  135. C     IN ORDER THAT COMPARISON OF EIGENVALUES BE POSSIBLE.
  136. C
  137. C
  138. C     CGEWZ
  139. C
  140.    90    ICT = 1
  141.          CALL  CBAL(NM,N,AR,AI,LOW,UPP,SCALE)
  142.          CALL  COMHES(NM,N,LOW,UPP,AR,AI,INT)
  143.          CALL  COMLR2(NM,N,LOW,UPP,INT,AR,AI,WR,WI,ZR,ZI,ERROR)
  144.          IERR(ICT) = ERROR
  145.          IF( ERROR .NE. 0 ) GO TO 220
  146.          DO  91  I = 1,N
  147.             WR1(I) = WR(I)
  148.    91       WI1(I) = WI(I)
  149.          CALL  CBABK2(NM,N,LOW,UPP,SCALE,N,ZR,ZI)
  150.          GO TO  190
  151. C
  152. C     CGEW
  153. C
  154.   100    IF( IERR(1) .NE. 0 ) GO TO 300
  155.          CALL  CBAL(NM,N,AR,AI,LOW,UPP,SCALE)
  156.          CALL  COMHES(NM,N,LOW,UPP,AR,AI,INT)
  157.          CALL  COMLR(NM,N,LOW,UPP,AR,AI,WR,WI,ERROR)
  158.          IERR(1) = ERROR
  159.          MATCH(1) = DIFF
  160.          DO  101  I = 1,N
  161.             IF( WR(I) .NE. WR1(I) .OR. WI(I) .NE. WI1(I) )  GO TO  300
  162.   101    CONTINUE
  163.          MATCH(1) = SAME
  164.          GO TO  300
  165. C
  166. C     CGEW1Z
  167. C
  168.   110    ICT = 2
  169.          CALL  CBAL(NM,N,AR,AI,LOW,UPP,SCALE)
  170.          CALL  COMHES(NM,N,LOW,UPP,AR,AI,INT)
  171.          DO  112  I = 1,N
  172.             DO  112  J = 1,N
  173.               ASAVER(I,J) = AR(I,J)
  174.   112         ASAVEI(I,J) = AI(I,J)
  175.          CALL  COMLR(NM,N,LOW,UPP,AR,AI,WR,WI,ERROR)
  176.          IERR(ICT) = ERROR
  177.          IF( ERROR .NE. 0 )  GO TO  220
  178.          CALL  CINVIT(NM,N,ASAVER,ASAVEI,WR,WI,SELECT,MM,M,ZR,ZI,
  179.      X                ERROR,RM1,RM2,RV1,RV2)
  180.          IERR(ICT) = IABS(ERROR)
  181.          CALL  COMBAK(NM,LOW,UPP,ASAVER,ASAVEI,INT,M,ZR,ZI)
  182.          CALL  CBABK2(NM,N,LOW,UPP,SCALE,M,ZR,ZI)
  183.          GO TO  190
  184. C
  185. C     CGNEWZ
  186. C
  187.   120    ICT = 3
  188.          CALL  COMHES(NM,N,1,N,AR,AI,INT)
  189.          CALL  COMLR2(NM,N,1,N,INT,AR,AI,WR,WI,ZR,ZI,ERROR)
  190.          IERR(ICT) = ERROR
  191.          IF( ERROR .NE. 0 ) GO TO 220
  192.          DO  121  I = 1,N
  193.             WR1(I) = WR(I)
  194.   121       WI1(I) = WI(I)
  195.          GO TO  190
  196. C
  197. C     CGNEW
  198. C
  199.   130    IF( IERR(3) .NE. 0 ) GO TO 300
  200.          CALL  COMHES(NM,N,1,N,AR,AI,INT)
  201.          CALL  COMLR(NM,N,1,N,AR,AI,WR,WI,ERROR)
  202.          IERR(3) = ERROR
  203.          MATCH(2) = DIFF
  204.          DO  131  I = 1,N
  205.             IF( WR(I) .NE. WR1(I) .OR. WI(I) .NE. WI1(I) )  GO TO  300
  206.   131    CONTINUE
  207.          MATCH(2) = SAME
  208.          GO TO  300
  209. C
  210. C     CGNEW1Z
  211. C
  212.   140    ICT = 4
  213.          CALL  COMHES(NM,N,1,N,AR,AI,INT)
  214.          DO  142  I = 1,N
  215.             DO  142  J = 1,N
  216.               ASAVER(I,J) = AR(I,J)
  217.   142         ASAVEI(I,J) = AI(I,J)
  218.          CALL  COMLR(NM,N,1,N,AR,AI,WR,WI,ERROR)
  219.          IERR(ICT) = ERROR
  220.          IF( ERROR .NE. 0 )  GO TO  220
  221.          CALL  CINVIT(NM,N,ASAVER,ASAVEI,WR,WI,SELECT,MM,M,ZR,ZI,
  222.      X                ERROR,RM1,RM2,RV1,RV2)
  223.          IERR(ICT) = IABS(ERROR)
  224.          CALL  COMBAK(NM,1,N,ASAVER,ASAVEI,INT,M,ZR,ZI)
  225.          GO TO 190
  226. C
  227. C     CGUWZ
  228. C     INVOKED FROM DRIVER SUBROUTINE  CG.
  229. C
  230.   160    ICT = 5
  231.          CALL  CG(NM,N,AR,AI,WR,WI,1,ZR,ZI,SCALE,ORTR,ORTI,ERROR)
  232.          IERR(ICT) = ERROR
  233.          IF( ERROR .NE. 0 ) GO TO 220
  234.          GO TO  190
  235. C
  236. C     CGUW
  237. C     INVOKED FROM DRIVER SUBROUTINE  CG.
  238. C
  239.   165    IF( IERR(5) .NE. 0 ) GO TO 300
  240.          CALL  CG(NM,N,AR,AI,WR1,WI1,0,AR,AI,SCALE,ORTR,ORTI,ERROR)
  241.          IERR(5) = ERROR
  242.          MATCH(3) = DIFF
  243.          DO  166  I = 1,N
  244.             IF( WR(I) .NE. WR1(I) .OR. WI(I) .NE. WI1(I) )  GO TO  300
  245.   166    CONTINUE
  246.          MATCH(3) = SAME
  247.          GO TO  300
  248. C
  249. C     CGUW1Z
  250. C
  251.   170    ICT = 6
  252.          CALL  CBAL(NM,N,AR,AI,LOW,UPP,SCALE)
  253.          CALL  CORTH(NM,N,LOW,UPP,AR,AI,ORTR,ORTI)
  254.          DO  172  I = 1,N
  255.             DO  172  J = 1,N
  256.               ASAVER(I,J) = AR(I,J)
  257.   172         ASAVEI(I,J) = AI(I,J)
  258.          CALL  COMQR(NM,N,LOW,UPP,AR,AI,WR,WI,ERROR)
  259.          IERR(ICT) = ERROR
  260.          IF( ERROR .NE. 0 )  GO TO  220
  261.          CALL  CINVIT(NM,N,ASAVER,ASAVEI,WR,WI,SELECT,MM,M,ZR,ZI,
  262.      X                        ERROR,RM1,RM2,RV1,RV2)
  263.          IERR(ICT) = IABS(ERROR)
  264.          CALL  CORTB(NM,LOW,UPP,ASAVER,ASAVEI,ORTR,ORTI,M,ZR,ZI)
  265.          CALL  CBABK2(NM,N,LOW,UPP,SCALE,M,ZR,ZI)
  266.          GO TO  190
  267. C
  268. C     CGNUWZ
  269. C
  270.   175    ICT = 7
  271.          CALL  CORTH(NM,N,1,N,AR,AI,ORTR,ORTI)
  272.          CALL  COMQR2(NM,N,1,N,ORTR,ORTI,AR,AI,WR,WI,ZR,ZI,ERROR)
  273.          IERR(ICT) = ERROR
  274.          IF( ERROR .NE. 0 ) GO TO 220
  275.          DO  176  I = 1,N
  276.             WR1(I) = WR(I)
  277.   176       WI1(I) = WI(I)
  278.          GO TO  190
  279. C
  280. C     CGNUW
  281. C
  282.   180    IF( IERR(7) .NE. 0 ) GO TO 300
  283.          CALL  CORTH(NM,N,1,N,AR,AI,ORTR,ORTI)
  284.          CALL  COMQR(NM,N,1,N,AR,AI,WR,WI,ERROR)
  285.          IERR(7) = ERROR
  286.          MATCH(4) = DIFF
  287.          DO  181  I = 1,N
  288.             IF( WR(I) .NE. WR1(I) .OR. WI(I) .NE. WI1(I) )  GO TO  300
  289.   181    CONTINUE
  290.          MATCH(4) = SAME
  291.          GO TO  300
  292. C
  293. C     CGNUW1Z
  294. C
  295.   185    ICT = 8
  296.          CALL  CORTH(NM,N,1,N,AR,AI,ORTR,ORTI)
  297.          DO  187  I = 1,N
  298.             DO  187  J = 1,N
  299.               ASAVER(I,J) = AR(I,J)
  300.   187         ASAVEI(I,J) = AI(I,J)
  301.          CALL  COMQR(NM,N,1,N,AR,AI,WR,WI,ERROR)
  302.          IERR(ICT) = ERROR
  303.          IF( ERROR .NE. 0 )  GO TO  220
  304.          CALL  CINVIT(NM,N,ASAVER,ASAVEI,WR,WI,SELECT,MM,M,ZR,ZI,
  305.      X                ERROR,RM1,RM2,RV1,RV2)
  306.          IERR(ICT) = IABS(ERROR)
  307.          CALL  CORTB(NM,1,N,ASAVER,ASAVEI,ORTR,ORTI,M,ZR,ZI)
  308. C
  309.   190    CALL  CMATIN(NM,N,AR,AI,ARHOLD,AIHOLD,1)
  310.          IF( ICALL2 .EQ. 0 ) CALL CGW1ZR(NM,N,AR,AI,WR,WI,SELECT,
  311.      X                                   M,ZR,ZI,RV1,RESDUL)
  312.          IF( ICALL2 .EQ. 1 ) CALL CGWZR(NM,N,AR,AI,WR,WI,ZR,ZI,
  313.      X                                  RV1,RESDUL)
  314.          DFL = 10 * N
  315.          TCRIT(ICT) = RESDUL/EPSLON(DFL)
  316.          GO TO 300
  317.   220    TCRIT(ICT) = -1.0D0
  318.   300 CONTINUE
  319. C
  320.       IF( MOD(LCOUNT,16) .EQ. 0 ) WRITE(IWRITE,5)
  321.       LCOUNT = LCOUNT + 1
  322.       WRITE(IWRITE,520) N,LOW,UPP,LR,(IERR(2*I-1),TCRIT(2*I-1),MATCH(I),
  323.      X              IERR(2*I),TCRIT(2*I),I=1,2),(SELECT(I),I=1,N)
  324.   520 FORMAT(I4,1X,2I4,2X,A1,2(I4,F6.3,1X,A4,I4,F6.3),
  325.      X       4X,21L1/(72X,20L1))
  326.       WRITE(IWRITE,530) QR,(IERR(2*I-1),TCRIT(2*I-1),MATCH(I),
  327.      X              IERR(2*I),TCRIT(2*I),I=3,4)
  328.   530 FORMAT(15X,A1,2(I4,F6.3,1X,A4,I4,F6.3))
  329.       GO TO  10
  330.       END
  331.       SUBROUTINE CMATIN(NM,N,AR,AI,ARHOLD,AIHOLD,INITIL)
  332. C
  333. C     THIS INPUT SUBROUTINE READS A COMPLEX MATRIX  A = (AR,AI)
  334. C     FROM SYSIN OF ORDER N.
  335. C     TO GENERATE THE MATRIX  A  INITIALLY,  INITIL  IS TO BE 0.
  336. C     TO REGENERATE THE MATRIX  A  FOR THE PURPOSE OF THE RESIDUAL
  337. C     CALCULATION,  INITIL  IS TO BE  1.
  338. C
  339. C     THIS ROUTINE IS CATALOGUED AS  EISPDRV4(CGREADI).
  340. C
  341.       DOUBLE PRECISION AR(NM,NM),AI(NM,NM),ARHOLD(NM,NM),AIHOLD(NM,NM)
  342.       INTEGER  IAR( 20), IAI( 20)
  343.       DATA IREADA/1/,IWRITE/6/
  344. C
  345.       IF( INITIL .EQ. 1 )  GO TO  30
  346.       READ(IREADA,5) N
  347.     5 FORMAT(I6)
  348.       IF( N .EQ. 0 )  GO TO  70
  349.       DO  15  I = 1,N
  350.          READ(IREADA,10) (IAR(J),IAI(J),J=1,N)
  351.    10    FORMAT(2I18)
  352.          DO  15  J = 1,N
  353.            AR(I,J) = IAR(J)
  354.    15      AI(I,J) = IAI(J)
  355.       DO  20  I = 1,N
  356.          DO  20  J = 1,N
  357.            ARHOLD(I,J) = AR(I,J)
  358.    20      AIHOLD(I,J) = AI(I,J)
  359.       RETURN
  360.    30 DO  40  I = 1,N
  361.          DO  40  J = 1,N
  362.            AR(I,J) = ARHOLD(I,J)
  363.    40      AI(I,J) = AIHOLD(I,J)
  364.       RETURN
  365.    70 WRITE(IWRITE,80)
  366.    80 FORMAT(44H0END OF DATA FOR SUBROUTINE CMATIN(CGREADI).  /1H1)
  367.       STOP
  368.       END
  369.       SUBROUTINE CGWZR(NM,N,AR,AI,WR,WI,ZR,ZI,NORM,RESDUL)
  370. C
  371.       DOUBLE PRECISION NORM(N),WR(N),WI(N),AR(NM,N),AI(NM,N),
  372.      X       ZR(NM,N), ZI(NM,N), NORMA, XR, XI, S, SUMA, SUMZ,
  373.      X       SUMR, SUMI, RESDUL, PYTHAG
  374. C
  375. C     THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
  376. C     A*Z-Z*DIAG(W)  WHERE  A  IS A COMPLEX GENERAL MATRIX,  Z  IS A
  377. C     MATRIX WHOSE COLUMNS ARE APPROXIMATE EIGENVECTORS OF  A,
  378. C     AND  W  IS A VECTOR WHOSE ELEMENTS ARE APPROXIMATE EIGENVALUES
  379. C     CORRESPONDING TO THE VECTORS IN  Z.  ALL NORMS USED IN THE
  380. C     COMMENTS BELOW ARE 1-NORMS.
  381. C
  382. C     THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(CGWZR)
  383. C
  384. C     INPUT.
  385. C
  386. C        NM IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
  387. C           AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT.
  388. C
  389. C        N IS THE ORDER OF THE MATRIX  A;
  390. C
  391. C        AR(NM,N), AI(NM,N) ARE ARRAYS CONTAINING THE REAL AND
  392. C           IMAGINARY PARTS OF THE ELEMENTS OF  A;
  393. C
  394. C        ZR(NM,N), ZI(NM,N) ARE ARRAYS CONTAINING THE REAL AND
  395. C           IMAGINARY PARTS OF THE ELEMENTS OF  Z;
  396. C
  397. C        WR(N), WI(N) ARE ARRAYS CONTAINING THE REAL AND IMAGINARY
  398. C           PARTS OF  W.
  399. C
  400. C     OUTPUT.
  401. C
  402. C        ZR(NM,N), ZI(NM,N) ARE ARRAYS WHOSE COLUMNS CONTAIN THE
  403. C           REAL AND IMAGINARY PARTS OF THE NORMALIZED APPROXIMATE
  404. C           EIGENVECTORS OF  A.  THE EIGENVECTORS ARE NORMALIZED BY
  405. C           THE 1-NORM IN SUCH A WAY THAT THE FIRST ELEMENT WHOSE
  406. C           MAGNITUDE IS LARGER THAN THE NORM OF THE EIGENVECTOR
  407. C           DIVIDED BY  N  IS REAL AND POSITIVE;
  408. C
  409. C        NORM(N) IS AN ARRAY SUCH THAT FOR EACH  K,
  410. C           NORM(K) = !!A*Z(K)-Z(K)*(W(K))!!/(!!A!!*!!Z(K)!!)
  411. C           WHERE  Z(K)  IS THE K-TH EIGENVECTOR;
  412. C
  413. C        RESDUL IS THE REAL NUMBER
  414. C           !!A*Z-Z*DIAG(W)!!/(!!A!!*!!Z!!).
  415. C
  416. C     ----------------------------------------------------------------
  417. C
  418.       RESDUL = 0.0D0
  419.       NORMA = 0.0D0
  420. C
  421.       DO 20 I=1,N
  422.          SUMA = 0.0D0
  423. C
  424.          DO 10 L=1,N
  425.    10       SUMA = SUMA + PYTHAG(AR(L,I),AI(L,I))
  426. C
  427.    20    NORMA = DMAX1(NORMA,SUMA)
  428. C
  429.       IF(NORMA .EQ. 0.0D0) NORMA = 1.0D0
  430. C
  431.       DO 80 I=1,N
  432.          SUMZ = 0.0D0
  433.          S = 0.0D0
  434. C
  435.          DO 40 L=1,N
  436.             SUMZ = SUMZ + PYTHAG(ZR(L,I),ZI(L,I))
  437.             SUMR = -WR(I)*ZR(L,I) + WI(I)*ZI(L,I)
  438.             SUMI = -WR(I)*ZI(L,I) - WI(I)*ZR(L,I)
  439. C
  440.             DO 30 K=1,N
  441.                SUMR = SUMR + AR(L,K)*ZR(K,I) - AI(L,K)*ZI(K,I)
  442.    30          SUMI = SUMI + AR(L,K)*ZI(K,I) + AI(L,K)*ZR(K,I)
  443. C
  444.    40       S = S + PYTHAG(SUMR,SUMI)
  445. C
  446.          NORM(I) = SUMZ
  447. C        ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE
  448. C                  WILL ALWAYS EXIST AN ELEMENT IN THE VECTOR Z(I)
  449. C                  LARGER THAN !!Z(I)!!/N..........
  450.          DO 50 L=1,N
  451.             IF(PYTHAG(ZR(L,I),ZI(L,I)) .GE. NORM(I)/N) GO
  452.      1      TO 60
  453.    50        CONTINUE
  454. C
  455.    60    XR = NORM(I)*ZR(L,I)/PYTHAG(ZR(L,I),ZI(L,I))
  456.          XI = NORM(I)*ZI(L,I)/PYTHAG(ZR(L,I),ZI(L,I))
  457. C
  458.          DO 70 L=1,N
  459.             CALL CDIV(ZR(L,I),ZI(L,I),XR,XI,ZR(L,I),ZI(L,I))
  460.    70    CONTINUE
  461. C
  462.          NORM(I) = S/(NORM(I)*NORMA)
  463.          RESDUL = DMAX1(NORM(I),RESDUL)
  464.    80    CONTINUE
  465. C
  466.       RETURN
  467.       END
  468.       SUBROUTINE CGW1ZR(NM,N,AR,AI,WR,WI,SELECT,M,ZR,ZI,NORM,RESDUL)
  469. C
  470.       DOUBLE PRECISION NORM(N),WR(N),WI(N),AR(NM,N),AI(NM,N),
  471.      X       ZR(NM,M), ZI(NM,M), NORMA, XR, XI, S, SUMA, SUMZ,
  472.      X       SUMR, SUMI, RESDUL, PYTHAG
  473.       LOGICAL SELECT(N)
  474. C
  475. C     THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
  476. C     A*Z-Z*DIAG(W)   WHERE  A  IS A COMPLEX GENERAL MATRIX,  Z  IS A
  477. C     MATRIX WHOSE COLUMNS ARE APPROXIMATE EIGENVECTORS OF  A,
  478. C     AND  W  IS A VECTOR WHOSE ELEMENTS ARE APPROXIMATE EIGENVALUES
  479. C     CORRESPONDING TO THE VECTORS IN  Z.  ALL NORMS USED IN THE
  480. C     COMMENTS BELOW ARE 1-NORMS.
  481. C
  482. C     THIS SUBROUTINE IS CATALOGUED AS  EISPDRV4(CGW1ZR).
  483. C
  484. C     INPUT.
  485. C
  486. C        NM  IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
  487. C           AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT;
  488. C
  489. C        N  IS THE ORDER OF THE MATRIX A;
  490. C
  491. C        AR(NM,N), AI(NM,N) ARE ARRAYS CONTAINING THE REAL AND
  492. C           IMAGINARY PARTS OF THE ELEMENTS OF  A;
  493. C
  494. C        ZR(NM,M), ZI(NM,M) ARE ARRAYS CONTAINING THE REAL AND
  495. C           IMAGINARY PARTS OF THE ELEMENTS OF  Z;
  496. C
  497. C        WR(N), WI(N) ARE ARRAYS CONTAINING THE REAL AND IMAGINARY
  498. C           PARTS OF  W;
  499. C
  500. C        SELECT(N)  IS A BOOLEAN ARRAY WHICH SPECIFIES THE EIGENVALUES
  501. C           IN THE VECTOR  W  THAT ARE ASSOCIATED WITH EIGENVECTORS
  502. C           IN THE MATRIX  Z.  THE I-TH VALUE IN  W  IS ASSOCIATED WITH
  503. C           A VECTOR IN  Z  IF SELECT(I) = .TRUE..  FURTHERMORE IF
  504. C           SELECT(I) IS THE J-TH  .TRUE.  ELEMENT, THEN THE I-TH
  505. C           ELEMENT OF  W  IS ASSOCIATED WITH THE J-TH VECTOR OF  Z;
  506. C
  507. C        M  IS THE NUMBER OF ELEMENTS IN SELECT WHOSE VALUE IS  .TRUE..
  508. C
  509. C     OUTPUT.
  510. C
  511. C        ZR(NM,M), ZI(NM,M) ARE ARRAYS WHOSE COLUMNS CONTAIN THE
  512. C           REAL AND IMAGINARY PARTS OF THE NORMALIZED APPROXIMATE
  513. C           EIGENVECTORS OF  A.  THE EIGENVECTORS ARE NORMALIZED BY
  514. C           THE 1-NORM IN SUCH A WAY THAT THE FIRST ELEMENT WHOSE
  515. C           MAGNITUDE IS LARGER THAN THE NORM OF THE EIGENVECTOR
  516. C           DIVIDED BY  N  IS REAL AND POSITIVE;
  517. C
  518. C        NORM(N)  IS AN ARRAY SUCH THAT WHENEVER SELECT(K) =  .TRUE.
  519. C           NORM(K) = !!A*Z(K)-Z(K)*W(K)!!/(!!A!!*!!Z(K)!!)
  520. C           WHERE  Z(K)  IS THE EIGENVECTOR CORRESPONDING TO THE
  521. C           K-TH EIGENVALUE (WR(K),WI(K));
  522. C
  523. C        RESDUL  IS A REAL NUMBER DEFINED BY THE FOLLOWING.
  524. C           !!A*Z-Z*DIAG(W)!!/(!!A!!*!!Z!!).
  525. C
  526. C     ----------------------------------------------------------------
  527. C
  528.       RESDUL = 0.0D0
  529.       IF( M .EQ. 0 ) RETURN
  530.       NORMA = 0.0D0
  531.       II=0
  532. C
  533.       DO 20 I=1,N
  534.          SUMA = 0.0D0
  535. C
  536.          DO 10 L=1,N
  537.    10       SUMA = SUMA + PYTHAG(AR(L,I),AI(L,I))
  538. C
  539.          NORMA = DMAX1(NORMA,SUMA)
  540.          IF(NORMA .EQ. 0.0D0) NORMA = 1.0D0
  541.    20    CONTINUE
  542. C
  543.       IF(NORMA .EQ. 0.0D0) NORMA = 1.0D0
  544. C
  545.       DO 80 I=1,N
  546.          SUMZ = 0.0D0
  547.          S = 0.0D0
  548.          IF( .NOT. SELECT(I) ) GO TO 80
  549.          II = II+1
  550. C
  551.          DO 40 L=1,N
  552.             SUMZ = SUMZ + PYTHAG(ZR(L,II),ZI(L,II))
  553.             SUMR = -WR(I)*ZR(L,II) + WI(I)*ZI(L,II)
  554.             SUMI = -WR(I)*ZI(L,II) - WI(I)*ZR(L,II)
  555. C
  556.             DO 30 K=1,N
  557.                SUMR = SUMR + AR(L,K)*ZR(K,II) - AI(L,K)*ZI(K,II)
  558.    30          SUMI = SUMI + AR(L,K)*ZI(K,II) + AI(L,K)*ZR(K,II)
  559. C
  560.    40       S = S + PYTHAG(SUMR,SUMI)
  561. C
  562.          IF(SUMZ .NE. 0.0D0) GO TO 45
  563.          NORM(I) = -1
  564.          GO TO 80
  565.    45    NORM(I) = SUMZ
  566. C        ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE
  567. C                  WILL ALWAYS EXIST AN ELEMENT IN THE VECTOR  Z(II)
  568. C                  LARGER THAN !!Z(II)!!/N..........
  569.          DO 50 L=1,N
  570.             IF(PYTHAG(ZR(L,II),ZI(L,II)) .GE. NORM(I)/N)
  571.      1      GO TO 60
  572.    50        CONTINUE
  573. C
  574.    60    XR = NORM(I)*ZR(L,II)/PYTHAG(ZR(L,II),ZI(L,II))
  575.          XI = NORM(I)*ZI(L,II)/PYTHAG(ZR(L,II),ZI(L,II))
  576. C
  577.          DO 70 L=1,N
  578.             CALL CDIV(ZR(L,II),ZI(L,II),XR,XI,ZR(L,II),ZI(L,II))
  579.    70    CONTINUE
  580. C
  581.          NORM(I) = S/(NORM(I)*NORMA)
  582.          RESDUL = DMAX1(NORM(I),RESDUL)
  583.    80    CONTINUE
  584. C
  585.       RETURN
  586.       END
  587.