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 / rgtest.f < prev    next >
Text File  |  1996-09-28  |  23KB  |  658 lines
  1. C
  2. C     THIS DRIVER TESTS  EISPACK  FOR THE CLASS OF REAL GENERAL
  3. C     MATRICES SUMMARIZING THE FIGURES OF MERIT FOR ALL PATHS.
  4. C
  5. C     THIS DRIVER IS CATALOGUED AS  EISPDRV4(RGSUMARY).
  6. C
  7. C     THE DIMENSION OF  A,Z,RM1, AND  ASAVE  SHOULD BE  NM  BY  NM.
  8. C     THE DIMENSION OF  WR,WI,WR1,WI1,SELECT,SLHOLD,INT,SCALE,ORT,RV1,
  9. C     AND  RV2  SHOULD BE  NM.
  10. C     THE DIMENSION OF  AHOLD  SHOULD BE  NM  BY  NM.
  11. C     HERE NM = 20.
  12. C
  13.       DOUBLE PRECISION A(20,20),Z(20,20),ASAVE(20,20),RM1(20,20),
  14.      X        AHOLD( 20, 20),ORT( 20),WR( 20),WI(20),
  15.      X        WR1( 20),WI1( 20),SCALE( 20),RV1( 20),RV2( 20),
  16.      X        TCRIT( 8),EPSLON,RESDUL,DFL
  17.       INTEGER  INT( 20),IERR( 8),MATCH( 4),LOW,UPP,ERROR,SAME,DIFF,
  18.      X         ELEM,ORTH
  19.       LOGICAL  SELECT( 20),SLHOLD( 20)
  20.       DATA IREAD1/1/,IREADC/5/,IWRITE/6/
  21.       DATA  SAME,DIFF,ELEM,ORTH/4HSAME,4HDIFF,1HE,1HO/
  22. C
  23.       OPEN(UNIT=IREAD1,FILE='FILE33')
  24.       OPEN(UNIT=IREADC,FILE='FILE34')
  25.       REWIND IREAD1
  26.       REWIND IREADC
  27. C
  28.       NM = 20
  29.       LCOUNT = 0
  30.       WRITE(IWRITE,1)
  31.     1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S
  32.      XTATISTICS//1H ,95(1H-)/ 26X,9HBALANCING,15X,12HNO BALANCING,10X,
  33.      X16HVECTORS COMPUTED/ 18X,3(24H------------------------,1X)/
  34.      X16H ORDER LOW UPP T,2X,2(25H   HQR2   HQR     INVIT      )/
  35.      X1H ,95(1H-)//95H 'BALANCING' IS THE OPTION THAT EMPLOYS SUBROUTINE
  36.      X  BALANC  TO BALANCE THE MATRIX BEFORE THE   /
  37.      X95H EIGENVALUES ARE COMPUTED AND  BALBAK  TO BACK TRANSFORM THE EI
  38.      XGENVECTORS.                     //
  39.      X95H 'VECTORS COMPUTED' IDENTIFIES BY POSITION, IN THE SET RETURNED
  40.      X BY  HQR, THOSE EIGENVALUES     /
  41.      X95H FOR WHICH  INVIT  COMPUTED THE ASSOCIATED EIGENVECTORS.  T  IN
  42.      XDEXES AN EIGENVALUE FOR WHICH   /
  43.      X66H THE EIGENVECTOR WAS COMPUTED AND  F  INDEXES A PASSED EIGENVAL
  44.      XUE.    // 51H UNDER 'ORDER' IS THE ORDER OF EACH TEST MATRIX.   //
  45.      X95H UNDER 'LOW' AND 'UPP' ARE INTEGERS INDICATING THE BOUNDARY IND
  46.      XICES FOR THE BALANCED MATRIX.   //
  47.      X95H UNDER 'T' IS THE LETTER E OR O INDICATING THE USE OF ELEMENTAR
  48.      XY OR ORTHOGONAL TRANSFORMATIONS.//
  49.      X95H UNDER 'HQR2   HQR' ARE TWO NUMBERS AND A KEYWORD.  THE FIRST N
  50.      XUMBER, AN INTEGER, IS THE       )
  51.       WRITE(IWRITE,2)
  52.     2 FORMAT(
  53.      X95H ABSOLUTE SUM OF THE ERROR FLAGS RETURNED SEPARATELY FROM  HQR2
  54.      X  AND  HQR.  THE SECOND         /
  55.      X95H NUMBER IS THE MEASURE OF PERFORMANCE BASED UPON THE RESIDUAL C
  56.      XOMPUTED FOR THE  HQR2  PATH.    /
  57.      X95H THE KEYWORD REPORTS THE DUPLICATION OF THE COMPUTED EIGENVALUE
  58.      XS FROM  HQR2  AND  HQR.         /
  59.      X95H 'SAME' MEANS THAT THE EIGENVALUES ARE EXACT DUPLICATES.  'DIFF
  60.      X' MEANS THAT FOR AT LEAST ONE   /
  61.      X95H PAIR OF CORRESPONDING EIGENVALUES, THE MEMBERS OF THE PAIR ARE
  62.      X NOT IDENTICAL.                 //
  63.      X62H UNDER 'INVIT' ARE TWO NUMBERS.  THE FIRST NUMBER, AN INTEGER,,
  64.      X27H IS THE ABSOLUTE SUM OF THE      /
  65.      X62H ERROR FLAGS RETURNED FROM THE PATH.  THE SECOND NUMBER IS THE,
  66.      X29H MEASURE OF PERFORMANCE BASED    /
  67.      X41H UPON THE RESIDUAL COMPUTED FOR THE PATH.//
  68.      X62H -1.0  AS THE MEASURE OF PERFORMANCE IS PRINTED IF AN ERROR IN,
  69.      X27H THE CORRESPONDING PATH HAS      /
  70.      X47H PREVENTED THE COMPUTATION OF THE EIGENVECTORS./)
  71.       WRITE(IWRITE,3)
  72.     3 FORMAT(45H0FOR REDUCTIONS BY ELEMENTARY TRANSFORMATIONS/
  73.      X95H THE  HQR2   PATH WITH    BALANCING USES THE EISPACK CODES  BAL
  74.      XANC-ELMHES-ELTRAN-HQR2  -BALBAK,/
  75.      X38H AS CALLED FROM DRIVER SUBROUTINE  RG. /
  76.      X95H THE   HQR   PATH WITH    BALANCING USES THE EISPACK CODES  BAL
  77.      XANC-ELMHES-HQR   ,              /
  78.      X38H AS CALLED FROM DRIVER SUBROUTINE  RG. /
  79.      X103H THE  INVIT  PATH WITH    BALANCING USES THE EISPACK CODES  BA
  80.      XLANC-ELMHES-HQR   -INVIT -ELMBAK-BALBAK.  /
  81.      X95H THE  HQR2   PATH WITH NO BALANCING USES THE EISPACK CODES  ELM
  82.      XHES-ELTRAN-HQR2  .              /
  83.      X95H THE   HQR   PATH WITH NO BALANCING USES THE EISPACK CODES  ELM
  84.      XHES-HQR   .                     /
  85.      X95H THE  INVIT  PATH WITH NO BALANCING USES THE EISPACK CODES  ELM
  86.      XHES-HQR   -INVIT -ELMBAK.       )
  87.       WRITE(IWRITE,4)
  88.     4 FORMAT(45H0FOR REDUCTIONS BY ORTHOGONAL TRANSFORMATIONS/
  89.      X95H THE  HQR2   PATH WITH    BALANCING USES THE EISPACK CODES  BAL
  90.      XANC-ORTHES-ORTRAN-HQR2  -BALBAK./
  91.      X95H THE   HQR   PATH WITH    BALANCING USES THE EISPACK CODES  BAL
  92.      XANC-ORTHES-HQR   .              /
  93.      X103H THE  INVIT  PATH WITH    BALANCING USES THE EISPACK CODES  BA
  94.      XLANC-ORTHES-HQR   -INVIT -ORTBAK-BALBAK.  /
  95.      X95H THE  HQR2   PATH WITH NO BALANCING USES THE EISPACK CODES  ORT
  96.      XHES-ORTRAN-HQR2  .              /
  97.      X95H THE   HQR   PATH WITH NO BALANCING USES THE EISPACK CODES  ORT
  98.      XHES-HQR   .                     /
  99.      X95H THE  INVIT  PATH WITH NO BALANCING USES THE EISPACK CODES  ORT
  100.      XHES-HQR   -INVIT- ORTBAK.       )
  101.       WRITE(IWRITE,15)
  102.    15 FORMAT(1X,21HD.P. VERSION 04/15/83 )
  103.     5 FORMAT( 53H1       TABULATION OF THE ERROR FLAG  ERROR  AND THE ,
  104.      X    31HMEASURE OF PERFORMANCE  Y  FOR /5X,
  105.      X    56HTHE  EISPACK  CODES.  THIS RUN DISPLAYS THESE STATISTICS ,
  106.      X    30H FOR REAL GENERAL MATRICES.      //
  107.      X    26X,9HBALANCING,15X,12HNO BALANCING,10X,16HVECTORS SELECTED/
  108.      X    18X,3(24H------------------------,1X)/
  109.      X    16H ORDER LOW UPP T,2X,2(25H   HQR2   HQR     INVIT      ))
  110.    10 CALL RMATIN(NM,N,A,AHOLD,0)
  111.       READ(IREADC,50) MM,(SELECT(I),I=1,N)
  112.    50 FORMAT(I4/(72L1))
  113.       DO  60  I = 1,N
  114.    60 SLHOLD(I) = SELECT(I)
  115. C
  116. C     MM  AND  SELECT  ARE READ FROM SYSIN AFTER THE MATRIX IS
  117. C     GENERATED.  MM  SPECIFIES TO  INVIT  THE MAXIMUM NUMBER OF
  118. C     EIGENVECTORS THAT WILL BE COMPUTED.  SELECT  CONTAINS INFORMATION
  119. C     ABOUT WHICH EIGENVECTORS ARE DESIRED.
  120. C
  121.       DO  300  ICALL = 1,12
  122.          ICALL2 = MOD(ICALL,3)
  123.          IF( ICALL2 .NE. 0 ) GO TO 80
  124.          DO 70 I = 1,N
  125.    70    SELECT(I) = SLHOLD(I)
  126.    80    IF( ICALL .NE. 1 )  CALL  RMATIN(NM,N,A,AHOLD,1)
  127.          GO TO  (90,100,110,120,130,140,160,165,170,175,180,185), ICALL
  128. C
  129. C     IF  HQR  PATH IS TAKEN THEN  HQR2  PATH MUST ALSO BE TAKEN
  130. C     IN ORDER THAT COMPARISON OF EIGENVALUES BE POSSIBLE.
  131. C
  132. C
  133. C     RGEWZ
  134. C     INVOKED FROM DRIVER SUBROUTINE  RG.
  135. C
  136.    90    ICT = 1
  137.          CALL  RG(NM,N,A,WR,WI,1,Z,INT,SCALE,ERROR)
  138.          IERR(ICT) = ERROR
  139.          IF( ERROR .NE. 0 ) GO TO 220
  140.          GO TO  190
  141. C
  142. C     RGEW
  143. C     INVOKED FROM DRIVER SUBROUTINE  RG.
  144. C
  145.   100    IF( IERR(1) .NE. 0 ) GO TO 300
  146.          CALL  RG(NM,N,A,WR1,WI1,0,A,INT,SCALE,ERROR)
  147.          IERR(1) = ERROR
  148.          MATCH(1) = DIFF
  149.          DO  101  I = 1,N
  150.             IF( WR(I) .NE. WR1(I) .OR. WI(I) .NE. WI1(I) )  GO TO  300
  151.   101    CONTINUE
  152.          MATCH(1) = SAME
  153.          GO TO  300
  154. C
  155. C     RGEW1Z
  156. C
  157.   110    ICT = 2
  158.          CALL  BALANC(NM,N,A,LOW,UPP,SCALE)
  159.          CALL  ELMHES(NM,N,LOW,UPP,A,INT)
  160.          DO  112  I = 1,N
  161.             DO  112  J = 1,N
  162.   112         ASAVE(I,J) = A(I,J)
  163.          CALL  HQR(NM,N,LOW,UPP,A,WR,WI,ERROR)
  164.          IERR(ICT) = ERROR
  165.          IF( ERROR .NE. 0 )  GO TO  220
  166.          CALL  INVIT(NM,N,ASAVE,WR,WI,SELECT,MM,M,Z,ERROR,RM1,RV1,RV2)
  167.          IERR(ICT) = IABS(ERROR)
  168.          CALL  ELMBAK(NM,LOW,UPP,ASAVE,INT,M,Z)
  169.          CALL  BALBAK(NM,N,LOW,UPP,SCALE,M,Z)
  170.          GO TO  190
  171. C
  172. C     RGNEWZ
  173. C
  174.   120    ICT = 3
  175.          CALL  ELMHES(NM,N,1,N,A,INT)
  176.          CALL  ELTRAN(NM,N,1,N,A,INT,Z)
  177.          CALL  HQR2(NM,N,1,N,A,WR,WI,Z,ERROR)
  178.          IERR(ICT) = ERROR
  179.          IF( ERROR .NE. 0 ) GO TO 220
  180.          DO  121  I = 1,N
  181.             WR1(I) = WR(I)
  182.   121       WI1(I) = WI(I)
  183.          GO TO  190
  184. C
  185. C     RGNEW
  186. C
  187.   130    IF( IERR(3) .NE. 0 ) GO TO 300
  188.          CALL  ELMHES(NM,N,1,N,A,INT)
  189.          CALL  HQR(NM,N,1,N,A,WR,WI,ERROR)
  190.          IERR(3) = ERROR
  191.          MATCH(2) = DIFF
  192.          DO  131  I = 1,N
  193.             IF( WR(I) .NE. WR1(I) .OR. WI(I) .NE. WI1(I) )  GO TO  300
  194.   131    CONTINUE
  195.          MATCH(2) = SAME
  196.          GO TO  300
  197. C
  198. C     RGNEW1Z
  199. C
  200.   140    ICT = 4
  201.          CALL  ELMHES(NM,N,1,N,A,INT)
  202.          DO  142  I = 1,N
  203.             DO  142  J = 1,N
  204.   142         ASAVE(I,J) = A(I,J)
  205.          CALL  HQR(NM,N,1,N,A,WR,WI,ERROR)
  206.          IERR(ICT) = ERROR
  207.          IF( ERROR .NE. 0 )  GO TO  220
  208.          CALL  INVIT(NM,N,ASAVE,WR,WI,SELECT,MM,M,Z,ERROR,RM1,RV1,RV2)
  209.          IERR(ICT) = IABS(ERROR)
  210.          CALL  ELMBAK(NM,1,N,ASAVE,INT,M,Z)
  211.          GO TO 190
  212. C
  213. C     RGOWZ
  214. C
  215.   160    ICT = 5
  216.          CALL  BALANC(NM,N,A,LOW,UPP,SCALE)
  217.          CALL  ORTHES(NM,N,LOW,UPP,A,ORT)
  218.          CALL  ORTRAN(NM,N,LOW,UPP,A,ORT,Z)
  219.          CALL  HQR2(NM,N,LOW,UPP,A,WR,WI,Z,ERROR)
  220.          IERR(ICT) = ERROR
  221.          IF( ERROR .NE. 0 ) GO TO 220
  222.          DO  161  I = 1,N
  223.             WR1(I) = WR(I)
  224.   161       WI1(I) = WI(I)
  225.          CALL  BALBAK(NM,N,LOW,UPP,SCALE,N,Z)
  226.          GO TO  190
  227. C
  228. C     RGOW
  229. C
  230.   165    IF( IERR(5) .NE. 0 ) GO TO 300
  231.          CALL  BALANC(NM,N,A,LOW,UPP,SCALE)
  232.          CALL  ORTHES(NM,N,LOW,UPP,A,ORT)
  233.          CALL  HQR(NM,N,LOW,UPP,A,WR,WI,ERROR)
  234.          IERR(5) = ERROR
  235.          MATCH(3) = DIFF
  236.          DO  166  I = 1,N
  237.             IF( WR(I) .NE. WR1(I) .OR. WI(I) .NE. WI1(I) )  GO TO  300
  238.   166    CONTINUE
  239.          MATCH(3) = SAME
  240.          GO TO  300
  241. C
  242. C     RGOW1Z
  243. C
  244.   170    ICT = 6
  245.          CALL  BALANC(NM,N,A,LOW,UPP,SCALE)
  246.          CALL  ORTHES(NM,N,LOW,UPP,A,ORT)
  247.          DO  172  I = 1,N
  248.             DO  172  J = 1,N
  249.   172         ASAVE(I,J) = A(I,J)
  250.          CALL  HQR(NM,N,LOW,UPP,A,WR,WI,ERROR)
  251.          IERR(ICT) = ERROR
  252.          IF( ERROR .NE. 0 )  GO TO  220
  253.          CALL  INVIT(NM,N,ASAVE,WR,WI,SELECT,MM,M,Z,ERROR,RM1,RV1,RV2)
  254.          IERR(ICT) = IABS(ERROR)
  255.          CALL  ORTBAK(NM,LOW,UPP,ASAVE,ORT,M,Z)
  256.          CALL  BALBAK(NM,N,LOW,UPP,SCALE,M,Z)
  257.          GO TO  190
  258. C
  259. C     RGNOWZ
  260. C
  261.   175    ICT = 7
  262.          CALL  ORTHES(NM,N,1,N,A,ORT)
  263.          CALL  ORTRAN(NM,N,1,N,A,ORT,Z)
  264.          CALL  HQR2(NM,N,1,N,A,WR,WI,Z,ERROR)
  265.          IERR(ICT) = ERROR
  266.          IF( ERROR .NE. 0 ) GO TO 220
  267.          DO  176  I = 1,N
  268.             WR1(I) = WR(I)
  269.   176       WI1(I) = WI(I)
  270.          GO TO  190
  271. C
  272. C     RGNOW
  273. C
  274.   180    IF( IERR(7) .NE. 0 ) GO TO 300
  275.          CALL  ORTHES(NM,N,1,N,A,ORT)
  276.          CALL  HQR(NM,N,1,N,A,WR,WI,ERROR)
  277.          IERR(7) = ERROR
  278.          MATCH(4) = DIFF
  279.          DO  181  I = 1,N
  280.             IF( WR(I) .NE. WR1(I) .OR. WI(I) .NE. WI1(I) )  GO TO  300
  281.   181    CONTINUE
  282.          MATCH(4) = SAME
  283.          GO TO  300
  284. C
  285. C     RGNOW1Z
  286. C
  287.   185    ICT = 8
  288.          CALL  ORTHES(NM,N,1,N,A,ORT)
  289.          DO  187  I = 1,N
  290.             DO  187  J = 1,N
  291.   187         ASAVE(I,J) = A(I,J)
  292.          CALL  HQR(NM,N,1,N,A,WR,WI,ERROR)
  293.          IERR(ICT) = ERROR
  294.          IF( ERROR .NE. 0 )  GO TO  220
  295.          CALL  INVIT(NM,N,ASAVE,WR,WI,SELECT,MM,M,Z,ERROR,RM1,RV1,RV2)
  296.          IERR(ICT) = IABS(ERROR)
  297.          CALL  ORTBAK(NM,1,N,ASAVE,ORT,M,Z)
  298. C
  299.   190    CALL  RMATIN(NM,N,A,AHOLD,1)
  300.          IF( ICALL2 .EQ. 0 ) CALL RGW1ZR(NM,N,A,WR,WI,SELECT,
  301.      X                                   M,Z,RV1,RESDUL)
  302.          IF( ICALL2 .EQ. 1 ) CALL RGWZR(NM,N,A,WR,WI,Z,RV1,RESDUL)
  303.          DFL = 10 * N
  304.          TCRIT(ICT) = RESDUL/EPSLON(DFL)
  305.          GO TO 300
  306.   220    TCRIT(ICT) = -1.0D0
  307.   300 CONTINUE
  308. C
  309.       IF( MOD(LCOUNT,16) .EQ. 0 ) WRITE(IWRITE,5)
  310.       LCOUNT = LCOUNT + 1
  311.       WRITE(IWRITE,520)
  312.      X              N,LOW,UPP,ELEM,(IERR(2*I-1),TCRIT(2*I-1),MATCH(I),
  313.      X              IERR(2*I),TCRIT(2*I),I=1,2),(SELECT(I),I=1,N)
  314.   520 FORMAT(I4,1X,2I4,2X,A1,2(I4,F6.3,1X,A4,I4,F6.3),
  315.      X       4X,21L1/(72X,20L1))
  316.       WRITE(IWRITE,530) ORTH,(IERR(2*I-1),TCRIT(2*I-1),MATCH(I),
  317.      X              IERR(2*I),TCRIT(2*I),I=3,4)
  318.   530 FORMAT(15X,A1,2(I4,F6.3,1X,A4,I4,F6.3))
  319.       GO TO  10
  320.       END
  321.       SUBROUTINE RGWZR(NM,N,A,WR,WI,Z,NORM,RESDUL)
  322. C
  323.       DOUBLE PRECISION NORM(N),WR(N),WI(N),A(NM,N),Z(NM,N),
  324.      X       NORMA,TNORM,XR,XI,S,SUM,SUMA,SUMZ,SUMR,SUMI,RESDUL,PYTHAG
  325.       LOGICAL CPLX
  326. C
  327. C     THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
  328. C     A*Z-Z*DIAG(W) WHERE A IS A REAL GENERAL MATRIX, Z IS
  329. C     A MATRIX WHICH CONTAINS THE EIGENVECTORS OF A, AND W
  330. C     IS A VECTOR WHICH CONTAINS THE EIGENVALUES OF A. ALL
  331. C     NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS.
  332. C
  333. C     THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RGWZR).
  334. C
  335. C     INPUT.
  336. C
  337. C        NM IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
  338. C           AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT;
  339. C
  340. C        N IS THE ORDER OF THE MATRIX  A;
  341. C
  342. C        A(NM,N) IS A REAL GENERAL MATRIX;
  343. C
  344. C        Z(NM,N) IS A MATRIX WHICH CONTAINS THE REAL AND IMAGINARY PARTS
  345. C           OF THE EIGENVECTORS OF A.  THE I-TH COLUMN OF Z IS A REAL
  346. C           EIGENVECTOR IF THE I-TH EIGENVALUE IS REAL.  IF THE I-TH
  347. C           EIGENVALUE IS COMPLEX AND THE FIRST OF A CONJUGATE PAIR, THE
  348. C           I-TH AND (I+1)-TH COLUMNS OF  Z  CONTAIN THE REAL AND
  349. C           IMAGINARY PARTS OF ITS EIGENVECTORS;
  350. C
  351. C        WR(N), WI(N) ARE ARRAYS CONTAINING THE REAL AND IMAGINARY
  352. C           PARTS OF THE EIGENVALUES OF  A.
  353. C
  354. C     OUTPUT.
  355. C
  356. C        Z(NM,N) IS AN ARRAY WHICH CONTAINS THE NORMALIZED
  357. C           APPROXIMATE EIGENVECTORS OF A.  THE EIGENVECTORS ARE
  358. C           NORMALIZED USING THE 1-NORM IN SUCH A WAY THAT THE FIRST
  359. C           ELEMENT WHOSE MAGNITUDE IS LARGER THAN THE NORM OF THE
  360. C           EIGENVECTOR DIVIDED BY N IS REAL AND POSITIVE;
  361. C
  362. C        NORM(N) IS AN ARRAY SUCH THAT FOR EACH K,
  363. C           NORM(K) = !!A*Z(K)-Z(K)*W(K)!!/(!!A!!*!!Z(K)!!)
  364. C           WHERE Z(K) IS THE K-TH EIGENVECTOR;
  365. C
  366. C        RESDUL IS THE REAL NUMBER
  367. C           !!A*Z-Z*DIAG(W)!!/(!!A!!*!!Z!!).
  368. C
  369. C     ----------------------------------------------------------------
  370. C
  371.       NORMA = 0.0D0
  372.       RESDUL = 0.0D0
  373.       CPLX = .FALSE.
  374. C
  375.       DO 20 I=1,N
  376.          SUMA = 0.0D0
  377. C
  378.          DO 10 L=1,N
  379.    10       SUMA = SUMA + DABS(A(L,I))
  380. C
  381.    20    NORMA = DMAX1(NORMA,SUMA)
  382. C
  383.       IF(NORMA .EQ. 0.0D0) NORMA = 1.0D0
  384. C
  385.       DO 160 I=1,N
  386.          IF(CPLX) GO TO 80
  387.          S = 0.0D0
  388.          SUMZ = 0.0D0
  389.          IF(WI(I) .NE. 0.0D0) GO TO 90
  390. C        ..........THIS IS FOR REAL EIGENVALUES..........
  391.          DO 40 L=1,N
  392.             SUMZ = SUMZ + DABS(Z(L,I))
  393.             SUM = -WR(I)*Z(L,I)
  394. C
  395.             DO 30 K=1,N
  396.    30          SUM = SUM + A(L,K)*Z(K,I)
  397. C
  398.    40       S = S + DABS(SUM)
  399. C
  400.          NORM(I) = SUMZ
  401. C        ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE
  402. C                  WILL ALWAYS EXIST AN ELEMENT IN THE VECTOR Z(I)
  403. C                  LARGER THAN !!Z(I)!!/N..........
  404.          DO 50 L=1,N
  405.             IF(DABS(Z(L,I)) .GE. NORM(I)/N) GO TO 60
  406.    50       CONTINUE
  407. C
  408.    60    TNORM = DSIGN(NORM(I),Z(L,I))
  409. C
  410.          DO 70 L=1,N
  411.    70       Z(L,I) = Z(L,I)/TNORM
  412. C
  413.          NORM(I) = S/(NORM(I)*NORMA)
  414.    80    CPLX = .FALSE.
  415.          GO TO 150
  416. C        ..........THIS IS FOR COMPLEX EIGENVALUES..........
  417.    90    DO 110 L = 1,N
  418.             SUMZ = SUMZ + PYTHAG(Z(L,I),Z(L,I+1))
  419.             SUMR = -WR(I)*Z(L,I) + WI(I)*Z(L,I+1)
  420.             SUMI = -WR(I)*Z(L,I+1) - WI(I)*Z(L,I)
  421. C
  422.             DO 100 K=1,N
  423.                SUMR = SUMR + A(L,K)*Z(K,I)
  424.   100          SUMI = SUMI + A(L,K)*Z(K,I+1)
  425. C
  426.   110       S = S + PYTHAG(SUMR,SUMI)
  427. C
  428.          NORM(I) = SUMZ
  429. C        ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE WILL
  430. C                  ALWAYS EXIST AN ELEMENT IN THE VECTOR Z(I) LARGER
  431. C                  THAN !!Z(I)!!/N..........
  432.          DO 120 L=1,N
  433.             IF(PYTHAG(Z(L,I),Z(L,I+1)) .GE. NORM(I)/N)
  434.      1      GO TO 130
  435.   120       CONTINUE
  436.   130    CONTINUE
  437. C
  438.          XR = PYTHAG(Z(L,I),Z(L,I+1))            
  439.          XR = NORM(I)*Z(L,I)/PYTHAG(Z(L,I),Z(L,I+1))
  440.          XI = NORM(I)*Z(L,I+1)/PYTHAG(Z(L,I),Z(L,I+1))
  441. C
  442.          DO 140 L= 1,N
  443.             CALL CDIV(Z(L,I),Z(L,I+1),XR,XI,Z(L,I),Z(L,I+1))
  444.   140    CONTINUE
  445. C
  446.          NORM(I) = S/(NORM(I)*NORMA)
  447.          NORM(I+1) = NORM(I)
  448.          CPLX = .TRUE.
  449.   150    RESDUL = DMAX1(NORM(I),RESDUL)
  450.   160    CONTINUE
  451. C
  452.       RETURN
  453.       END
  454.       SUBROUTINE RGW1ZR(NM,N,A,WR,WI,SELECT,M,Z,NORM,RESDUL)
  455. C
  456.       DOUBLE PRECISION NORM(N),WR(N),WI(N),A(NM,N),Z(NM,M),
  457.      X       NORMA,TNORM,XR,XI,S,SUM,SUMA,SUMZ,SUMR,SUMI,RESDUL,PYTHAG
  458.       LOGICAL SELECT(N), CPLX
  459. C
  460. C     THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
  461. C     A*Z-Z*DIAG(W)   WHERE  A  IS A REAL GENERAL MATRIX,  Z  IS
  462. C     A MATRIX WHICH CONTAINS SOME OF THE EIGENVECTORS OF  A,
  463. C     AND  W  IS A VECTOR WHICH CONTAINS THE EIGENVALUES OF  A.
  464. C     ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS.
  465. C
  466. C     THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RGW1ZR).
  467. C
  468. C     INPUT.
  469. C
  470. C        NM  IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
  471. C           AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT.
  472. C
  473. C        N  IS THE ORDER OF THE MATRIX  A.
  474. C
  475. C        A(NM,N)  IS A REAL GENERAL MATRIX.
  476. C
  477. C        Z(NM,M)  IS A MATRIX WHICH CONTAINS THE REAL AND IMAGINARY
  478. C           OF THE EIGENVECTORS OF  A.  THE EIGENVECTORS ARE STORED
  479. C           PARTS IN PACKED FORM, THAT IS, USING ONE COLUMN FOR EACH
  480. C           REAL EIGENVECTOR AND USING TWO COLUMNS FOR EACH COMPLEX
  481. C           EIGENVECTOR, THE EIGENVECTORS CORRESPONDING TO THE FLAGGED
  482. C           EIGENVALUES ARE STORED IN CONSECUTIVE COLUMNS, STARTING IN
  483. C           COLUMN  1.  (NOTE. FOR THE DEFINITION OF FLAGGED EIGENVALUES
  484. C           SEE THE DESCRIPTION OF THE SELECT ARRAY.)
  485. C
  486. C        WR(N),WI(N)  ARE ARRAYS CONTAINING THE REAL AND IMAGINARY
  487. C           PARTS OF THE EIGENVALUES OF  A.
  488. C
  489. C        SELECT(N)  IS A LOGICAL ARRAY WHICH SPECIFIES THE
  490. C           FLAGGED EIGENVALUES WHOSE EIGENVECTORS ARE
  491. C           CONTAINED IN THE MATRIX  Z.  THE I-TH EIGENVALUE
  492. C           (WR(I),WI(I))  IS FLAGGED IF  SELECT(I) = .TRUE.
  493. C
  494. C        M  IS THE NUMBER OF COLUMNS OF  Z  CONTAINING EIGENVECTORS
  495. C           OF  A.
  496. C
  497. C     OUTPUT.
  498. C
  499. C        Z(NM,M)  IS AN ARRAY WHICH CONTAINS THE NORMALIZED
  500. C           APPROXIMATE EIGENVECTORS OF  A  IN PACKED FORM AS
  501. C           DESCRIBED ABOVE.  THE EIGENVECTORS ARE NORMALIZED
  502. C           USING THE 1-NORM IN SUCH A WAY THAT THE FIRST ELEMENT
  503. C           WHOSE MAGNITUDE IS LARGER THAN THE NORM OF THE
  504. C           EIGENVECTOR DIVIDED BY  N  IS REAL AND POSITIVE.
  505. C
  506. C        NORM(N)  IS AN ARRAY SUCH THAT WHENEVER SELECT(K)  = .TRUE.
  507. C           NORM(K) = !!A*Z(K)-Z(K)*W(K)!!/(!!A!!*!!Z(K)!!)
  508. C           WHERE  Z(K)  IS THE EIGENVECTOR CORRESPONDING TO THE
  509. C           K-TH EIGENVALUE  (WR(K),WI(K)).
  510. C
  511. C        RESDUL IS THE REAL NUMBER
  512. C           !!A*Z-Z*DIAG(W)!!/(!!A!!*!!Z!!)
  513. C
  514. C     ----------------------------------------------------------------
  515. C
  516.       KCOM = 0
  517.       NORMA = 0.0D0
  518.       RESDUL = 0.0D0
  519.       IF( M .EQ. 0 ) RETURN
  520.       CPLX = .FALSE.
  521.       II=1
  522. C
  523.       DO 20 I=1,N
  524.          SUMA = 0.0D0
  525. C
  526.          DO 10 L=1,N
  527.    10       SUMA = SUMA + DABS(A(L,I))
  528. C
  529.    20    NORMA = DMAX1(NORMA,SUMA)
  530. C
  531.       IF(NORMA .EQ. 0.0D0) NORMA = 1.0D0
  532. C
  533.       DO 160 I=1,N
  534.          IF(CPLX) GO TO 80
  535.          S = 0.0D0
  536.          SUMZ = 0.0D0
  537.          IF(WI(I) .NE. 0.0D0) GO TO 90
  538.          IF( .NOT. SELECT(I)) GO TO 160
  539. C        ..........THIS IS FOR REAL EIGENVALUES..........
  540.          DO 40 L=1,N
  541.             SUMZ = SUMZ + DABS(Z(L,II))
  542.             SUM = -WR(I)*Z(L,II)
  543. C
  544.             DO 30 K=1,N
  545.    30          SUM = SUM + A(L,K)*Z(K,II)
  546. C
  547.    40       S = S + DABS(SUM)
  548. C
  549.          IF(SUMZ .NE. 0.0D0) GO TO 45
  550.          NORM(I) = -1.0D0
  551.          GO TO 150
  552.    45    NORM(I) = SUMZ
  553. C        ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE
  554. C                  WILL ALWAYS EXIST AN ELEMENT IN THE VECTOR Z(II)
  555. C                  LARGER THAN !!Z(II)!!/N..........
  556.          DO 50 L=1,N
  557.             IF(DABS(Z(L,II)) .GE. NORM(I)/N) GO TO 60
  558.    50       CONTINUE
  559. C
  560.    60    TNORM = DSIGN(NORM(I),Z(L,II))
  561. C
  562.          DO 70 J=1,N
  563.    70       Z(J,II) = Z(J,II)/TNORM
  564. C
  565.          NORM(I) = S/(NORM(I)*NORMA)
  566.          GO TO 150
  567.    80    CPLX = .FALSE.
  568.          GO TO 155
  569. C        ..........THIS IS FOR COMPLEX EIGENVALUES..........
  570.    90    IF(SELECT(I)) GO TO 95
  571.          KCOM = MOD(KCOM+1,2)
  572.          GO TO 160
  573. C
  574.    95    DO 110 L = 1,N
  575.             SUMZ = SUMZ + PYTHAG(Z(L,II),Z(L,II+1))
  576.             SUMR = -WR(I)*Z(L,II) + WI(I)*Z(L,II+1)
  577.             SUMI = -WR(I)*Z(L,II+1) - WI(I)*Z(L,II)
  578. C
  579.             DO 100 K=1,N
  580.                SUMR = SUMR + A(L,K)*Z(K,II)
  581.   100          SUMI = SUMI + A(L,K)*Z(K,II+1)
  582. C
  583.   110       S = S + PYTHAG(SUMR,SUMI)
  584. C
  585.          IF(SUMZ .NE. 0.0D0) GO TO 115
  586.          NORM(I) = -1.0D0
  587.          GO TO 145
  588.   115    NORM(I) = SUMZ
  589. C        ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE WILL
  590. C                  ALWAYS EXIST AN ELEMENT IN THE VECTOR Z(II) LARGER
  591. C                  THAN !!Z(II)!!/N..........
  592.          DO 120 L=1,N
  593.             IF(PYTHAG(Z(L,II),Z(L,II+1)) .GE. NORM(I)/N)
  594.      1      GO TO 130
  595.   120       CONTINUE
  596. C
  597.   130    XR = NORM(I)*Z(L,II)/PYTHAG(Z(L,II),Z(L,II+1))
  598.          XI = NORM(I)*Z(L,II+1)/PYTHAG(Z(L,II),Z(L,II+1))
  599. C
  600.          DO 140 J= 1,N
  601.             CALL CDIV(Z(J,II),Z(J,II+1),XR,XI,Z(J,II),Z(J,II+1))
  602.   140    CONTINUE
  603. C
  604.          NORM(I) = S/(NORM(I)*NORMA)
  605.   145    CPLX = .TRUE.
  606.          IF(KCOM .EQ. 0) GO TO 150
  607.          CPLX = .FALSE.
  608.          II = II+1
  609.   150    RESDUL = DMAX1(NORM(I),RESDUL)
  610.   155    II = II+1
  611.   160    CONTINUE
  612. C
  613.       RETURN
  614.       END
  615.       SUBROUTINE RMATIN(NM,N,A,AHOLD,INITIL)
  616. C
  617. C     THIS INPUT SUBROUTINE READS A REAL MATRIX FROM SYSIN OF
  618. C     ORDER N.
  619. C     TO GENERATE THE MATRIX  A  INITIALLY,  INITIL  IS TO BE 0.
  620. C     TO REGENERATE THE MATRIX  A  FOR THE PURPOSE OF THE RESIDUAL
  621. C     CALCULATION,  INITIL  IS TO BE  1.
  622. C
  623. C     THIS ROUTINE IS CATALOGUED AS  EISPDRV4(RGREADI).
  624. C
  625.       DOUBLE PRECISION A(NM,NM),AHOLD(NM,NM)
  626.       INTEGER  IA( 20)
  627.       DATA IREADA/1/,IWRITE/6/
  628. C
  629.       IF( INITIL .EQ. 1 )  GO TO  30
  630.       READ(IREADA,5)N, M
  631.     5 FORMAT(I6,6X,I6)
  632.       IF( N .EQ. 0 )  GO TO  70
  633.       IF (M .NE. 1) GO TO 14
  634.       DO  13  I = 1,N
  635.          READ(IREADA,16) (IA(J), J=I,N)
  636.          DO  12  J = I,N
  637.            A(I,J) = IA(J)
  638.    12      A(J,I) = A(I,J)
  639.    13 CONTINUE
  640.       GO TO 19
  641.    14 DO  17  I = 1,N
  642.          READ(IREADA,16) (IA(J), J=1,N)
  643.    16    FORMAT(6I12)
  644.          DO  17  J = 1,N
  645.    17      A(I,J) = IA(J)
  646.    19 DO  20  I = 1,N
  647.          DO  20  J = 1,N
  648.    20      AHOLD(I,J) = A(I,J)
  649.       RETURN
  650.    30 DO  40  I = 1,N
  651.          DO  40  J = 1,N
  652.    40      A(I,J) = AHOLD(I,J)
  653.       RETURN
  654.    70 WRITE(IWRITE,80)
  655.    80 FORMAT(46H0END OF DATA FOR SUBROUTINE  RMATIN(RGREADI). /1H1)
  656.       STOP
  657.       END
  658.