C C THIS DRIVER TESTS EISPACK FOR THE CLASS OF REAL GENERAL C MATRICES SUMMARIZING THE FIGURES OF MERIT FOR ALL PATHS. C C THIS DRIVER IS CATALOGUED AS EISPDRV4(RGSUMARY). C C THE DIMENSION OF A,Z,RM1, AND ASAVE SHOULD BE NM BY NM. C THE DIMENSION OF WR,WI,WR1,WI1,SELECT,SLHOLD,INT,SCALE,ORT,RV1, C AND RV2 SHOULD BE NM. C THE DIMENSION OF AHOLD SHOULD BE NM BY NM. C HERE NM = 20. C REAL A(20,20),Z(20,20),ASAVE(20,20),RM1(20,20), X AHOLD( 20, 20),ORT( 20),WR( 20),WI(20), X WR1( 20),WI1( 20),SCALE( 20),RV1( 20),RV2( 20), X TCRIT( 8),EPSLON,RESDUL,DFL INTEGER INT( 20),IERR( 8),MATCH( 4),LOW,UPP,ERROR LOGICAL SELECT( 20),SLHOLD( 20) CHARACTER*4 SAME,DIFF,ELEM,ORTH DATA IREAD1/1/,IREADC/5/,IWRITE/6/ DATA SAME,DIFF,ELEM,ORTH/'SAME','DIFF','E','O'/ C OPEN(IREAD1,FILE='FILE33') OPEN(IREADC,FILE='FILE34') REWIND IREAD1 REWIND IREADC C NM = 20 LCOUNT = 0 WRITE(*,1) 1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S XTATISTICS//1H ,95(1H-)/ 26X,9HBALANCING,15X,12HNO BALANCING,10X, X16HVECTORS COMPUTED/ 18X,3(24H------------------------,1X)/ X16H ORDER LOW UPP T,2X,2(25H HQR2 HQR INVIT )/ X1H ,95(1H-)//95H 'BALANCING' IS THE OPTION THAT EMPLOYS SUBROUTINE X BALANC TO BALANCE THE MATRIX BEFORE THE / X95H EIGENVALUES ARE COMPUTED AND BALBAK TO BACK TRANSFORM THE EI XGENVECTORS. // X95H 'VECTORS COMPUTED' IDENTIFIES BY POSITION, IN THE SET RETURNED X BY HQR, THOSE EIGENVALUES / X95H FOR WHICH INVIT COMPUTED THE ASSOCIATED EIGENVECTORS. T IN XDEXES AN EIGENVALUE FOR WHICH / X66H THE EIGENVECTOR WAS COMPUTED AND F INDEXES A PASSED EIGENVAL XUE. // 51H UNDER 'ORDER' IS THE ORDER OF EACH TEST MATRIX. // X95H UNDER 'LOW' AND 'UPP' ARE INTEGERS INDICATING THE BOUNDARY IND XICES FOR THE BALANCED MATRIX. // X95H UNDER 'T' IS THE LETTER E OR O INDICATING THE USE OF ELEMENTAR XY OR ORTHOGONAL TRANSFORMATIONS.// X95H UNDER 'HQR2 HQR' ARE TWO NUMBERS AND A KEYWORD. THE FIRST N XUMBER, AN INTEGER, IS THE ) WRITE(*,2) 2 FORMAT( X95H ABSOLUTE SUM OF THE ERROR FLAGS RETURNED SEPARATELY FROM HQR2 X AND HQR. THE SECOND / X95H NUMBER IS THE MEASURE OF PERFORMANCE BASED UPON THE RESIDUAL C XOMPUTED FOR THE HQR2 PATH. / X95H THE KEYWORD REPORTS THE DUPLICATION OF THE COMPUTED EIGENVALUE XS FROM HQR2 AND HQR. / X95H 'SAME' MEANS THAT THE EIGENVALUES ARE EXACT DUPLICATES. 'DIFF X' MEANS THAT FOR AT LEAST ONE / X95H PAIR OF CORRESPONDING EIGENVALUES, THE MEMBERS OF THE PAIR ARE X NOT IDENTICAL. // X62H UNDER 'INVIT' ARE TWO NUMBERS. THE FIRST NUMBER, AN INTEGER,, X27H IS THE ABSOLUTE SUM OF THE / X62H ERROR FLAGS RETURNED FROM THE PATH. THE SECOND NUMBER IS THE, X29H MEASURE OF PERFORMANCE BASED / X41H UPON THE RESIDUAL COMPUTED FOR THE PATH.// X62H -1.0 AS THE MEASURE OF PERFORMANCE IS PRINTED IF AN ERROR IN, X27H THE CORRESPONDING PATH HAS / X47H PREVENTED THE COMPUTATION OF THE EIGENVECTORS./) WRITE(*,3) 3 FORMAT(45H0FOR REDUCTIONS BY ELEMENTARY TRANSFORMATIONS/ X95H THE HQR2 PATH WITH BALANCING USES THE EISPACK CODES BAL XANC-ELMHES-ELTRAN-HQR2 -BALBAK,/ X38H AS CALLED FROM DRIVER SUBROUTINE RG. / X95H THE HQR PATH WITH BALANCING USES THE EISPACK CODES BAL XANC-ELMHES-HQR , / X38H AS CALLED FROM DRIVER SUBROUTINE RG. / X103H THE INVIT PATH WITH BALANCING USES THE EISPACK CODES BA XLANC-ELMHES-HQR -INVIT -ELMBAK-BALBAK. / X95H THE HQR2 PATH WITH NO BALANCING USES THE EISPACK CODES ELM XHES-ELTRAN-HQR2 . / X95H THE HQR PATH WITH NO BALANCING USES THE EISPACK CODES ELM XHES-HQR . / X95H THE INVIT PATH WITH NO BALANCING USES THE EISPACK CODES ELM XHES-HQR -INVIT -ELMBAK. ) WRITE(*,4) 4 FORMAT(45H0FOR REDUCTIONS BY ORTHOGONAL TRANSFORMATIONS/ X95H THE HQR2 PATH WITH BALANCING USES THE EISPACK CODES BAL XANC-ORTHES-ORTRAN-HQR2 -BALBAK./ X95H THE HQR PATH WITH BALANCING USES THE EISPACK CODES BAL XANC-ORTHES-HQR . / X103H THE INVIT PATH WITH BALANCING USES THE EISPACK CODES BA XLANC-ORTHES-HQR -INVIT -ORTBAK-BALBAK. / X95H THE HQR2 PATH WITH NO BALANCING USES THE EISPACK CODES ORT XHES-ORTRAN-HQR2 . / X95H THE HQR PATH WITH NO BALANCING USES THE EISPACK CODES ORT XHES-HQR . / X95H THE INVIT PATH WITH NO BALANCING USES THE EISPACK CODES ORT XHES-HQR -INVIT- ORTBAK. ) 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 30H FOR REAL GENERAL MATRICES. // X 26X,9HBALANCING,15X,12HNO BALANCING,10X,16HVECTORS SELECTED/ X 18X,3(24H------------------------,1X)/ X 16H ORDER LOW UPP T,2X,2(25H HQR2 HQR INVIT )) 10 CALL RMATIN(NM,N,A,AHOLD,0) READ(IREADC,50) MM,(SELECT(I),I=1,N) 50 FORMAT(I4/(72L1)) DO 60 I = 1,N 60 SLHOLD(I) = SELECT(I) C C MM AND SELECT ARE READ FROM SYSIN AFTER THE MATRIX IS C GENERATED. MM SPECIFIES TO INVIT THE MAXIMUM NUMBER OF C EIGENVECTORS THAT WILL BE COMPUTED. SELECT CONTAINS INFORMATION C ABOUT WHICH EIGENVECTORS ARE DESIRED. C DO 300 ICALL = 1,12 ICALL2 = MOD(ICALL,3) IF( ICALL2 .NE. 0 ) GO TO 80 DO 70 I = 1,N 70 SELECT(I) = SLHOLD(I) 80 IF( ICALL .NE. 1 ) CALL RMATIN(NM,N,A,AHOLD,1) GO TO (90,100,110,120,130,140,160,165,170,175,180,185), ICALL C C IF HQR PATH IS TAKEN THEN HQR2 PATH MUST ALSO BE TAKEN C IN ORDER THAT COMPARISON OF EIGENVALUES BE POSSIBLE. C C C RGEWZ C INVOKED FROM DRIVER SUBROUTINE RG. C 90 ICT = 1 CALL RG(NM,N,A,WR,WI,1,Z,INT,SCALE,ERROR) IERR(ICT) = ERROR IF( ERROR .NE. 0 ) GO TO 220 GO TO 190 C C RGEW C INVOKED FROM DRIVER SUBROUTINE RG. C 100 IF( IERR(1) .NE. 0 ) GO TO 300 CALL RG(NM,N,A,WR1,WI1,0,A,INT,SCALE,ERROR) IERR(1) = ERROR MATCH(1) = DIFF DO 101 I = 1,N IF( WR(I) .NE. WR1(I) .OR. WI(I) .NE. WI1(I) ) GO TO 300 101 CONTINUE MATCH(1) = SAME GO TO 300 C C RGEW1Z C 110 ICT = 2 CALL BALANC(NM,N,A,LOW,UPP,SCALE) CALL ELMHES(NM,N,LOW,UPP,A,INT) DO 112 I = 1,N DO 112 J = 1,N 112 ASAVE(I,J) = A(I,J) CALL HQR(NM,N,LOW,UPP,A,WR,WI,ERROR) IERR(ICT) = ERROR IF( ERROR .NE. 0 ) GO TO 220 CALL INVIT(NM,N,ASAVE,WR,WI,SELECT,MM,M,Z,ERROR,RM1,RV1,RV2) IERR(ICT) = IABS(ERROR) CALL ELMBAK(NM,LOW,UPP,ASAVE,INT,M,Z) CALL BALBAK(NM,N,LOW,UPP,SCALE,M,Z) GO TO 190 C C RGNEWZ C 120 ICT = 3 CALL ELMHES(NM,N,1,N,A,INT) CALL ELTRAN(NM,N,1,N,A,INT,Z) CALL HQR2(NM,N,1,N,A,WR,WI,Z,ERROR) IERR(ICT) = ERROR IF( ERROR .NE. 0 ) GO TO 220 DO 121 I = 1,N WR1(I) = WR(I) 121 WI1(I) = WI(I) GO TO 190 C C RGNEW C 130 IF( IERR(3) .NE. 0 ) GO TO 300 CALL ELMHES(NM,N,1,N,A,INT) CALL HQR(NM,N,1,N,A,WR,WI,ERROR) IERR(3) = ERROR MATCH(2) = DIFF DO 131 I = 1,N IF( WR(I) .NE. WR1(I) .OR. WI(I) .NE. WI1(I) ) GO TO 300 131 CONTINUE MATCH(2) = SAME GO TO 300 C C RGNEW1Z C 140 ICT = 4 CALL ELMHES(NM,N,1,N,A,INT) DO 142 I = 1,N DO 142 J = 1,N 142 ASAVE(I,J) = A(I,J) CALL HQR(NM,N,1,N,A,WR,WI,ERROR) IERR(ICT) = ERROR IF( ERROR .NE. 0 ) GO TO 220 CALL INVIT(NM,N,ASAVE,WR,WI,SELECT,MM,M,Z,ERROR,RM1,RV1,RV2) IERR(ICT) = IABS(ERROR) CALL ELMBAK(NM,1,N,ASAVE,INT,M,Z) GO TO 190 C C RGOWZ C 160 ICT = 5 CALL BALANC(NM,N,A,LOW,UPP,SCALE) CALL ORTHES(NM,N,LOW,UPP,A,ORT) CALL ORTRAN(NM,N,LOW,UPP,A,ORT,Z) CALL HQR2(NM,N,LOW,UPP,A,WR,WI,Z,ERROR) IERR(ICT) = ERROR IF( ERROR .NE. 0 ) GO TO 220 DO 161 I = 1,N WR1(I) = WR(I) 161 WI1(I) = WI(I) CALL BALBAK(NM,N,LOW,UPP,SCALE,N,Z) GO TO 190 C C RGOW C 165 IF( IERR(5) .NE. 0 ) GO TO 300 CALL BALANC(NM,N,A,LOW,UPP,SCALE) CALL ORTHES(NM,N,LOW,UPP,A,ORT) CALL HQR(NM,N,LOW,UPP,A,WR,WI,ERROR) IERR(5) = ERROR MATCH(3) = DIFF DO 166 I = 1,N IF( WR(I) .NE. WR1(I) .OR. WI(I) .NE. WI1(I) ) GO TO 300 166 CONTINUE MATCH(3) = SAME GO TO 300 C C RGOW1Z C 170 ICT = 6 CALL BALANC(NM,N,A,LOW,UPP,SCALE) CALL ORTHES(NM,N,LOW,UPP,A,ORT) DO 172 I = 1,N DO 172 J = 1,N 172 ASAVE(I,J) = A(I,J) CALL HQR(NM,N,LOW,UPP,A,WR,WI,ERROR) IERR(ICT) = ERROR IF( ERROR .NE. 0 ) GO TO 220 CALL INVIT(NM,N,ASAVE,WR,WI,SELECT,MM,M,Z,ERROR,RM1,RV1,RV2) IERR(ICT) = IABS(ERROR) CALL ORTBAK(NM,LOW,UPP,ASAVE,ORT,M,Z) CALL BALBAK(NM,N,LOW,UPP,SCALE,M,Z) GO TO 190 C C RGNOWZ C 175 ICT = 7 CALL ORTHES(NM,N,1,N,A,ORT) CALL ORTRAN(NM,N,1,N,A,ORT,Z) CALL HQR2(NM,N,1,N,A,WR,WI,Z,ERROR) IERR(ICT) = ERROR IF( ERROR .NE. 0 ) GO TO 220 DO 176 I = 1,N WR1(I) = WR(I) 176 WI1(I) = WI(I) GO TO 190 C C RGNOW C 180 IF( IERR(7) .NE. 0 ) GO TO 300 CALL ORTHES(NM,N,1,N,A,ORT) CALL HQR(NM,N,1,N,A,WR,WI,ERROR) IERR(7) = ERROR MATCH(4) = DIFF DO 181 I = 1,N IF( WR(I) .NE. WR1(I) .OR. WI(I) .NE. WI1(I) ) GO TO 300 181 CONTINUE MATCH(4) = SAME GO TO 300 C C RGNOW1Z C 185 ICT = 8 CALL ORTHES(NM,N,1,N,A,ORT) DO 187 I = 1,N DO 187 J = 1,N 187 ASAVE(I,J) = A(I,J) CALL HQR(NM,N,1,N,A,WR,WI,ERROR) IERR(ICT) = ERROR IF( ERROR .NE. 0 ) GO TO 220 CALL INVIT(NM,N,ASAVE,WR,WI,SELECT,MM,M,Z,ERROR,RM1,RV1,RV2) IERR(ICT) = IABS(ERROR) CALL ORTBAK(NM,1,N,ASAVE,ORT,M,Z) C 190 CALL RMATIN(NM,N,A,AHOLD,1) IF( ICALL2 .EQ. 0 ) CALL RGW1ZR(NM,N,A,WR,WI,SELECT, X M,Z,RV1,RESDUL) IF( ICALL2 .EQ. 1 ) CALL RGWZR(NM,N,A,WR,WI,Z,RV1,RESDUL) DFL = 10 * N TCRIT(ICT) = RESDUL/EPSLON(DFL) GO TO 300 220 TCRIT(ICT) = -1.0E0 300 CONTINUE C IF( MOD(LCOUNT,16) .EQ. 0 ) WRITE(*,5) LCOUNT = LCOUNT + 1 WRITE(*,520) X N,LOW,UPP,ELEM,(IERR(2*I-1),TCRIT(2*I-1),MATCH(I), X IERR(2*I),TCRIT(2*I),I=1,2),(SELECT(I),I=1,N) 520 FORMAT(I4,1X,2I4,2X,A1,2(I4,F6.3,1X,A4,I4,F6.3), X 4X,21L1/(72X,20L1)) WRITE(*,530) ORTH,(IERR(2*I-1),TCRIT(2*I-1),MATCH(I), X IERR(2*I),TCRIT(2*I),I=3,4) 530 FORMAT(15X,A1,2(I4,F6.3,1X,A4,I4,F6.3)) GO TO 10 END SUBROUTINE RGWZR(NM,N,A,WR,WI,Z,NORM,RESDUL) C REAL NORM(N),WR(N),WI(N),A(NM,N),Z(NM,N), X NORMA,TNORM,XR,XI,S,SUM,SUMA,SUMZ,SUMR,SUMI,RESDUL,PYTHAG LOGICAL CPLX C C THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX C A*Z-Z*DIAG(W) WHERE A IS A REAL GENERAL MATRIX, Z IS C A MATRIX WHICH CONTAINS THE EIGENVECTORS OF A, AND W C IS A VECTOR WHICH CONTAINS THE EIGENVALUES OF A. ALL C NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS. C C THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RGWZR). 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 MATRIX A; C C A(NM,N) IS A REAL GENERAL MATRIX; C C Z(NM,N) IS A MATRIX WHICH CONTAINS THE REAL AND IMAGINARY PARTS C OF THE EIGENVECTORS OF A. THE I-TH COLUMN OF Z IS A REAL C EIGENVECTOR IF THE I-TH EIGENVALUE IS REAL. IF THE I-TH C EIGENVALUE IS COMPLEX AND THE FIRST OF A CONJUGATE PAIR, THE C I-TH AND (I+1)-TH COLUMNS OF Z CONTAIN THE REAL AND C IMAGINARY PARTS OF ITS EIGENVECTORS; C C WR(N), WI(N) ARE ARRAYS CONTAINING THE REAL AND IMAGINARY C PARTS OF THE EIGENVALUES OF A. C C OUTPUT. C C Z(NM,N) IS AN ARRAY WHICH CONTAINS THE NORMALIZED C APPROXIMATE EIGENVECTORS OF A. THE EIGENVECTORS ARE C 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 NORM(K) = !!A*Z(K)-Z(K)*W(K)!!/(!!A!!*!!Z(K)!!) C WHERE Z(K) IS THE K-TH EIGENVECTOR; C C RESDUL IS THE REAL NUMBER C !!A*Z-Z*DIAG(W)!!/(!!A!!*!!Z!!). C C ---------------------------------------------------------------- C NORMA = 0.0E0 RESDUL = 0.0E0 CPLX = .FALSE. C DO 20 I=1,N SUMA = 0.0E0 C DO 10 L=1,N 10 SUMA = SUMA + ABS(A(L,I)) C 20 NORMA = AMAX1(NORMA,SUMA) C IF(NORMA .EQ. 0.0E0) NORMA = 1.0E0 C DO 160 I=1,N IF(CPLX) GO TO 80 S = 0.0E0 SUMZ = 0.0E0 IF(WI(I) .NE. 0.0E0) GO TO 90 C ..........THIS IS FOR REAL EIGENVALUES.......... DO 40 L=1,N SUMZ = SUMZ + ABS(Z(L,I)) SUM = -WR(I)*Z(L,I) C DO 30 K=1,N 30 SUM = SUM + A(L,K)*Z(K,I) C 40 S = S + ABS(SUM) C NORM(I) = SUMZ 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)*NORMA) 80 CPLX = .FALSE. GO TO 150 C ..........THIS IS FOR COMPLEX EIGENVALUES.......... 90 DO 110 L = 1,N SUMZ = SUMZ + PYTHAG(Z(L,I),Z(L,I+1)) SUMR = -WR(I)*Z(L,I) + WI(I)*Z(L,I+1) SUMI = -WR(I)*Z(L,I+1) - WI(I)*Z(L,I) C DO 100 K=1,N SUMR = SUMR + A(L,K)*Z(K,I) 100 SUMI = SUMI + A(L,K)*Z(K,I+1) C 110 S = S + PYTHAG(SUMR,SUMI) C NORM(I) = SUMZ 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 130 CONTINUE C XR = PYTHAG(Z(L,I),Z(L,I+1)) 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)*NORMA) NORM(I+1) = NORM(I) CPLX = .TRUE. 150 RESDUL = AMAX1(NORM(I),RESDUL) 160 CONTINUE C RETURN END SUBROUTINE RGW1ZR(NM,N,A,WR,WI,SELECT,M,Z,NORM,RESDUL) C REAL NORM(N),WR(N),WI(N),A(NM,N),Z(NM,M), X NORMA,TNORM,XR,XI,S,SUM,SUMA,SUMZ,SUMR,SUMI,RESDUL,PYTHAG LOGICAL SELECT(N), CPLX C C THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX C A*Z-Z*DIAG(W) WHERE A IS A REAL GENERAL MATRIX, Z IS C A MATRIX WHICH CONTAINS SOME OF THE EIGENVECTORS OF A, C AND W IS A VECTOR WHICH CONTAINS THE EIGENVALUES OF A. C ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS. C C THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RGW1ZR). 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 MATRIX A. C C A(NM,N) IS A REAL GENERAL MATRIX. C C Z(NM,M) IS A MATRIX WHICH CONTAINS THE REAL AND IMAGINARY C OF THE EIGENVECTORS OF A. THE EIGENVECTORS ARE STORED C PARTS IN PACKED FORM, THAT IS, USING ONE COLUMN FOR EACH C REAL EIGENVECTOR AND USING TWO COLUMNS FOR EACH COMPLEX C EIGENVECTOR, THE EIGENVECTORS CORRESPONDING TO THE FLAGGED C EIGENVALUES ARE STORED IN CONSECUTIVE COLUMNS, STARTING IN C COLUMN 1. (NOTE. FOR THE DEFINITION OF FLAGGED EIGENVALUES C SEE THE DESCRIPTION OF THE SELECT ARRAY.) C C WR(N),WI(N) ARE ARRAYS CONTAINING THE REAL AND IMAGINARY C PARTS OF THE EIGENVALUES OF A. C C SELECT(N) IS A LOGICAL ARRAY WHICH SPECIFIES THE C FLAGGED EIGENVALUES WHOSE EIGENVECTORS ARE C CONTAINED IN THE MATRIX Z. THE I-TH EIGENVALUE C (WR(I),WI(I)) IS FLAGGED IF SELECT(I) = .TRUE. C C M IS THE NUMBER OF COLUMNS OF Z CONTAINING EIGENVECTORS C OF A. C C OUTPUT. C C Z(NM,M) IS AN ARRAY WHICH CONTAINS THE NORMALIZED C APPROXIMATE EIGENVECTORS OF A IN PACKED FORM AS C DESCRIBED ABOVE. THE EIGENVECTORS ARE NORMALIZED C USING THE 1-NORM IN SUCH A WAY THAT THE FIRST ELEMENT C 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 WHENEVER SELECT(K) = .TRUE. C NORM(K) = !!A*Z(K)-Z(K)*W(K)!!/(!!A!!*!!Z(K)!!) C WHERE Z(K) IS THE EIGENVECTOR CORRESPONDING TO THE C K-TH EIGENVALUE (WR(K),WI(K)). C C RESDUL IS THE REAL NUMBER C !!A*Z-Z*DIAG(W)!!/(!!A!!*!!Z!!) C C ---------------------------------------------------------------- C KCOM = 0 NORMA = 0.0E0 RESDUL = 0.0E0 IF( M .EQ. 0 ) RETURN CPLX = .FALSE. II=1 C DO 20 I=1,N SUMA = 0.0E0 C DO 10 L=1,N 10 SUMA = SUMA + ABS(A(L,I)) C 20 NORMA = AMAX1(NORMA,SUMA) C IF(NORMA .EQ. 0.0E0) NORMA = 1.0E0 C DO 160 I=1,N IF(CPLX) GO TO 80 S = 0.0E0 SUMZ = 0.0E0 IF(WI(I) .NE. 0.0E0) GO TO 90 IF( .NOT. SELECT(I)) GO TO 160 C ..........THIS IS FOR REAL EIGENVALUES.......... DO 40 L=1,N SUMZ = SUMZ + ABS(Z(L,II)) SUM = -WR(I)*Z(L,II) C DO 30 K=1,N 30 SUM = SUM + A(L,K)*Z(K,II) C 40 S = S + ABS(SUM) C IF(SUMZ .NE. 0.0E0) GO TO 45 NORM(I) = -1.0E0 GO TO 150 45 NORM(I) = SUMZ C ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE C WILL ALWAYS EXIST AN ELEMENT IN THE VECTOR Z(II) C LARGER THAN !!Z(II)!!/N.......... DO 50 L=1,N IF(ABS(Z(L,II)) .GE. NORM(I)/N) GO TO 60 50 CONTINUE C 60 TNORM = SIGN(NORM(I),Z(L,II)) C DO 70 J=1,N 70 Z(J,II) = Z(J,II)/TNORM C NORM(I) = S/(NORM(I)*NORMA) GO TO 150 80 CPLX = .FALSE. GO TO 155 C ..........THIS IS FOR COMPLEX EIGENVALUES.......... 90 IF(SELECT(I)) GO TO 95 KCOM = MOD(KCOM+1,2) GO TO 160 C 95 DO 110 L = 1,N SUMZ = SUMZ + PYTHAG(Z(L,II),Z(L,II+1)) SUMR = -WR(I)*Z(L,II) + WI(I)*Z(L,II+1) SUMI = -WR(I)*Z(L,II+1) - WI(I)*Z(L,II) C DO 100 K=1,N SUMR = SUMR + A(L,K)*Z(K,II) 100 SUMI = SUMI + A(L,K)*Z(K,II+1) C 110 S = S + PYTHAG(SUMR,SUMI) C IF(SUMZ .NE. 0.0E0) GO TO 115 NORM(I) = -1.0E0 GO TO 145 115 NORM(I) = SUMZ C ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE WILL C ALWAYS EXIST AN ELEMENT IN THE VECTOR Z(II) LARGER C THAN !!Z(II)!!/N.......... DO 120 L=1,N IF(PYTHAG(Z(L,II),Z(L,II+1)) .GE. NORM(I)/N) 1 GO TO 130 120 CONTINUE C 130 XR = NORM(I)*Z(L,II)/PYTHAG(Z(L,II),Z(L,II+1)) XI = NORM(I)*Z(L,II+1)/PYTHAG(Z(L,II),Z(L,II+1)) C DO 140 J= 1,N CALL CDIV(Z(J,II),Z(J,II+1),XR,XI,Z(J,II),Z(J,II+1)) 140 CONTINUE C NORM(I) = S/(NORM(I)*NORMA) 145 CPLX = .TRUE. IF(KCOM .EQ. 0) GO TO 150 CPLX = .FALSE. II = II+1 150 RESDUL = AMAX1(NORM(I),RESDUL) 155 II = II+1 160 CONTINUE C RETURN END SUBROUTINE RMATIN(NM,N,A,AHOLD,INITIL) C C THIS INPUT SUBROUTINE READS A REAL MATRIX FROM SYSIN OF C ORDER N. C TO GENERATE THE MATRIX A INITIALLY, INITIL IS TO BE 0. C TO REGENERATE THE MATRIX A FOR THE PURPOSE OF THE RESIDUAL C CALCULATION, INITIL IS TO BE 1. C C THIS ROUTINE IS CATALOGUED AS EISPDRV4(RGREADI). C REAL A(NM,NM),AHOLD(NM,NM) INTEGER IA( 20) DATA IREADA/1/,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 14 DO 13 I = 1,N READ(IREADA,16) (IA(J), J=I,N) DO 12 J = I,N A(I,J) = IA(J) 12 A(J,I) = A(I,J) 13 CONTINUE GO TO 19 14 DO 17 I = 1,N READ(IREADA,16) (IA(J), J=1,N) 16 FORMAT(6I12) DO 17 J = 1,N 17 A(I,J) = IA(J) 19 DO 20 I = 1,N DO 20 J = 1,N 20 AHOLD(I,J) = A(I,J) RETURN 30 DO 40 I = 1,N DO 40 J = 1,N 40 A(I,J) = AHOLD(I,J) RETURN 70 WRITE(*,80) 80 FORMAT(46H0END OF DATA FOR SUBROUTINE RMATIN(RGREADI). /1H1) STOP END