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 / rggtest.f < prev    next >
Text File  |  1996-09-28  |  11KB  |  320 lines
  1. C
  2. C     THIS DRIVER TESTS  EISPACK  FOR THE REAL GENERAL
  3. C     GENERALIZED EIGENPROBLEM  A*X = (LAMBDA)*B*X, SUMMARIZING THE
  4. C     FIGURES OF MERIT FOR ALL PATHS.
  5. C
  6. C     THIS DRIVER IS CATALOGUED AS  EISPDRV4(RGGSUMAR).
  7. C
  8. C     THE DIMENSION OF  A,  B, AND  Z  SHOULD BE  NM  BY  NM.
  9. C     THE DIMENSION OF  ALFR,ALFI,BETA, AND  NORM  SHOULD BE  NM.
  10. C     THE DIMENSION OF  ALFR1,ALFI1, AND  BETA1  SHOULD BE  NM.
  11. C     THE DIMENSION OF  AHOLD  AND  BHOLD  SHOULD BE  NM  BY  NM.
  12. C     HERE NM = 20.
  13. C
  14.       DOUBLE PRECISION A( 20, 20),B( 20, 20),Z( 20, 20),
  15.      X        ALFR( 20),ALFI( 20),BETA( 20),NORM( 20),TCRIT( 1),
  16.      X        EPSLON,RESDUL,DFL
  17.       DOUBLE PRECISION AHOLD( 20, 20),BHOLD( 20, 20)
  18.       DOUBLE PRECISION ALFR1( 20),ALFI1( 20),BETA1( 20)
  19.       INTEGER IERR( 1),ERROR
  20.       INTEGER MATCH( 1),SAME,DIFF
  21.       DATA IREAD1/1/,IREAD2/2/,IWRITE/6/
  22.       DATA SAME,DIFF/4HSAME,4HDIFF/
  23. C
  24.       OPEN(UNIT=IREAD1,FILE='FILE45')
  25.       OPEN(UNIT=IREAD2,FILE='FILE46')
  26.       REWIND IREAD1
  27.       REWIND IREAD2
  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-)/24H ORDER  QZIT(T)  QZIT(F) /
  34.      X1H ,95(1H-)//
  35.      X55H UNDER 'ORDER' IS THE ORDER OF EACH TEST MATRIX SYSTEM. //
  36.      X95H UNDER 'QZIT(T)  QZIT(F)' ARE TWO NUMBERS AND A KEYWORD.  THE F
  37.      XIRST NUMBER, AN INTEGER, IS THE /
  38.      X95H ERROR FLAG RETURNED FROM  QZIT.  THE SECOND NUMBER IS THE MEAS
  39.      XURE OF PERFORMANCE BASED UPON   /
  40.      X95H THE RESIDUAL COMPUTED FOR THE  QZIT(T) (BOTH EIGENVALUES AND E
  41.      XIGENVECTORS) PATH.  THE KEYWORD /
  42.      X95H REPORTS THE DUPLICATION OF COMPUTED EIGENVALUES FROM THE  QZIT
  43.      X(T)  AND  QZIT(F) (EIGENVALUES  /
  44.      X95H ONLY) PATHS.  'SAME' MEANS THAT THE EIGENVALUES ARE EXACT DUPL
  45.      XICATES.  'DIFF' MEANS THAT FOR  /
  46.      X95H AT LEAST ONE PAIR OF CORRESPONDING EIGENVALUES, THE MEMBERS OF
  47.      X THE PAIR ARE NOT IDENTICAL.    )
  48.       WRITE(IWRITE,2)
  49.     2 FORMAT(81H0-1.0  AS THE MEASURE OF PERFORMANCE IS PRINTED IF AN ER
  50.      XROR IN THE PATH HAS        /
  51.      X58H PREVENTED THE COMPUTATION OF THE EIGENVALUES AND VECTORS. //
  52.      X60H THE  QZIT(T)  PATH USES THE EISPACK CODES  QZHES-QZIT-QZVAL
  53.      X,8H-QZVEC, /
  54.      X39H AS CALLED FROM DRIVER SUBROUTINE  RGG. /
  55.      X61H THE  QZIT(F)  PATH USES THE EISPACK CODES  QZHES-QZIT-QZVAL, /
  56.      X39H AS CALLED FROM DRIVER SUBROUTINE  RGG. )
  57.       WRITE(IWRITE,15)
  58.    15 FORMAT(1X,21HD.P. VERSION 04/15/83 )
  59.     5 FORMAT( 53H1       TABULATION OF THE ERROR FLAG  ERROR  AND THE ,
  60.      X    31HMEASURE OF PERFORMANCE  Y  FOR /5X,
  61.      X    56HTHE  EISPACK  CODES.  THIS RUN DISPLAYS THESE STATISTICS ,
  62.      X    21H FOR THE REAL GENERAL/4X,30H GENERALIZED EIGENPROBLEM  A*X,
  63.      X    15H= (LAMBDA)*B*X./
  64.      X    24H0ORDER  QZIT(T)  QZIT(F) )
  65.    10 CALL  RMATIN(NM,N,A,B,AHOLD,BHOLD,0)
  66.       DO  510  ICALL = 1,2
  67.          IF( ICALL .NE. 1 )  CALL  RMATIN(NM,N,A,B,AHOLD,BHOLD,1)
  68.          GO TO  (70,100), ICALL
  69. C
  70. C     RGGWZ
  71. C     INVOKED FROM DRIVER SUBROUTINE  RGG.
  72. C
  73.    70    ICT = 1
  74.          CALL  RGG(NM,N,A,B,ALFR,ALFI,BETA,1,Z,ERROR)
  75.          IERR(ICT) = ERROR
  76.          IF( ERROR .NE. 0 ) GO TO 500
  77.          GO TO 190
  78. C
  79. C     RGGW
  80. C     INVOKED FROM DRIVER SUBROUTINE  RGG.
  81. C
  82.   100    IF( IERR(1) .NE. 0 ) GO TO 510
  83.          CALL  RGG(NM,N,A,B,ALFR1,ALFI1,BETA1,0,A,ERROR)
  84.          IERR(1) = ERROR
  85.          MATCH(1) = DIFF
  86.          DO  101  I = 1,N
  87.             IF( ALFR(I) .NE. ALFR1(I) .OR. ALFI(I) .NE. ALFI1(I) .OR.
  88.      X          BETA(I) .NE. BETA1(I) ) GO TO 510
  89.   101    CONTINUE
  90.          MATCH(1) = SAME
  91.          GO TO 510
  92. C
  93.   190    CALL  RMATIN(NM,N,A,B,AHOLD,BHOLD,1)
  94.          CALL  RGGWZR(NM,N,A,B,ALFR,ALFI,BETA,Z,NORM,RESDUL)
  95.          DFL = 10 * N
  96.          TCRIT(ICT) = RESDUL/EPSLON(DFL)
  97.          GO TO 510
  98.   500    TCRIT(ICT) = -1.0D0
  99.   510 CONTINUE
  100. C
  101.       IF( MOD(LCOUNT,35) .EQ. 0 ) WRITE(IWRITE,5)
  102.       LCOUNT = LCOUNT + 1
  103.       WRITE(IWRITE,520) N,IERR(1),TCRIT(1),MATCH(1)
  104.   520 FORMAT(I4,I5,F6.3,4X,A4)
  105.       GO TO 10
  106.       END
  107.       SUBROUTINE RGGWZR(NM,N,A,B,ALFR,ALFI,BETA,Z,NORM,RESDUL)
  108. C
  109.       DOUBLE PRECISION NORM(N),ALFR(N),ALFI(N),BETA(N),A(NM,N),
  110.      X       B(NM,N), Z(NM,N), TNORM, XR, XI, S, SUM,  SUMZ, SUMR, SUMI,
  111.      X       SUMQ,RESDUL,NORMAB,SUM1,SUMA,SUMB,NORMA,NORMB,
  112.      X       SUMR2,SUMI2,SUMR3,SUMI3,PYTHAG
  113.       LOGICAL CPLX
  114. C
  115. C     THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
  116. C     A*Z-B*Z*DIAG(W) WHERE A AND B ARE REAL GENERAL MATRICES, Z IS
  117. C     A MATRIX WHICH CONTAINS THE EIGENVECTORS OF THE EIGENPROBLEM
  118. C     A*Z - B*Z*DIAG(W), AND  W  STANDS FOR A VECTOR OF CORRESPONDING
  119. C     EIGENVALUES OF THE EIGENPROBLEM OBTAINED FROM THE VECTORS  ALFR,
  120. C     ALFI,  AND  BETA  BY THE CORRESPONDENCES
  121. C                    W(J) = (ALFR(J) + I*ALFI(J)) / BETA(J) .
  122. C     ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS.
  123. C
  124. C     THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RGGWZR).
  125. C
  126. C     INPUT.
  127. C
  128. C        NM IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
  129. C           AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT;
  130. C
  131. C        N IS THE ORDER OF THE MATRICES  A  AND  B;
  132. C
  133. C        A(NM,N)  AND  B(NM,N)  ARE REAL GENERAL MATRICES;
  134. C
  135. C        Z(NM,N) IS A MATRIX WHICH CONTAINS THE REAL AND IMAGINARY PARTS
  136. C           OF THE EIGENVECTORS OF THE SYSTEM.  THE I-TH COLUMN OF Z IS
  137. C           A REAL EIGENVECTOR IF THE I-TH EIGENVALUE IS REAL.  IF THE
  138. C           I-TH EIGENVALUE IS COMPLEX AND THE FIRST OF A CONJUGATE
  139. C           PAIR, THE I-TH AND (I+1)-TH COLUMNS OF  Z  CONTAIN THE REAL
  140. C           AND IMAGINARY PARTS OF ITS EIGENVECTORS;
  141. C
  142. C        ALFR(N), ALFI(N)  AND BETA(N)  ARE ARRAYS CONTAINING THE
  143. C           COMPONENTS OF THE EIGENVALUES OF THE SYSTEM.
  144. C
  145. C     OUTPUT.
  146. C
  147. C        Z(NM,N) IS AN ARRAY WHICH CONTAINS THE NORMALIZED
  148. C           APPROXIMATE EIGENVECTORS OF THE SYSTEM.  THE EIGENVECTORS
  149. C           ARE NORMALIZED USING THE 1-NORM IN SUCH A WAY THAT THE FIRST
  150. C           ELEMENT WHOSE MAGNITUDE IS LARGER THAN THE NORM OF THE
  151. C           EIGENVECTOR DIVIDED BY N IS REAL AND POSITIVE;
  152. C
  153. C        NORM(N) IS AN ARRAY SUCH THAT FOR EACH K,
  154. C
  155. C                          !!BETA(K)*A*Z(K)-ALFA*B*Z(K)!!
  156. C         NORM(K) = -----------------------------------------------
  157. C                     !!Z(K)!!*(BETA(K)*!!A!!+!ALFA(K)!*!!B!!)
  158. C
  159. C         WHERE Z(K) IS THE K-TH EIGENVECTOR AND ALFA = (ALFR + I*ALFI);
  160. C
  161. C        RESDUL IS THE REAL NUMBER
  162. C           !!BETA*A*Z-ALFA*B*Z!!/(!!Z!!*(BETA*!!A!!+!ALFA!*!!B!!))
  163. C
  164. C     ----------------------------------------------------------------
  165. C
  166.       NORMB = 0.0D0
  167.       NORMA = 0.0D0
  168.       RESDUL = 0.0D0
  169.       CPLX = .FALSE.
  170. C
  171.       DO 20 I = 1,N
  172.         SUMA = 0.0D0
  173.         SUMB = 0.0D0
  174.         DO 10 L = 1,N
  175.           SUMA = SUMA + DABS(A(L,I))
  176.    10     SUMB = SUMB + DABS(B(L,I))
  177.         NORMA = DMAX1(NORMA,SUMA)
  178.    20   NORMB = DMAX1(NORMB,SUMB)
  179. C
  180.       DO 160 I=1,N
  181.          IF(CPLX) GO TO 80
  182.          S = 0.0D0
  183.          SUMZ = 0.0D0
  184.          IF(ALFI(I) .NE. 0.0D0) GO TO 90
  185. C        ..........THIS IS FOR REAL EIGENVALUES..........
  186.          DO 40 L=1,N
  187.             SUMQ = 0.0D0
  188.             SUM1 = 0.0D0
  189.             SUMZ = SUMZ + DABS(Z(L,I))
  190. C
  191.             DO 35 K=1,N
  192.                SUMQ = SUMQ + B(L,K)*Z(K,I)
  193.    35          SUM1 = SUM1 + A(L,K)*Z(K,I)
  194. C
  195.             SUM = -ALFR(I)*SUMQ
  196.    40       S = S + DABS(SUM + SUM1*BETA(I))
  197. C
  198.          NORMAB = NORMA*BETA(I) + NORMB*DABS(ALFR(I))
  199.          IF( NORMAB .EQ. 0.0D0 ) NORMAB = 1.0D0
  200.          NORM(I) = SUMZ
  201.          IF( SUMZ .EQ. 0.0D0 ) GO TO 150
  202. C        ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE
  203. C                  WILL ALWAYS EXIST AN ELEMENT IN THE VECTOR Z(I)
  204. C                  LARGER THAN !!Z(I)!!/N..........
  205.          DO 50 L=1,N
  206.             IF(DABS(Z(L,I)) .GE. NORM(I)/N) GO TO 60
  207.    50       CONTINUE
  208. C
  209.    60    TNORM = DSIGN(NORM(I),Z(L,I))
  210. C
  211.          DO 70 L=1,N
  212.    70       Z(L,I) = Z(L,I)/TNORM
  213. C
  214.          NORM(I) = S/(NORM(I)*NORMAB)
  215.    80    CPLX = .FALSE.
  216.          GO TO 150
  217. C        ..........THIS IS FOR COMPLEX EIGENVALUES..........
  218.    90    DO 110 L = 1,N
  219.             SUMR = 0.0D0
  220.             SUMI = 0.0D0
  221.             SUMR2 = 0.0D0
  222.             SUMI2 = 0.0D0
  223.             SUMZ = SUMZ + PYTHAG(Z(L,I),Z(L,I+1))
  224. C
  225.             DO 100 K=1,N
  226.                SUMR = SUMR + B(L,K)*Z(K,I)
  227.                SUMI = SUMI + B(L,K)*Z(K,I+1)
  228.                SUMR2 = SUMR2 + A(L,K)*Z(K,I)
  229.   100          SUMI2 = SUMI2 + A(L,K)*Z(K,I+1)
  230.             SUMR3 = -ALFR(I)*SUMR + ALFI(I)*SUMI
  231.             SUMI3 = -ALFI(I)*SUMR - ALFR(I)*SUMI
  232. C
  233.   110       S = S+PYTHAG(SUMR3+SUMR2*BETA(I),SUMI3+SUMI2*BETA(I))
  234.             NORMAB = NORMA*BETA(I)+NORMB*PYTHAG(ALFR(I),ALFI(I))
  235.             IF( NORMAB .EQ. 0.0D0 ) NORMAB = 1.0D0
  236. C
  237.          NORM(I) = SUMZ
  238.          IF( SUMZ .EQ. 0.0D0 ) GO TO 150
  239. C        ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE WILL
  240. C                  ALWAYS EXIST AN ELEMENT IN THE VECTOR Z(I) LARGER
  241. C                  THAN !!Z(I)!!/N..........
  242.          DO 120 L=1,N
  243.             IF(PYTHAG(Z(L,I),Z(L,I+1)) .GE. NORM(I)/N)
  244.      1      GO TO 130
  245.   120       CONTINUE
  246. C
  247.   130    XR = NORM(I)*Z(L,I)/PYTHAG(Z(L,I),Z(L,I+1))
  248.          XI = NORM(I)*Z(L,I+1)/PYTHAG(Z(L,I),Z(L,I+1))
  249. C
  250.          DO 140 L= 1,N
  251.             CALL CDIV(Z(L,I),Z(L,I+1),XR,XI,Z(L,I),Z(L,I+1))
  252.   140    CONTINUE
  253. C
  254.          NORM(I) = S/(NORM(I)*NORMAB)
  255.          NORM(I+1) = NORM(I)
  256.          CPLX = .TRUE.
  257.   150    RESDUL = DMAX1(NORM(I),RESDUL)
  258.   160    CONTINUE
  259. C
  260.       RETURN
  261.       END
  262.       SUBROUTINE RMATIN(NM,N,A,B,AHOLD,BHOLD,INITIL)
  263. C
  264. C     THIS INPUT SUBROUTINE READS TWO REAL MATRICES  A  AND  B  FROM
  265. C     SYSIN OF ORDER N.
  266. C     TO GENERATE THE MATRICES INITIALLY,  INITIL  IS TO BE 0.
  267. C     TO REGENERATE THE MATRICES FOR THE PURPOSE OF THE RESIDUAL
  268. C     CALCULATION,  INITIL  IS TO BE  1.
  269. C
  270. C     THIS ROUTINE IS CATALOGUED AS  EISPDRV4(RSGREADI).
  271. C
  272.       DOUBLE PRECISION A(NM,NM),B(NM,NM),AHOLD(NM,NM),BHOLD(NM,NM) 
  273.       INTEGER  IA( 20), IB( 20)
  274.       DATA IREADA/1/,IREADB/2/,IWRITE/6/
  275. C
  276.       IF( INITIL .EQ. 1 )  GO TO  30
  277.       READ(IREADA,5) N, M
  278.     5 FORMAT(I6,6X,I6)
  279.       IF( N .EQ. 0 )  GO TO  70
  280.       IF (M .NE. 1) GO TO 16
  281.       DO  10  I = 1,N
  282.          READ(IREADA,17) (IA(J), J=I,N)  
  283.          DO   9  J = I,N
  284.            A(I,J) = IA(J)
  285.     9      A(J,I) =  A(I,J)
  286.    10 CONTINUE
  287.    11 READ(IREADB,5) N,M
  288.       IF( M .NE. 1 ) GO TO 20
  289.       DO 15 I = 1,N
  290.          READ(IREADB,17) (IB(J), J=I,N)   
  291.          DO 14 J = I,N
  292.             B(I,J)=IB(J)
  293.    14       B(J,I)=B(I,J)
  294.    15 CONTINUE
  295.       GO TO 22
  296.    16 DO  18  I = 1,N
  297.          READ(IREADA,17) (IA(J), J=1,N)  
  298.    17    FORMAT(6I12)
  299.          DO  18  J = 1,N
  300.    18      A(I,J) = IA(J)
  301.       GO TO 11
  302.    20 DO 25 I = 1,N
  303.          READ(IREADB,17) (IB(J),J=1,N) 
  304.          DO 25 J = 1,N
  305.    25    B(I,J) = IB(J)
  306.    22 DO  23  I = 1,N
  307.          DO  23  J = 1,N
  308.            BHOLD(I,J) = B(I,J)
  309.    23      AHOLD(I,J) = A(I,J)
  310.       RETURN
  311.    30 DO  40  I = 1,N
  312.          DO  40  J = 1,N
  313.            B(I,J) = BHOLD(I,J)
  314.    40      A(I,J) = AHOLD(I,J)
  315.       RETURN
  316.    70 WRITE(IWRITE,80)
  317.    80 FORMAT(47H0END OF DATA FOR SUBROUTINE  RMATIN(RSGREADI). /1H1)
  318.       STOP
  319.       END
  320.