PROGRAM DC IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NMAX=160,N0=2000) REAL*8 M DIMENSION U(3,NMAX,NMAX),M(3,NMAX,NMAX), 2NS(NMAX),XT(3,NMAX),OMEGA(NMAX,NMAX),RRV(3,N0) real*8 E(N0),UM(N0,3),UA(N0,3),R(N0,3), 1COEF(N0,N0),EVAL(N0),DM(N0,3),AMM(N0,3),AMMI(N0,3), 1V(N0,N0) CHARACTER*1 OK logical*4 lex pi=3.141592265d0 bohr=0.529177d0 debye=2.541765d0 cm=219470.0d0 WRITE(6,66666)2005 66666 FORMAT(80(1H-),/,15X,'TRANSITION DIPOLE COUPLING',/,80(1H-),/,/, 25X,I10, ' version, UOCHB, Prague',/,/, 35X,'INPUT : DC.SC ... dipole definition',/, 35X,' DC.PAR ... options (Optional for uvcd)',/, 35X,'OUTPUT: DC.TAB ... spectrum',/, 55X,' DIP.MOM ... dipole moments (optional)',/, 65X,' EIGEN.VEC ... eigenvectors',/) c c Read DC.SC: OPEN(22,FILE='DC.SC',STATUS='OLD') READ(22,*)NG if(NG.gt.NMAX)call report('NG>NMAX') DO 1 I=1,NG READ(22,*) READ(22,*)(XT(j,I),j=1,3) READ(22,*)NS(I) if(NS(I).gt.NMAX)call report('NS>NMAX') READ(22,*) READ(22,*)IDUMM,OMEGA(I,1) DO 2 J=2,NS(I) READ(22,*)IDUMM,OMEGA(I,J) 2 OMEGA(I,J)=OMEGA(I,J)-OMEGA(I,1) READ(22,*) DO 1 J=2,NS(I) 1 READ(22,*)IDUMM,(U(K,I,J),K=1,3),(M(K,I,J),K=1,3) CLOSE(22) C SC INPUT FILE USES NANOMETERS, OMEGA in eV,debyes c C WORKING UNITS : ATOMIC UNITS !!!!! c c OUTPUT UNITS : cm-1 or nm, debyes c c Rewrite input and put everything into aus: c ============================================================ IOSC=0 DO 4 I=1,NG DO 4 J=2,NS(I) IOSC=IOSC+1 if(IOSC.gt.N0)call report('IOSC>N0') c c eV -> hartrees E(IOSC)=OMEGA(I,J)*8065.54093D0/cm DO 4 K=1,3 c c nm -> bohrs R(IOSC,K)=XT(K,I)*10.0D0/bohr c c debyes -> aus UA(IOSC,K)=U(K,I,J)/debye 4 UM(IOSC,K)=M(K,I,J)/debye c ============================================================ N=IOSC WRITE(6,*)' ',N,' OSCILLATORS ',NG,' GROUPS' inquire(file='DC.PAR',exist=lex) if(lex)then open(45,file='DC.PAR') read(45,*)EPS0,iwre close(45) else WRITE(6,*) WRITE(6,*)' Relative permitivity:' READ(5,*)EPS0 endif c c setup potential matrix DO 2101 I = 1,N 2101 V(I,I)=E(I) c DO 210 I = 1,N xi=R(I,1) yi=R(I,2) zi=R(I,3) uix=UA(I,1) uiy=UA(I,2) uiz=UA(I,3) DO 210 J = I+1,N xj=R(J,1) yj=R(J,2) zj=R(J,3) xji=xj-xi yji=yj-yi zji=zj-zi ujx=UA(J,1) ujy=UA(J,2) ujz=UA(J,3) D2=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2 rij=dsqrt(D2) D5=D2*D2*rij c c allow to interact only oscillators in different group (rij>0): VIJ=0.0d0 IF(rij.GT.1.0d-4)THEN uij=uix*ujx+uiy*ujy+uiz*ujz uir=uix*xji+uiy*yji+uiz*zji ujr=ujx*xji+ujy*yji+ujz*zji EPSR=EPS0 IF (EPS0.LT.0.0D0) EPSR=-EPS0*rij VIJ=(uij*D2-3.0d0*uir*ujr)/D5/EPSR ENDIF V(I,J)=VIJ 210 V(J,I)=VIJ vimax=abs(v(1,1)) imax=1 vijmax=abs(v(1,2)) ijimax=1 ijjmax=2 do i=1,N if(vimax.lt.abs(v(i,i)))then vimax=abs(v(i,i)) imax=i endif do j=i+1,N if(vijmax.lt.abs(v(i,j)))then vijmax=abs(v(i,j)) ijimax=i ijjmax=j endif enddo enddo write(6,659)imax,vimax,ijimax,ijjmax,vijmax 659 format( 'Maximal diagonal element ',i5,5x, 1f12.6,/,'Maximal off-diagonal element ',i5,i5,f12.6) if(iwre.ne.0)then open(56,file='V.TXT') do 101 I=1,N write(56,609)I 609 format(' Vector ',i5,':') 101 write(56,619)(V(I,J)*cm,J=1,I) 619 format(10f8.0) close(56) endif 212 WRITE(6,*)' I,J FOR THE H(I,J) {[0,0] TO CONTINUE, [-100,0]', 1',[-200,0] OPTIONS}:' if(lex)then I=0 J=0 else READ(5,*)I,J endif IF(I.EQ.-100)THEN WRITE(6,*)' 1-2,3-4,5-6 ETC. INTERACTIONS CANCELED' DO 67 I=1,N-1,2 V(I,I+1)=0.0d0 67 V(I+1,I)=0.0d0 GOTO 212 ENDIF IF(I.EQ.-200)THEN WRITE(6,*)' ODD-EVEN AND EVEN-ODD INTERACTIONS CANCELED' DO 68 I=1,N DO 68 J=1,I IF(MOD(I-J,2).NE.0)THEN V(I,J)=0.0d0 V(J,I)=0.0d0 ENDIF 68 CONTINUE GOTO 212 ENDIF IF(I*J.EQ.0)GOTO 211 WRITE(6,666)V(IABS(I),IABS(J)),' au' WRITE(6,666)V(IABS(I),IABS(J))*cm,' cm-1' WRITE(6,666)1.0d7/(V(IABS(I),IABS(J))*cm),' nm' 666 format(f12.6,a5) IF (I.LT.0) THEN I=IABS(I) J=IABS(J) WRITE(6,*)' NEW V(I,J):' READ(5,*)V(I,J) V(J,I)=V(I,J) ENDIF GOTO 212 211 WRITE(6,*)' DIAGONALIZING THE HAMILTONIAN' IERR=0 ICODE=2 CALL TRED12(N,V,COEF,EVAL,ICODE,IERR) IF (IERR.NE.0)call report('DIAGONALIZATION ERROR') WRITE(6,*)' DIAGONALIZATION SUCCESSFUL' c if(lex)then OK='N' else WRITE(6,*)' NORMALIZE THE RESULTS (Y/N/S) ? ' READ(5,'(A)')OK endif NN=1 IF (OK.EQ.'Y'.OR.OK.EQ.'y')NN=NG IF (OK.EQ.'S'.OR.OK.EQ.'s')THEN WRITE(6,*)' NN: ' READ(5,*)NN ENDIF Z=1.0d0/dble(NN) OPEN(23,FILE='DC.TAB') WRITE(23,2323) 2323 FORMAT( 1 ' MODE D RTOT RINT', 2/,' E DEBYE^2 debye^2=1e-36 cgs (esu**2 cm**2)', 3/,80(1H-)) DTOT=0.0d0 RTOT=0.0d0 DO 8 I=1,N c do 81 ix=1,3 AMMI(I,ix)=0.0d0 81 DM(I,ix)=0.0d0 ax=0.0d0 ay=0.0d0 az=0.0d0 DO 370 j = 1,N cij=COEF(j,I) ax=ax+cij*0.5d0*E(j)*(R(j,2)*ua(j,3)-R(j,3)*ua(j,2)) ay=ay+cij*0.5d0*E(j)*(R(j,3)*ua(j,1)-R(j,1)*ua(j,3)) az=az+cij*0.5d0*E(j)*(R(j,1)*ua(j,2)-R(j,2)*ua(j,1)) DO 370 ix=1,3 AMMI(I,ix)=AMMI(I,ix)+cij*UM(j,ix) 370 DM(I,ix)=DM(I,ix)+cij*UA(j,ix) AMM(I,1)=ax AMM(I,2)=ay AMM(I,3)=az D =DM(I,1)* DM(I,1)+DM(I,2)* DM(I,2)+DM(I,3)* DM(I,3) RI=DM(I,1)*AMMI(I,1)+DM(I,2)*AMMI(I,2)+DM(I,3)*AMMI(I,3) RS=DM(I,1)* AMM(I,1)+DM(I,2)* AMM(I,2)+DM(I,3)* AMM(I,3) c normally output in cm-1: c trasform to nanometers when configuration file present: EI=EVAL(I)*cm if(lex)EI=1.0d7/EI c D in debyes^2 D=D*debye**2 c c R in debye^2=1e-36 cgs - this is magic con=pi*bohr*cm*debye**2*2.0d0*1.0d-8 RS=RS*con RI=RI*con DTOT=DTOT+D*Z RTOT=RTOT+(RS+RI)*Z c 8 WRITE(23,2324)I,EI,D*Z,Z*(RS+RI),Z*RI 2324 FORMAT(I4,F10.3,4G16.6) WRITE(23,2325) 2325 FORMAT(80(1H-)) WRITE (23,450) DTOT,RTOT,NN 450 FORMAT ('DSUM=',G16.6,' ','RSUM=',G16.6,/, 1'SPECTRUM OF ALL SYSTEM DEVIDED BY ',I4) CLOSE(23) WRITE(6,*)' DC.TAB written' if(iwre.eq.1)CALL WREI(N0,N,COEF,DM,AMM,AMMI,EVAL,RRV,lex) WRITE(6,*)' >>> PROGRAM TERMINATED <<<' STOP END SUBROUTINE WREI(NM,N,COEF,DM,AMM,AMMI,EVAL,R,lex) IMPLICIT REAL*8(A-H,O-Z) DIMENSION COEF(NM,NM),DM(NM,3),EVAL(NM),AMM(NM,3),AMMI(NM,3), 1R(3,NM) REAL*8 ME(3),MM(3),AI(3,3),TI(3) logical*4 lex WRITE(6,*)' EIGEN.VEC ... FILE WITH THE EIGENVECTORS' OPEN(26,FILE='EIGEN.VEC') WRITE(26,*)' THE NORMALIZED EIGENVECTORS (BASE MODE):' NQ=N M=8 MI=-7 207 MI=MI+8 MJ=M-7 IF(M.GT.NQ)M=NQ WRITE(26,28)(I,I=MJ,M) 28 FORMAT(8X,8(I7,2X)) DO 206 I=1,NQ 206 WRITE(26,27)I,(COEF(J,I),J=MJ,M) 27 FORMAT(I5,3X,8F9.5) IF (M.LT.NQ)THEN M=M+8 GOTO 207 ENDIF CLOSE(26) write(6,*)'EIGEN.VEC written' DO 281 I=1,N DO 281 J=1,3 281 AMM(I,J)=AMM(I,J)+AMMI(I,J) OPEN(3,FILE='FILE.X',STATUS='OLD') READ(3,*) READ(3,*)NAT write(6,*)NAT if(NAT.gt.NM)call report('too many atoms') do 1 IA=1,NAT 1 READ(3,*)IDUMM,(R(ix,IA),ix=1,3) CLOSE(3) write(6,*)'Coordinates read from FILE.X' CALL INERTIA(R,AI,TI,NAT,NM) OPEN(26,FILE='DIP.MOM') WRITE(26,*)' PROJECTIONS OF TRANSITION DIPOLES:' WRITE(26,2627)(TI(I),I=1,3) 2627 FORMAT(3F15.5,/, 1'MODE DX DY DZ [%] / MX MY MZ [%]') DO 2280 I = 1,N EI=EVAL(I)*219470.0d0 if(lex)EI=1.0d7/EI DO 39 II=1,3 ME(II)=0.0d0 MM(II)=0.0d0 DO 39 J=1,3 MM(II)=MM(II)+AI(J,II)*AMM(I,J) 39 ME(II)=ME(II)+AI(J,II)*DM (I,J) ANORM=ME(1)**2+ME(2)**2+ME(3)**2 BNORM=MM(1)**2+MM(2)**2+MM(3)**2 DO 404 IX=1,3 IF(BNORM.GT.1.0D-38)MM(IX)=MM(IX)**2/BNORM*100.0d0 404 IF(ANORM.GT.1.0D-38)ME(IX)=ME(IX)**2/ANORM*100.0d0 2280 WRITE (26,2626)I,EI,(ME(IX),IX=1,3),(MM(IX),IX=1,3) 2626 FORMAT(I4,F8.1,3F6.1,' /',3F6.1) CLOSE(26) RETURN END SUBROUTINE INERTIA(R,A,TI,N,NMAX) IMPLICIT REAL*8(A-H,O-Z) DIMENSION R(3,NMAX),A(3,3),T(3,3),TI(3),CM(3) DO 2 I=1,3 CM(I)=0.0d0 DO 3 J=1,N 3 CM(I)=CM(I)+R(I,J) 2 CM(I)=CM(I)/N DO 1 I=1,3 DO 1 J=1,3 T(I,J)=0.0d0 DO 1 IS=1,N 1 T(I,J)=T(I,J)+(R(I,IS)-CM(I))*(R(J,IS)-CM(J)) CALL TRED(3,3,T,A,TI,2,IERR) RETURN C DEL(I,J)=SUM(OVER JS) OF A(I,JS)*A(J,JS) C DIAGONAL TD(I,J) = A(IS,I)*T(IS,JS)*A(JS,j) END SUBROUTINE TRED(NMAX,N,A,Z,D,IEIG,IERR) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(NMAX,NMAX),Z(NMAX,NMAX),D(NMAX) PARAMETER (N0=2000) COMMON/DIA/E(N0) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) JP1 = J + 1 IF (L .LT. JP1) GO TO 220 DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE HH = F / (H + H) DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) 290 D(I) = H 300 CONTINUE 320 D(1) = 0.0D0 E(1) = 0.0D0 IF (ICODE .NE. 2) GO TO 600 DO 500 I = 1, N L = I - 1 IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380 DO 360 J = 1, L G = 0.0D0 DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE 500 CONTINUE 620 CALL TQL (NMAX,N,Z,D,IERR,ICODE) RETURN 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END SUBROUTINE TQL(NMAX,N,Z,D,IERR,ICODE) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Z(NMAX,NMAX),D(NMAX) REAL*8 MACHEP PARAMETER (N0=2000) COMMON/DIA/E(N0) EPS = 1.0D0 10 EPS = 0.50D0*EPS TOL1 = EPS + 1.0D0 IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10 IF(TOL1-EPS.EQ.1.0D0) EPS = EPS + EPS MACHEP = EPS IERR = 0 IF (N .EQ. 1) GO TO 1001 DO 100 I = 2, N 100 E(I-1) = E(I) F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 110 CONTINUE 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 L1 = L + 1 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) R = DSQRT(P*P+1.0D0) D(L) = E(L) / (P + DSIGN(R,P)) H = G - D(L) DO 140 I = L1, N 140 D(I) = D(I) - H F = F + H P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L DO 200 II = 1, MML I = M - II G = C * E(I) H = C * P IF (DABS(P) .LT. DABS(E(I))) GO TO 150 C = E(I) / P R = DSQRT(C*C+1.0D0) E(I+1) = S * P * R S = C / R C = 1.0D0 / R GO TO 160 150 C = P / E(I) R = DSQRT(C*C+1.0D0) E(I+1) = S * E(I) * R S = 1.0D0 / R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) IF (ICODE .NE. 2) GO TO 200 DO 180 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 180 CONTINUE 200 CONTINUE E(L) = S * P D(L) = C * P IF (DABS(E(L)) .GT. B) GO TO 130 220 D(L) = D(L) + F 240 CONTINUE GO TO 1001 1000 IERR = L 1001 RETURN END subroutine report(s) character*(*) s write(6,60) write(6,*)s write(6,60) 60 format(60(1H*)) stop end SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (N0=2000) DIMENSION A(N0,N),Z(N0,N),D(N0),E(N0) COMMON/DIAG/E IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N C DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE C 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 C :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- :::::::::: DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 C :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) :::::::::: DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) C IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE C G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 C DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 C :::::::::: FORM ELEMENT OF A*U :::::::::: DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C :::::::::: FORM ELEMENT OF P :::::::::: 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C :::::::::: FORM REDUCED A :::::::::: DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0D0 E(1) = 0.0D0 C :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES :::::::::: IF (ICODE .NE. 2) GO TO 600 DO 500 I = 1, N L = I - 1 IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380 C DO 360 J = 1, L G = 0.0D0 C DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 C DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE C 500 CONTINUE C 620 CALL TQL12 (N,Z,D,IERR,ICODE) RETURN C :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY :::::::::: 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END C SUBROUTINE TQL12(N,Z,D,IERR,ICODE) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (N0=2000) DIMENSION Z(N0,N),D(N0),E(N0) REAL*8 MACHEP COMMON/DIAG/E EPS = 1.0D0 10 EPS = 0.50D0*EPS TOL1 = EPS + 1.0D0 IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10 IF(TOL1-EPS.EQ.1.0D0) EPS = EPS + EPS MACHEP = EPS C C MACHEP=16.0D0**(-13) C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 C DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT :::::::::: DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 C :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP :::::::::: 110 CONTINUE C 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C :::::::::: FORM SHIFT :::::::::: L1 = L + 1 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) R = DSQRT(P*P+1.0D0) D(L) = E(L) / (P + DSIGN(R,P)) H = G - D(L) C DO 140 I = L1, N 140 D(I) = D(I) - H C F = F + H C :::::::::: QL TRANSFORMATION :::::::::: P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L C :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- :::::::::: DO 200 II = 1, MML I = M - II G = C * E(I) H = C * P IF (DABS(P) .LT. DABS(E(I))) GO TO 150 C = E(I) / P R = DSQRT(C*C+1.0D0) E(I+1) = S * P * R S = C / R C = 1.0D0 / R GO TO 160 150 C = P / E(I) R = DSQRT(C*C+1.0D0) E(I+1) = S * E(I) * R S = 1.0D0 / R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C :::::::::: FORM VECTOR :::::::::: IF (ICODE .NE. 2) GO TO 200 DO 180 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 180 CONTINUE C 200 CONTINUE C E(L) = S * P D(L) = C * P IF (DABS(E(L)) .GT. B) GO TO 130 220 D(L) = D(L) + F 240 CONTINUE C :::::::::: ORDER EIGENVALUES AND EIGENVECTORS :::::::::: DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .LE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C IF (ICODE .NE. 2) GO TO 300 DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C :::::::::: SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS :::::::::: 1000 IERR = L 1001 RETURN END