C C THIS DRIVER TESTS EISPACK FOR THE CLASS OF REAL SYMMETRIC C BAND MATRICES SUMMARIZING THE FIGURES OF MERIT FOR ALL PATHS. C C THIS DRIVER IS CATALOGUED AS EISPDRV4(RSBSUMAR). C C THE DIMENSION OF ST AND AA SHOULD BE NM BY MMB. C THE DIMENSION OF Z SHOULD BE NM BY NM. C THE DIMENSION OF W,W1,W2,D,E,E2,RV3,RV4,RV5,RV6, C AND IND SHOULD BE NM. C THE DIMENSION OF STHOLD SHOULD BE NM BY MMB. C THE DIMENSION OF RV1 SHOULD BE 2*MMB**2+4*MMB-3. C THE DIMENSION OF RV2 SHOULD BE NM*(2*MMB-1). C HERE NM = 44 AND MMB = 5. C DOUBLE PRECISION ST( 44, 5),Z( 44, 44),AA( 44, 5), X STHOLD( 44, 5),W( 44),D( 44),E( 44), X E2( 44),W1( 44),W2( 44),RV1( 67),RV2( 396), X RV4( 44),RV5( 44),RV6( 44),TCRIT( 8),EPSLON,RESDUL, X MAXDIF,MAXEIG,U,LB,UB,EPS1,T,R,TSTART,DFL REAL XUB,XLB INTEGER IND( 44),IERR( 6),ERROR DATA IREAD1/1/,IREADC/5/,IWRITE/6/ C OPEN(UNIT=IREAD1,FILE='FILE37') OPEN(UNIT=IREADC,FILE='FILE38') REWIND IREAD1 REWIND IREADC C NM = 44 MMB = 5 NV1 = 67 NV2 = 396 LCOUNT = 0 WRITE(IWRITE,1) 1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S XTATISTICS//1H ,95(1H-)/ 96H ORDER B.W. TQL2 TQLRAT IMTQL2 IMTQL X1 LB UB M IMTQLV BISECT M1 NO TRIDIB BQR /1H , X95(1H-)//48H UNDER 'ORDER' IS THE ORDER OF EACH TEST MATRIX. // X51H UNDER 'B.W.' IS THE HALF BAND WIDTH OF THE MATRIX. // X95H UNDER 'TQL2 TQLRAT' ARE THREE NUMBERS. THE FIRST NUMBER, AN X INTEGER, IS THE ABSOLUTE SUM OF/ X61H THE ERROR FLAGS RETURNED SEPARATELY FROM TQL2 AND TQLRAT., X34H THE SECOND NUMBER IS THE MEASURE/ X62H OF PERFORMANCE BASED UPON THE RESIDUAL COMPUTED FOR THE TQL2, X25H PATH. THE THIRD NUMBER / X62H MEASURES THE AGREEMENT OF THE EIGENVALUES FROM THE TQL2 AND, X16H TQLRAT PATHS. // X95H UNDER 'IMTQL2 IMTQL1' ARE THREE NUMBERS WITH MEANING LIKE THOS XE UNDER 'TQL2 TQLRAT'. // X95H UNDER 'LB' AND 'UB' ARE THE INPUT VARIABLES SPECIFYING THE INT XERVAL TO BISECT. // X61H UNDER 'M' IS THE NUMBER OF EIGENVALUES DETERMINED BY BISECT, X17H THAT LIE IN THE /18H INTERVAL (LB,UB)./) WRITE(IWRITE,2) 2 FORMAT( X62H UNDER EACH OF 'IMTQLV', 'BISECT', 'TRIDIB', AND 'BQR' ARE TWO, X28H NUMBERS. THE FIRST NUMBER, / X95H AN INTEGER, IS THE ABSOLUTE SUM OF THE ERROR FLAGS RETURNED FR XOM THE RESPECTIVE PATH. / X95H THE SECOND NUMBER IS THE MEASURE OF PERFORMANCE BASED UPON THE X RESIDUAL COMPUTED FOR THE PATH.// X95H UNDER 'M1' AND 'NO' ARE THE VARIABLES SPECIFYING THE LOWER BOU XNDARY INDEX AND THE NUMBER / X28H OF EIGENVALUES TO TRIDIB. // 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. // X62H THE TQL2 PATH USES THE EISPACK CODES BANDR-TQL2 , / X39H AS CALLED FROM DRIVER SUBROUTINE RSB. / X62H THE TQLRAT PATH USES THE EISPACK CODES BANDR-TQLRAT, / X39H AS CALLED FROM DRIVER SUBROUTINE RSB. / X62H THE IMTQL2 PATH USES THE EISPACK CODES BANDR-IMTQL2. ) WRITE(IWRITE,3) 3 FORMAT( X62H THE IMTQL1 PATH USES THE EISPACK CODES BANDR-IMTQL1. / X64H THE IMTQLV PATH USES THE EISPACK CODES BANDR-IMTQLV-BANDV. X / X64H THE BISECT PATH USES THE EISPACK CODES BANDR-BISECT-BANDV. X / X64H THE TRIDIB PATH USES THE EISPACK CODES BANDR-TRIDIB-BANDV. X / X62H THE BQR PATH USES THE EISPACK CODES BQR -BANDV. ) WRITE(IWRITE,15) 15 FORMAT(1X,21HD.P. VERSION 04/15/83 ) 5 FORMAT( 53H1 TABULATION OF THE ERROR FLAG ERROR AND THE , X 30HMEASURE OF PERFORMANCE Y FOR /5X, X 56HTHE EISPACK CODES. THIS RUN DISPLAYS THESE STATISTICS , X 38H FOR REAL SYMMETRIC BAND MATRICES. / X 58H0ORDER B.W. TQL2 TQLRAT IMTQL2 IMTQL1 LB UB M, X 38H IMTQLV BISECT M1 NO TRIDIB BQR ) 10 CALL RMATIN(NM,MMB,N,MB,ST,STHOLD,0) READ(IREADC,50) MM,LB,UB,M11,NO,MBQR,TSTART 50 FORMAT(I4,2D24.16,3I4,F8.4) C C MM,LB,UB,M11,NO,MBQR, AND TSTART ARE READ FROM SYSIN AFTER THE C MATRIX IS GENERATED. MM,LB, AND UB SPECIFY TO BISECT THE C MAXIMUM NUMBER OF EIGENVALUES AND THE BOUNDS FOR THE INTERVAL C WHICH IS TO BE SEARCHED. M11 AND NO SPECIFY TO TRIDIB THE C LOWER BOUNDARY INDEX AND THE NUMBER OF DESIRED EIGENVALUES. C MBQR AND TSTART SPECIFY TO BQR THE NUMBER OF EIGENVALUES C AND THE VALUE TO WHICH THEY ARE DESIRED CLOSEST. C DO 230 ICALL = 1,10 IF( ICALL .NE. 1 ) CALL RMATIN(NM,MMB,N,MB,ST,STHOLD,1) C C IF TQLRAT PATH (LABEL 80) IS TAKEN THEN TQL2 PATH (LABEL 70) C MUST ALSO BE TAKEN IN ORDER THAT THE MEASURE OF PERFORMANCE BE C MEANINGFUL. C IF IMTQL1 PATH (LABEL 85) IS TAKEN THEN IMTQL2 PATH (LABEL 75) C MUST ALSO BE TAKEN IN ORDER THAT THE MEASURE OF PERFORMANCE BE C MEANINGFUL. C IF TQL2 (IMTQL2) PATH FAILS, THEN TQLRAT (IMTQL1) PATH IS C OMITTED AND PRINTOUT FLAGGED WITH -1.0. C GO TO (70,75,80,85,89,230,230,120,125,130), ICALL C C RSBWZ USING TQL2 C INVOKED FROM DRIVER SUBROUTINE RSB. C 70 ICT = 1 CALL RSB(NM,N,MB,ST,W,1,Z,E,E,ERROR) IERR(ICT) = ERROR M = ERROR - 1 IF( ERROR .NE. 0 ) GO TO 200 DO 71 I = 1,N W1(I) = W(I) 71 CONTINUE M = N GO TO 195 C C RSBWZ USING IMTQL2 C 75 ICT = 2 CALL BANDR(NM,N,MB,ST,W,E,E,.TRUE.,Z) CALL IMTQL2(NM,N,W,E,Z,ERROR) IERR(ICT) = ERROR M = ERROR - 1 IF( ERROR .NE. 0 ) GO TO 200 DO 78 I=1,N 78 W2(I) = W(I) M = N GO TO 195 C C RSBW USING TQLRAT C INVOKED FROM DRIVER SUBROUTINE RSB. C 80 ICT = 7 IF( IERR(1) .NE. 0 ) GO TO 200 CALL RSB(NM,N,MB,ST,W,0,Z,E,E2,ERROR) IF( ERROR .NE. 0 ) GO TO 200 IERR(1) = ERROR MAXEIG = 0.0D0 MAXDIF = 0.0D0 DO 81 I = 1,N IF( DABS(W(I)) .GT. MAXEIG ) MAXEIG = DABS(W(I)) U = DABS(W1(I) - W(I)) IF( U .GT. MAXDIF ) MAXDIF = U 81 CONTINUE IF( MAXEIG .EQ. 0.0D0 ) MAXEIG = 1.0D0 DFL = 10 * N TCRIT(ICT) = MAXDIF/EPSLON(MAXEIG*DFL) GO TO 230 C C RSBW USING IMTQL1 C 85 ICT = 8 IF( IERR(2) .NE. 0 ) GO TO 200 CALL BANDR(NM,N,MB,ST,W,E,E,.FALSE.,Z) CALL IMTQL1(N,W,E,ERROR) IERR(2) = ERROR MAXEIG = 0.0D0 MAXDIF = 0.0D0 DO 86 I = 1,N IF( DABS(W(I)) .GT. MAXEIG ) MAXEIG = DABS(W(I)) U = DABS(W2(I) - W(I)) IF( U .GT. MAXDIF ) MAXDIF = U 86 CONTINUE IF( MAXEIG .EQ. 0.0D0 ) MAXEIG = 1.0D0 DFL = 10 * N TCRIT(ICT) = MAXDIF/EPSLON(MAXEIG*DFL) GO TO 230 C C RSBW1Z ( USAGE HERE COMPUTES ALL THE EIGENVECTORS ) C 89 ICT = 3 DO 894 I = 1,MB DO 894 J = 1,N 894 AA(J,I) = ST(J,I) CALL BANDR(NM,N,MB,ST,D,E,E2,.FALSE.,Z) CALL IMTQLV(N,D,E,E2,W,IND,ERROR,RV1) IERR(ICT) = ERROR M = N IF( ERROR .NE. 0 ) M = ERROR - 1 CALL BANDV(NM,N,MB,AA,0.0D0,M,W,Z,ERROR,NV2,RV2,RV6) IERR(ICT) = IERR(ICT) + IABS(ERROR) GO TO 195 C C RSB1W1Z USING BISECT AND BANDV C 120 ICT = 4 EPS1 = 0.0D0 DO 122 I=1,MB DO 122 J=1,N 122 AA(J,I) = ST(J,I) CALL BANDR(NM,N,MB,ST,D,E,E2,.FALSE.,Z) CALL BISECT(N,EPS1,D,E,E2,LB,UB,MM,M,W,IND,ERROR,RV4,RV5) MBISCT = M XLB = LB XUB = UB IERR(ICT) = ERROR IF( ERROR .NE. 0 ) GO TO 200 CALL BANDV(NM,N,MB,AA,0.0D0,M,W,Z,ERROR,NV2,RV2,RV6) IERR(ICT) = IABS(ERROR) GO TO 195 C C RSB1W1Z USING TRIDIB AND BANDV C 125 ICT = 5 EPS1 = 0.0D0 DO 126 I=1,MB DO 126 J=1,N 126 AA(J,I) = ST(J,I) CALL BANDR(NM,N,MB,ST,D,E,E2,.FALSE.,Z) CALL TRIDIB(N,EPS1,D,E,E2,LB,UB,M11,NO,W,IND,ERROR,RV4,RV5) IERR(ICT) = ERROR IF( ERROR .NE. 0 ) GO TO 200 M = NO CALL BANDV(NM,N,MB,AA,0.0D0,M,W,Z,ERROR,NV2,RV2,RV6) IERR(ICT) = IABS(ERROR) GO TO 195 C C RSB1W1Z USING BQR AND BANDV C 130 ICT = 6 DO 132 I = 1,MB DO 132 J = 1,N 132 AA(J,I) = ST(J,I) N1 = N T = TSTART R = 0.0D0 DO 133 J = 1,N 133 ST(J,MB) = ST(J,MB) - T DO 138 I = 1,MBQR CALL BQR(NM,N1,MB,ST,T,R,ERROR,NV1,RV1) IF( ERROR .EQ. 0 ) GO TO 134 M = I - 1 GO TO 139 134 IF( I .EQ. 1 ) GO TO 136 DO 135 JJ = 2,I J = I + 2 - JJ IF( T .GT. W(J - 1) ) GO TO 137 W(J) = W(J - 1) 135 CONTINUE 136 J = 1 137 W(J) = T N1 = N1 - 1 138 CONTINUE M = MBQR 139 IERR(ICT) = ERROR CALL BANDV(NM,N,MB,AA,0.0D0,M,W,Z,ERROR,NV2,RV2,RV6) IERR(ICT) = IERR(ICT) + IABS(ERROR) GO TO 195 C 195 CALL RMATIN(NM,MMB,N,MB,ST,STHOLD,1) CALL RSBWZR(NM,N,M,MB,ST,W,Z,RV1,RESDUL) DFL = 10 * N TCRIT(ICT) = RESDUL/EPSLON(DFL) GO TO 230 200 TCRIT(ICT) = -1.0D0 230 CONTINUE C IF( MOD(LCOUNT,35) .EQ. 0 ) WRITE(IWRITE,5) LCOUNT = LCOUNT + 1 WRITE(IWRITE,520) N,MB,IERR(1),TCRIT(1),TCRIT(7),IERR(2),TCRIT(2), X TCRIT(8),XLB,XUB,MBISCT,(IERR(I),TCRIT(I),I=3,4), X M11,NO,(IERR(I),TCRIT(I),I=5,6) 520 FORMAT(I4,2X,I3,2(I3,2F6.3),2(1PE8.0),I3,2(I3,0PF6.3), X I2,1X,I2,2(I3,0PF6.3)) GO TO 10 END SUBROUTINE RSBWZR(NM,N,M,MB,A,W,Z,NORM,RESDUL) C INTEGER MB,MB1 DOUBLE PRECISION NORM(M), W(M), A(NM,MB), Z(NM,M), NORMA, TNORM, X S, SUM, SUMA, SUMZ, RESDUL C C THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX C A*Z-Z*DIAG(W) WHERE A IS A REAL SYMMETRIC BAND MATRIX, C W IS A VECTOR WHICH CONTAINS M EIGENVALUES OF A, AND Z C IS AN ARRAY WHICH CONTAINS THE M CORRESPONDING EIGENVECTORS OF C A. ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS. C C THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RSBWZR). 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 M IS THE NUMBERS OF EIGENVECTORS WHOSE RESIDUALS ARE DESIRED; C C MB IS THE BAND WIDTH OF THE INPUT MATRIX A . BAND WIDTH IS C DEFINED AS THE NUMBER OF ADJACENT DIAGONALS, INCLUDING THE C PRINCIPAL DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO C PORTION OF THE LOWER TRIANGLE OF THE MATRIX; C C A(N,MB) IS AN ARRAY WHICH CONTAINS IN ITS COLUMNS THE C SUBDIAGONAL AND DIAGONAL OF THE SYMMETRIC BAND C MATRIX; C C W(M) IS A VECTOR WHOSE FIRST M COMPONENTS CONTAIN EIGENVALUES C OF A; C C Z(NM,M) IS AN ARRAY WHOSE FIRST M COLUMNS CONTAIN THE C EIGENVECTORS OF A CORRESPONDING TO THE EIGENVALUES IN W. C C OUTPUT. C C Z(NM,M) IS AN ARRAY WHICH CONTAINS THE NORMALIZED C APPROXIMATE EIGENVECTORS OF A. THE EIGENVECTORS C ARE NORMALIZED USING THE 1-NORM IN SUCH A WAY C THAT THE FIRST ELEMENT WHOSE MAGNITUDE IS LARGER C THAN THE NORM OF THE EIGENVECTOR DIVIDED BY N IS C POSITIVE; C C NORM(M) 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 RESDUL = 0.0D0 IF( M .EQ. 0 ) RETURN NORMA = 0.0D0 C DO 40 I=1,N J=I SUMA = 0.0D0 IF(I .EQ. 1) GO TO 20 MB1 = MB-1 LSTART = MAX0(1,MB+1-I) C DO 10 L=LSTART,MB1 10 SUMA = SUMA + DABS(A(I,L)) C 20 LSTOP = MIN0(MB,N+1-J) C DO 30 L=1,LSTOP L1 = MB + 1 - L SUMA = SUMA + DABS(A(J,L1)) 30 J = J+1 C 40 NORMA = DMAX1(SUMA,NORMA) C IF(NORMA .EQ. 0.0D0) NORMA = 1.0D0 C DO 120 I=1,M S = 0.0D0 SUMZ = 0.0D0 C DO 80 L=1,N SUM = -W(I)*Z(L,I) SUMZ = SUMZ + DABS(Z(L,I)) J = MAX0(0,L-MB) IF(L .EQ. 1) GO TO 60 MB1 = MB-1 KSTART = MAX0(1,MB+1-L) C DO 50 K=KSTART,MB1 J = J+1 50 SUM = SUM + A(L,K)*Z(J,I) C 60 KSTOP = MIN0(MB,N+1-L) C DO 70 K=1,KSTOP J = J+1 K1 = MB + 1 - K 70 SUM = SUM + A(J,K1)*Z(J,I) C 80 S = S + DABS(SUM) C NORM(I) = SUMZ IF (SUMZ .EQ. 0.0D0) GO TO 120 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 90 L=1,N IF(DABS(Z(L,I)) .GE. NORM(I)/N) GO TO 100 90 CONTINUE C 100 TNORM = DSIGN(NORM(I),Z(L,I)) C DO 110 L=1,N 110 Z(L,I) = Z(L,I)/TNORM C NORM(I) = S/(NORM(I)*NORMA) 120 RESDUL = DMAX1(NORM(I),RESDUL) C RETURN END SUBROUTINE RMATIN(NM,MMB,N,MB,ST,STHOLD,INITIL) C C THIS INPUT SUBROUTINE READS A REAL SYMMETRIC BAND MATRIX C FROM SYSIN OF ORDER N, AND BAND WIDTH MB . C TO GENERATE THE MATRIX ST INITIALLY, INITIL IS TO BE 0. C TO REGENERATE THE MATRIX ST FOR THE PURPOSE OF THE RESIDUAL C CALCULATION, INITIL IS TO BE 1. C C THIS ROUTINE IS CATALOGUED AS EISPDRV4(RSBREADI). C DOUBLE PRECISION ST(NM,MMB),STHOLD(NM,MMB) INTEGER IA( 5) DATA IREADA/1/,IWRITE/6/ C IF( INITIL .EQ. 1 ) GO TO 30 READ(IREADA,5) N,MB 5 FORMAT(2I6) IF( N .EQ. 0 ) GO TO 70 DO 8 I = 1,N DO 7 J = 1,MB ST(I,J) = 0.0D0 7 CONTINUE 8 CONTINUE DO 15 I=1,N MBB = MIN0(MB,N-I+1) READ(IREADA,10) (IA(J),J=1,MBB) 10 FORMAT(6I12) DO 11 J=1,MBB M = MB+1-J K = I+J-1 11 ST(K,M) = IA(J) 15 CONTINUE DO 20 I = 1,N DO 20 J = 1,MB 20 STHOLD(I,J) = ST(I,J) RETURN 30 DO 40 I = 1,N DO 40 J = 1,MB 40 ST(I,J) = STHOLD(I,J) RETURN 70 WRITE(IWRITE,80) 80 FORMAT(46H0END OF DATA FOR SUBROUTINE RMATIN(RSBREADI). /1H1) STOP END