C C THIS DRIVER TESTS EISPACK FOR THE REAL GENERAL C GENERALIZED EIGENPROBLEM A*X = (LAMBDA)*B*X, SUMMARIZING THE C FIGURES OF MERIT FOR ALL PATHS. C C THIS DRIVER IS CATALOGUED AS EISPDRV4(RGGSUMAR). C C THE DIMENSION OF A, B, AND Z SHOULD BE NM BY NM. C THE DIMENSION OF ALFR,ALFI,BETA, AND NORM SHOULD BE NM. C THE DIMENSION OF ALFR1,ALFI1, AND BETA1 SHOULD BE NM. C THE DIMENSION OF AHOLD AND BHOLD SHOULD BE NM BY NM. C HERE NM = 20. C REAL A( 20, 20),B( 20, 20),Z( 20, 20), X ALFR( 20),ALFI( 20),BETA( 20),NORM( 20),TCRIT( 1), X EPSLON,RESDUL,DFL REAL AHOLD( 20, 20),BHOLD( 20, 20) REAL ALFR1( 20),ALFI1( 20),BETA1( 20) INTEGER IERR( 1),ERROR INTEGER MATCH( 1) CHARACTER*4 SAME,DIFF DATA IREAD1/1/,IREAD2/2/,IWRITE/6/ DATA SAME,DIFF/'SAME','DIFF'/ C OPEN(IREAD1,FILE='FILE45') OPEN(IREAD2,FILE='FILE46') REWIND IREAD1 REWIND IREAD2 C NM = 20 LCOUNT = 0 WRITE(*,1) 1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S XTATISTICS//1H ,95(1H-)/24H ORDER QZIT(T) QZIT(F) / X1H ,95(1H-)// X55H UNDER 'ORDER' IS THE ORDER OF EACH TEST MATRIX SYSTEM. // X95H UNDER 'QZIT(T) QZIT(F)' ARE TWO NUMBERS AND A KEYWORD. THE F XIRST NUMBER, AN INTEGER, IS THE / X95H ERROR FLAG RETURNED FROM QZIT. THE SECOND NUMBER IS THE MEAS XURE OF PERFORMANCE BASED UPON / X95H THE RESIDUAL COMPUTED FOR THE QZIT(T) (BOTH EIGENVALUES AND E XIGENVECTORS) PATH. THE KEYWORD / X95H REPORTS THE DUPLICATION OF COMPUTED EIGENVALUES FROM THE QZIT X(T) AND QZIT(F) (EIGENVALUES / X95H ONLY) PATHS. 'SAME' MEANS THAT THE EIGENVALUES ARE EXACT DUPL XICATES. 'DIFF' MEANS THAT FOR / X95H AT LEAST ONE PAIR OF CORRESPONDING EIGENVALUES, THE MEMBERS OF X THE PAIR ARE NOT IDENTICAL. ) WRITE(*,2) 2 FORMAT(81H0-1.0 AS THE MEASURE OF PERFORMANCE IS PRINTED IF AN ER XROR IN THE PATH HAS / X58H PREVENTED THE COMPUTATION OF THE EIGENVALUES AND VECTORS. // X60H THE QZIT(T) PATH USES THE EISPACK CODES QZHES-QZIT-QZVAL X,8H-QZVEC, / X39H AS CALLED FROM DRIVER SUBROUTINE RGG. / X61H THE QZIT(F) PATH USES THE EISPACK CODES QZHES-QZIT-QZVAL, / X39H AS CALLED FROM DRIVER SUBROUTINE RGG. ) WRITE(*,15) 15 FORMAT(1X,21HS.P. VERSION 04/15/83 ) 5 FORMAT( 53H1 TABULATION OF THE ERROR FLAG ERROR AND THE , X 31HMEASURE OF PERFORMANCE Y FOR /5X, X 56HTHE EISPACK CODES. THIS RUN DISPLAYS THESE STATISTICS , X 21H FOR THE REAL GENERAL/4X,30H GENERALIZED EIGENPROBLEM A*X, X 15H= (LAMBDA)*B*X./ X 24H0ORDER QZIT(T) QZIT(F) ) 10 CALL RMATIN(NM,N,A,B,AHOLD,BHOLD,0) DO 510 ICALL = 1,2 IF( ICALL .NE. 1 ) CALL RMATIN(NM,N,A,B,AHOLD,BHOLD,1) GO TO (70,100), ICALL C C RGGWZ C INVOKED FROM DRIVER SUBROUTINE RGG. C 70 ICT = 1 CALL RGG(NM,N,A,B,ALFR,ALFI,BETA,1,Z,ERROR) IERR(ICT) = ERROR IF( ERROR .NE. 0 ) GO TO 500 GO TO 190 C C RGGW C INVOKED FROM DRIVER SUBROUTINE RGG. C 100 IF( IERR(1) .NE. 0 ) GO TO 510 CALL RGG(NM,N,A,B,ALFR1,ALFI1,BETA1,0,A,ERROR) IERR(1) = ERROR MATCH(1) = DIFF DO 101 I = 1,N IF( ALFR(I) .NE. ALFR1(I) .OR. ALFI(I) .NE. ALFI1(I) .OR. X BETA(I) .NE. BETA1(I) ) GO TO 510 101 CONTINUE MATCH(1) = SAME GO TO 510 C 190 CALL RMATIN(NM,N,A,B,AHOLD,BHOLD,1) CALL RGGWZR(NM,N,A,B,ALFR,ALFI,BETA,Z,NORM,RESDUL) DFL = 10 * N TCRIT(ICT) = RESDUL/EPSLON(DFL) GO TO 510 500 TCRIT(ICT) = -1.0E0 510 CONTINUE C IF( MOD(LCOUNT,35) .EQ. 0 ) WRITE(*,5) LCOUNT = LCOUNT + 1 WRITE(*,520) N,IERR(1),TCRIT(1),MATCH(1) 520 FORMAT(I4,I5,F6.3,4X,A4) GO TO 10 END SUBROUTINE RGGWZR(NM,N,A,B,ALFR,ALFI,BETA,Z,NORM,RESDUL) C REAL NORM(N),ALFR(N),ALFI(N),BETA(N),A(NM,N), X B(NM,N), Z(NM,N), TNORM, XR, XI, S, SUM, SUMZ, SUMR, SUMI, X SUMQ,RESDUL,NORMAB,SUM1,SUMA,SUMB,NORMA,NORMB, X SUMR2,SUMI2,SUMR3,SUMI3,PYTHAG LOGICAL CPLX C C THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX C A*Z-B*Z*DIAG(W) WHERE A AND B ARE REAL GENERAL MATRICES, Z IS C A MATRIX WHICH CONTAINS THE EIGENVECTORS OF THE EIGENPROBLEM C A*Z - B*Z*DIAG(W), AND W STANDS FOR A VECTOR OF CORRESPONDING C EIGENVALUES OF THE EIGENPROBLEM OBTAINED FROM THE VECTORS ALFR, C ALFI, AND BETA BY THE CORRESPONDENCES C W(J) = (ALFR(J) + I*ALFI(J)) / BETA(J) . C ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS. C C THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RGGWZR). C C INPUT. C C NM IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS C AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT; C C N IS THE ORDER OF THE MATRICES A AND B; C C A(NM,N) AND B(NM,N) ARE REAL GENERAL MATRICES; C C Z(NM,N) IS A MATRIX WHICH CONTAINS THE REAL AND IMAGINARY PARTS C OF THE EIGENVECTORS OF THE SYSTEM. THE I-TH COLUMN OF Z IS C A REAL EIGENVECTOR IF THE I-TH EIGENVALUE IS REAL. IF THE C I-TH EIGENVALUE IS COMPLEX AND THE FIRST OF A CONJUGATE C PAIR, THE I-TH AND (I+1)-TH COLUMNS OF Z CONTAIN THE REAL C AND IMAGINARY PARTS OF ITS EIGENVECTORS; C C ALFR(N), ALFI(N) AND BETA(N) ARE ARRAYS CONTAINING THE C COMPONENTS OF THE EIGENVALUES OF THE SYSTEM. C C OUTPUT. C C Z(NM,N) IS AN ARRAY WHICH CONTAINS THE NORMALIZED C APPROXIMATE EIGENVECTORS OF THE SYSTEM. THE EIGENVECTORS C ARE NORMALIZED USING THE 1-NORM IN SUCH A WAY THAT THE FIRST C ELEMENT WHOSE MAGNITUDE IS LARGER THAN THE NORM OF THE C EIGENVECTOR DIVIDED BY N IS REAL AND POSITIVE; C C NORM(N) IS AN ARRAY SUCH THAT FOR EACH K, C C !!BETA(K)*A*Z(K)-ALFA*B*Z(K)!! C NORM(K) = ----------------------------------------------- C !!Z(K)!!*(BETA(K)*!!A!!+!ALFA(K)!*!!B!!) C C WHERE Z(K) IS THE K-TH EIGENVECTOR AND ALFA = (ALFR + I*ALFI); C C RESDUL IS THE REAL NUMBER C !!BETA*A*Z-ALFA*B*Z!!/(!!Z!!*(BETA*!!A!!+!ALFA!*!!B!!)) C C ---------------------------------------------------------------- C NORMB = 0.0E0 NORMA = 0.0E0 RESDUL = 0.0E0 CPLX = .FALSE. C DO 20 I = 1,N SUMA = 0.0E0 SUMB = 0.0E0 DO 10 L = 1,N SUMA = SUMA + ABS(A(L,I)) 10 SUMB = SUMB + ABS(B(L,I)) NORMA = AMAX1(NORMA,SUMA) 20 NORMB = AMAX1(NORMB,SUMB) C DO 160 I=1,N IF(CPLX) GO TO 80 S = 0.0E0 SUMZ = 0.0E0 IF(ALFI(I) .NE. 0.0E0) GO TO 90 C ..........THIS IS FOR REAL EIGENVALUES.......... DO 40 L=1,N SUMQ = 0.0E0 SUM1 = 0.0E0 SUMZ = SUMZ + ABS(Z(L,I)) C DO 35 K=1,N SUMQ = SUMQ + B(L,K)*Z(K,I) 35 SUM1 = SUM1 + A(L,K)*Z(K,I) C SUM = -ALFR(I)*SUMQ 40 S = S + ABS(SUM + SUM1*BETA(I)) C NORMAB = NORMA*BETA(I) + NORMB*ABS(ALFR(I)) IF( NORMAB .EQ. 0.0E0 ) NORMAB = 1.0E0 NORM(I) = SUMZ IF( SUMZ .EQ. 0.0E0 ) GO TO 150 C ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE C WILL ALWAYS EXIST AN ELEMENT IN THE VECTOR Z(I) C LARGER THAN !!Z(I)!!/N.......... DO 50 L=1,N IF(ABS(Z(L,I)) .GE. NORM(I)/N) GO TO 60 50 CONTINUE C 60 TNORM = SIGN(NORM(I),Z(L,I)) C DO 70 L=1,N 70 Z(L,I) = Z(L,I)/TNORM C NORM(I) = S/(NORM(I)*NORMAB) 80 CPLX = .FALSE. GO TO 150 C ..........THIS IS FOR COMPLEX EIGENVALUES.......... 90 DO 110 L = 1,N SUMR = 0.0E0 SUMI = 0.0E0 SUMR2 = 0.0E0 SUMI2 = 0.0E0 SUMZ = SUMZ + PYTHAG(Z(L,I),Z(L,I+1)) C DO 100 K=1,N SUMR = SUMR + B(L,K)*Z(K,I) SUMI = SUMI + B(L,K)*Z(K,I+1) SUMR2 = SUMR2 + A(L,K)*Z(K,I) 100 SUMI2 = SUMI2 + A(L,K)*Z(K,I+1) SUMR3 = -ALFR(I)*SUMR + ALFI(I)*SUMI SUMI3 = -ALFI(I)*SUMR - ALFR(I)*SUMI C 110 S = S+PYTHAG(SUMR3+SUMR2*BETA(I),SUMI3+SUMI2*BETA(I)) NORMAB = NORMA*BETA(I)+NORMB*PYTHAG(ALFR(I),ALFI(I)) IF( NORMAB .EQ. 0.0E0 ) NORMAB = 1.0E0 C NORM(I) = SUMZ IF( SUMZ .EQ. 0.0E0 ) GO TO 150 C ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE WILL C ALWAYS EXIST AN ELEMENT IN THE VECTOR Z(I) LARGER C THAN !!Z(I)!!/N.......... DO 120 L=1,N IF(PYTHAG(Z(L,I),Z(L,I+1)) .GE. NORM(I)/N) 1 GO TO 130 120 CONTINUE C 130 XR = NORM(I)*Z(L,I)/PYTHAG(Z(L,I),Z(L,I+1)) XI = NORM(I)*Z(L,I+1)/PYTHAG(Z(L,I),Z(L,I+1)) C DO 140 L= 1,N CALL CDIV(Z(L,I),Z(L,I+1),XR,XI,Z(L,I),Z(L,I+1)) 140 CONTINUE C NORM(I) = S/(NORM(I)*NORMAB) NORM(I+1) = NORM(I) CPLX = .TRUE. 150 RESDUL = AMAX1(NORM(I),RESDUL) 160 CONTINUE C RETURN END SUBROUTINE RMATIN(NM,N,A,B,AHOLD,BHOLD,INITIL) C C THIS INPUT SUBROUTINE READS TWO REAL MATRICES A AND B FROM C SYSIN OF ORDER N. C TO GENERATE THE MATRICES INITIALLY, INITIL IS TO BE 0. C TO REGENERATE THE MATRICES FOR THE PURPOSE OF THE RESIDUAL C CALCULATION, INITIL IS TO BE 1. C C THIS ROUTINE IS CATALOGUED AS EISPDRV4(RSGREADI). C REAL A(NM,NM),B(NM,NM),AHOLD(NM,NM),BHOLD(NM,NM) INTEGER IA( 20), IB( 20) DATA IREADA/1/,IREADB/2/,IWRITE/6/ C IF( INITIL .EQ. 1 ) GO TO 30 READ(IREADA,5) N, M 5 FORMAT(I6,6X,I6) IF( N .EQ. 0 ) GO TO 70 IF (M .NE. 1) GO TO 16 DO 10 I = 1,N READ(IREADA,17) (IA(J), J=I,N) DO 9 J = I,N A(I,J) = IA(J) 9 A(J,I) = A(I,J) 10 CONTINUE 11 READ(IREADB,5) N,M IF( M .NE. 1 ) GO TO 20 DO 15 I = 1,N READ(IREADB,17) (IB(J), J=I,N) DO 14 J = I,N B(I,J)=IB(J) 14 B(J,I)=B(I,J) 15 CONTINUE GO TO 22 16 DO 18 I = 1,N READ(IREADA,17) (IA(J), J=1,N) 17 FORMAT(6I12) DO 18 J = 1,N 18 A(I,J) = IA(J) GO TO 11 20 DO 25 I = 1,N READ(IREADB,17) (IB(J),J=1,N) DO 25 J = 1,N 25 B(I,J) = IB(J) 22 DO 23 I = 1,N DO 23 J = 1,N BHOLD(I,J) = B(I,J) 23 AHOLD(I,J) = A(I,J) RETURN 30 DO 40 I = 1,N DO 40 J = 1,N B(I,J) = BHOLD(I,J) 40 A(I,J) = AHOLD(I,J) RETURN 70 WRITE(*,80) 80 FORMAT(47H0END OF DATA FOR SUBROUTINE RMATIN(RSGREADI). /1H1) STOP END