C----------------------------------------------------------------------- C EMPIRICAL PSEUDOPOTENTIAL BANDSTRUCTURE CALCULATION C----------------------------------------------------------------------- C C----------------------------------------------------------------------- C C Written by: Umberto Ravaioli C C 3255 Beckman Institute C University of Illinois C at Urbana-Champaign C 405 N. Mathews Avenue C Urbana, IL 61801 C USA C C Tel: (217) 244-5765 C FAX: (217) 244-4333 C e-mail: ravaioli@ceg.uiuc.edu C C----------------------------------------------------------------------- C First version: July 26, 1994 C Last revision: August 17, 1994 C----------------------------------------------------------------------- C----------------------------------------------------------------------- C C This program calculates the bandstructure of face centered cubic C semiconductors using the local empirical pseudopotential method C as described in the references: C C [1] M.L. Cohen and T.K. Bergstresser, "Band Structures and C Pseudopotential Form Factors for Fourteen Semiconductors of the C Diamond and Zinc-blende Structures", Physical Review, vol. 141, C n. 2, pp. 789-796, 1966. C C [2] J.R. Chelikowsky and M.L. Cohen, "Nonlocal Pseudopotential C calculations for the electronic structure of eleven diamond and C zinc-blende semiconductors", Physical Review B, vol. 14, n. 5, C pp. 556-582, 1976. C C [3] M.L. Cohen and J.R. Chelikowsky, "Electronic Structure and Optical C Properties of Semiconductors", 2nd edition, Springer Verlag, 1989. C C This version uses the "local" pseudopotential approach and is C calibrated for silicon. C----------------------------------------------------------------------- C The present version generates a full matrix for the solution of C the eigenvalue problem. While the storage is not optimal, it is C convenient for linking commonly available solution packages. C For extensive applications, the use of specialized eigenvalue C solvers provides faster applications. A more efficient version C is under development. C----------------------------------------------------------------------- C This program is written in standard FORTRAN 77 and was developed C on a Hewlett-Packard 735 workstation. It should be portable to C most of the existing computers and workstations. C----------------------------------------------------------------------- C Type convention: Upper case is used for "permanent" FORTRAN C instructions. Lower case is used for "temporary" FORTRAN C instructions, comments, debugging additions. C----------------------------------------------------------------------- C This program is distributed by the National Center for C Computational Electronics (NCCE). Please acknowledge the C author and NCCE if the code is used for publications, reports, C or presentations. C----------------------------------------------------------------------- PROGRAM BANDS1 C----------------------------------------------------------------------- C DOUBLE PRECISION C----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) C----------------------------------------------------------------------- C PARAMETER LIST C----------------------------------------------------------------------- PARAMETER (IGSITES = 5000) PARAMETER (IGSIZE = 500) C----------------------------------------------------------------------- C COMMON BLOCK C----------------------------------------------------------------------- COMMON /MATHEM/ PI COMMON /PHYSIC/ BOHR,RYDBE COMMON /LATTIC/ CONSTA,NBASIS,NSITES COMMON /PSEUDO/ PSEUD1,PSEUD2,PSEUD3 COMMON /POTENT/ PSPOT1(IGSITES),POTVAL(IGSITES),IFLAG(IGSITES) COMMON /BSCAL1/ EIGMAT(IGSIZE,IGSIZE),IBAND(IGSIZE) COMMON /BSCAL2/ ENERGY(IGSIZE),E(IGSIZE),E2(IGSIZE) COMMON /CUTOFF/ EMAX,CUTOF1,CUTOF2,CUTOF3 COMMON /CRYSTA/ BASIS(3,3),RBASIS(3,3),ATOM1(3),ATOM2(3) COMMON /RECIP1/ GVECX(IGSITES),GVECY(IGSITES),GVECZ(IGSITES) COMMON /RECIP2/ GVMAG(IGSITES),IARRAY(IGSIZE,IGSIZE) C----------------------------------------------------------------------- C OPEN OUTPUT FILE C----------------------------------------------------------------------- OPEN(UNIT=20,file='CBAND.OUT') OPEN(UNIT=30,file='VBAND.OUT') C----------------------------------------------------------------------- C CONSTANTS C----------------------------------------------------------------------- C PI = 3.14 ... C BOHR = Bohr radius (in Angstroms) C RYDBE = 1 Rydberg (in electron Volts) C----------------------------------------------------------------------- PI = ACOS(-1.D0) BOHR = 0.52918D0 RYDBE = 13.606D0 C----------------------------------------------------------------------- C INPUTS C----------------------------------------------------------------------- C Inputs (Silicon considered here): C C CONSTA = Lattice constant (in Angstroms) C NBASIS = Atoms in the basis (here diamond lattice) C NSITES = Number of sites included in the calculation C PSEUDi = Pseudopotential Fourier components C EMAX = Maximum energy (cutoff) C NEIGEN = Number of eigenvalue branches considered C NSTEP = Number of step in momentum C STEPK = Mesh size in momentum space C VALREF = Energy reference on top of Si valence band (eV). C This value was computed separately with this code, and C is used to scale the energy eigenvalues, so that the C top of the valence band has zero energy. C CONREF = Additional energy amount to put zero refernce at the C bottom of the conduction band. The value was also C calculated using this code, looking at the energy C minimum near the X-point. This may be slightly C approximate. The minimum for the conduction band C is taken at (0.85,0,0) for the normalized momentum C----------------------------------------------------------------------- CONSTA = 5.43097D0 NBASIS = 2 NSITES = 16 PSEUD1 = -0.224D0 PSEUD2 = 0.055D0 PSEUD3 = 0.072D0 EMAX = 10.D0 NEIGEN = 10 NSTEP = 10 STEPK = 0.005 VALREF = 10.23761802283895D0 CONREF = 1.03919625763984D0 C----------------------------------------------------------------------- C NORMALIZE LATTICE CONSTANT C----------------------------------------------------------------------- CONSTA = CONSTA/BOHR HALFA = CONSTA/2.0D0 C----------------------------------------------------------------------- C DIRECT LATTICE C----------------------------------------------------------------------- C BASIS(i,j) = basis vectors of primitive cell C (face centered cubic lattice) C i = 1,2,3 vector index C j = 1,2,3 coordinate index (x,y,z) C C VOLDIR = Volume of primitive cell C C ATOM1(i) = Coordinates of first atom in the basis C ATOM2(i) = Coordinates of second atom in the basis C (The reference lattice site is at (0,0,0)) C C Here we assume that the origin of real space is taken half way C the two fcc lattices of the diamond lattice, so that the atoms C of the real space basis are identified by vectors (1/8,1/8,1/8) C and (-1/8,-1/8,-1/8). C C----------------------------------------------------------------------- BASIS(1,1) = - HALFA BASIS(1,2) = 0.0D0 BASIS(1,3) = HALFA BASIS(2,1) = 0.0D0 BASIS(2,2) = HALFA BASIS(2,3) = HALFA BASIS(3,1) = - HALFA BASIS(3,2) = HALFA BASIS(3,3) = 0.0D0 VOLDIR = 2.0D0 * HALFA*HALFA*HALFA CONST8 = CONSTA/8.0D0 ATOM1(1) = CONST8 ATOM1(2) = CONST8 ATOM1(3) = CONST8 ATOM2(1) = - CONST8 ATOM2(2) = - CONST8 ATOM2(3) = - CONST8 C----------------------------------------------------------------------- C RECIPROCAL LATTICE C----------------------------------------------------------------------- C CONSTR = Lattice constant of reciprocal lattice C FACTOR = Lattice constant / volume of primitive cell C C RBASIS(i,j) = basis vectors of reciprocal lattice primitive cell C (face centered cubic direct lattice) C i = 1,2,3 vector index C j = 1,2,3 coordinate index (x,y,z) C C IMAX = Maximum index for loops C----------------------------------------------------------------------- CONSTR = 2.0D0*PI/CONSTA FACTOR = CONSTA/VOLDIR CUTOF1 = EMAX/CONSTR/CONSTR CUTOF2 = (DSQRT(5.0D0)/2.0D0 + DSQRT(CUTOF1))**2 CUTOF3 = CUTOF2*4.0D0 RBASIS(1,1) = - FACTOR * BASIS(2,3) * BASIS(3,2) RBASIS(1,2) = FACTOR * BASIS(2,3) * BASIS(3,1) RBASIS(1,3) = - FACTOR * BASIS(2,2) * BASIS(3,1) RBASIS(2,1) = FACTOR * BASIS(3,2) * BASIS(1,3) RBASIS(2,2) = - FACTOR * BASIS(3,1) * BASIS(1,3) RBASIS(2,3) = - FACTOR * BASIS(3,2) * BASIS(1,1) RBASIS(3,1) = - FACTOR * BASIS(1,3) * BASIS(2,2) RBASIS(3,2) = - FACTOR * BASIS(1,1) * BASIS(2,3) RBASIS(3,3) = FACTOR * BASIS(1,1) * BASIS(2,2) IMAX = 2*NSITES -1 INDEX1 = 0 INDEX2 = 0 C----------------------------------------------------------------------- C Construct a "database" of reciprocal lattice vectors, considering C a region around (0,0,0) in momentum space, which extends along C the basis vector directions from -NSITES to NSITES . The C vectors contained in the array have magnitude smaller than C a specified safety margin (CUTOF3). The total number of vectors C is INDEX1. C C GX,GY,GZ = Components of the reciprocal vectors C GMAG = Magnitude of the reciprocal vectors C----------------------------------------------------------------------- DO 10 I=1,IMAX DO 10 J=1,IMAX DO 10 K=1,IMAX COEF1 = DBLE(I-NSITES) COEF2 = DBLE(J-NSITES) COEF3 = DBLE(K-NSITES) GX = COEF1*RBASIS(1,1)+COEF2*RBASIS(2,1)+ # COEF3*RBASIS(3,1) GY = COEF1*RBASIS(1,2)+COEF2*RBASIS(2,2)+ # COEF3*RBASIS(3,2) GZ = COEF1*RBASIS(1,3)+COEF2*RBASIS(2,3)+ # COEF3*RBASIS(3,3) GMAG = GX*GX + GY*GY + GZ*GZ IF (GMAG.LE.CUTOF3) THEN INDEX1 = INDEX1 + 1 GVECX(INDEX1) = GX GVECY(INDEX1) = GY GVECZ(INDEX1) = GZ GVMAG(INDEX1) = GMAG IF (GMAG.LE.CUTOF2) INDEX2 = INDEX2 + 1 ENDIF 10 CONTINUE C----------------------------------------------------------------------- C Reciprocal lattice vectors are ordered according to magnitude C----------------------------------------------------------------------- DO 20 I=1,INDEX1 DO 20 J=I,INDEX1 IF (GVMAG(J).LT.GVMAG(I)) THEN TEMPX = GVECX(I) TEMPY = GVECY(I) TEMPZ = GVECZ(I) TEMPM = GVMAG(I) GVECX(I) = GVECX(J) GVECY(I) = GVECY(J) GVECZ(I) = GVECZ(J) GVMAG(I) = GVMAG(J) GVECX(J) = TEMPX GVECY(J) = TEMPY GVECZ(J) = TEMPZ GVMAG(J) = TEMPM ENDIF 20 CONTINUE C----------------------------------------------------------------------- DO 30 I=2,INDEX2 DO 30 J=1,I-1 C --------------------------------------------------- IARRAY(I,J) = 0 TEMPX = GVECX(I) - GVECX(J) TEMPY = GVECY(I) - GVECY(J) TEMPZ = GVECZ(I) - GVECZ(J) TEMPM = TEMPX*TEMPX + TEMPY*TEMPY + TEMPZ*TEMPZ DO 40 K=2,INDEX1 IF ( ABS(TEMPM - GVMAG(K)).LT.(1.0D-3) ) THEN IF(( ABS(TEMPX - GVECX(K)).LT.(1.0D-3) ).AND. # ( ABS(TEMPY - GVECY(K)).LT.(1.0D-3) ).AND. # ( ABS(TEMPZ - GVECZ(K)).LT.(1.0D-3) )) THEN IARRAY(I,J) = K GO TO 30 ENDIF ENDIF 40 CONTINUE C --------------------------------------------------- 30 CONTINUE C----------------------------------------------------------------------- C In correspondence of G^2 = 3,8,11 the values for the C pseudopotential in reciprocal space are assigned. As indicated C in references [1,2,3], only these points are necessary for C covalent semiconductors to define the pseudopotential in C reciprocal space. C----------------------------------------------------------------------- DO 50 I=1,INDEX1 PSPOT1(I) = 0.0D0 IFLAG(I) = 0 50 CONTINUE DO 60 I=1,INDEX1 IF (GVMAG(I).GT.(2.99).AND.GVMAG(I).LT.(3.01)) THEN POTVAL(I) = PSEUD1 IFLAG(I) = 1 ELSE IF (GVMAG(I).GT.(7.99).AND.GVMAG(I).LT.(8.01)) THEN POTVAL(I) = PSEUD2 IFLAG(I) = 1 ELSE IF (GVMAG(I).GT.(10.99).AND.GVMAG(I).LT.(11.01)) THEN POTVAL(I) = PSEUD3 IFLAG(I) = 1 ENDIF ENDIF ENDIF 60 CONTINUE C----------------------------------------------------------------------- C The potentials for reciprocal lattice vectors are evaluated as a C product of the appropriate pseudopotential component and the C structure factor. For covalent semiconductors, the structure C factor is simply the cosine of the scalar product between the C reciprocal vector and the real space vector of the basis atom. C----------------------------------------------------------------------- DO 70 I=1,INDEX1 IF (IFLAG(I).EQ.1) THEN VALUE1 = (GVECX(I)*ATOM1(1) + GVECY(I)*ATOM1(2) + # GVECZ(I)*ATOM1(3))* CONSTR VALUE2 = (GVECX(I)*ATOM2(1) + GVECY(I)*ATOM2(2) + # GVECZ(I)*ATOM2(3))* CONSTR PSPOT1(I) = (DCOS(VALUE1) + DCOS(VALUE2)) * POTVAL(I)/2.0D0 ENDIF 70 CONTINUE C----------------------------------------------------------------------- C CALCULATION OF THE BANDSTRUCTURE C C The program is set up to calculate the bandstructure in a C rectangular region of momentum space, which contains the C irreducible wedge of the Brillouin zone, with some computational C redundancy. C To calculate specific regions or on specific orientations, simply C modify the loops DO 1000. C As set up below, BANDX goes from the Gamma to the X point. C C BANDX, BANDY and BANDZ are the components of the crystal momentum. C----------------------------------------------------------------------- BANDX = 0.0D0 BANDY = 0.0D0 BANDZ = 0.0D0 DO 1000 I1=1,NSTEP BANDX = DBLE(I1-1)*STEPK DO 1000 I2=1,NSTEP BANDY = DBLE(I2-1)*STEPK DO 1000 I3=1,NSTEP BANDZ = DBLE(I3-1)*STEPK C ------------------------------------------------------------------ C Generation of the matrix EIGMAT(I,J) for the eigenvalue problem C NOTE: The matrix is dense C ------------------------------------------------------------------ DO 80 I=1,INDEX2 DO 80 J=1,INDEX2 EIGMAT(I,J)=0.0D0 80 CONTINUE C ------------------------------------------------------------------ C Diagonal Elements C C Note: The size of the matrix is determined by the number of C reciprocal lattice sites within the specified cutoff and C is stored in the integer variable NBAND. C ------------------------------------------------------------------ NBAND=0 DO 100 I=1,INDEX2 TEMPX = BANDX + GVECX(I) TEMPY = BANDY + GVECY(I) TEMPZ = BANDZ + GVECZ(I) TEMPM = TEMPX*TEMPX + TEMPY*TEMPY + TEMPZ*TEMPZ IF (TEMPM.LE.CUTOF1) THEN NBAND = NBAND+1 EIGMAT(NBAND,NBAND) = TEMPM*CONSTR*CONSTR IBAND(NBAND) = I ENDIF 100 CONTINUE C ------------------------------------------------------------------ C Off-diagonal Elements C Note: For silicon the matrix is real and symmetric C ------------------------------------------------------------------ DO 110 I=2,NBAND DO 110 J=1,I-1 IF (IARRAY(IBAND(I),IBAND(J)).GT.0) THEN EIGMAT(I,J) = PSPOT1(IARRAY(IBAND(I),IBAND(J))) EIGMAT(J,I) = EIGMAT(I,J) ENDIF 110 CONTINUE C----------------------------------------------------------------------- C SOLVE THE EIGENVALUE PROBLEM C----------------------------------------------------------------------- C The following calling sequences are suggested for use with EISPACK C routines (FORTRAN IV). Details on EISPACK can be found on: C C [4] B.S. Garbow, "EISPACK-A Package of Matrix Eigensystem Routines", C Computer Physics Communications, North-Holland, Amsterdam, vol. 7, C pp. 179-184, 1974. C C [5] B.T. Smith, J.M. Boyle, J.J. Dongarra, B.S. Garbow, Y. Ikebe, C V.C. Klema, and C.B. Moler, "Matrix Eignsystem Routines - EISPACK C Guide", Lecture Notes in Computer Science, vol. 6, Second Edition, C Springer Verlag 1976. C C [6] B.S. Garbow, J.M. Boyle, J.J. Dongarra, C.B. Moler, "Matrix C Eigensystem Routines - EISPACK Guide Extension", Lecture Notes in C Computer Science, vol. 51, Springer Verlag 1977. C C----------------------------------------------------------------------- C If only eigenvalues are necessary, one can use the EISPACK C routines TRED1, TQLRAT. C C Subroutine TRED1 reduces a Real Symmetric Matrix to a Symmetric C Tridiagonal Matrix using orthogonal transformation. A discussion C of the method and a listing of the subroutine can be found on C page 302 of reference [6]. C C Subroutine TQLRAT determines the eigenvalues of a Symmetric C Tridiagonal Matrix using a rational variant of the QL method. C A discussion of the method and a listing of the subroutine can C be found on page 286 of reference [6]. C----------------------------------------------------------------------- C If eigenvalues and eigenvectors are necessary, one can use the C EISPACK routines TRED2 and TQL2. The matrix is not restored C at the end of the procedure. C C Subroutine TRED2 reduces a Real Symmetric Matrix to a Symmetric C Tridiagonal Matrix, and is similar to TRED1. A listing can be C found on page 305 of reference [6]. C C Subroutine TQL2 determines eigenvalues and eigenvectors of a C Tridiagonal Symmetric Matrix using the standard QL method. A C listing can be found on page 291 of reference [6]. C----------------------------------------------------------------------- C The present release uses TRED1 and TQLRAT C----------------------------------------------------------------------- C----------------------------------------------------------------------- C CALL EISPACK ROUTINES C----------------------------------------------------------------------- CALL TRED1(IGSIZE,NBAND,EIGMAT,ENERGY,E,E2) CALL TQLRAT(NBAND,ENERGY,E2,IERR) C----------------------------------------------------------------------- DO 120 I=1,NEIGEN ENERGY(I) = ENERGY(I)*RYDBE - VALREF - CONREF WRITE(20,2000) I1,I2,I3,ENERGY(5),ENERGY(6),ENERGY(7),ENERGY(8) WRITE(30,2000) I1,I2,I3,ENERGY(4),ENERGY(3),ENERGY(2),ENERGY(1) 2000 FORMAT(3(I3,1x),4(F8.6,1x)) 120 CONTINUE C----------------------------------------------------------------------- 1000 CONTINUE C----------------------------------------------------------------------- CLOSE(10) CLOSE(11) C----------------------------------------------------------------------- STOP END C----------------------------------------------------------------------- c----------------------------------------------------------------- SUBROUTINE TRED1(NM,N,A,D,E,E2) C INTEGER I,J,K,L,N,II,NM,JP1 DOUBLE PRECISION A(NM,N),D(N),E(N),E2(N) DOUBLE PRECISION F,G,H,SCALE C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED1, C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX C TO A SYMMETRIC TRIDIAGONAL MATRIX USING C ORTHOGONAL SIMILARITY TRANSFORMATIONS. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C A CONTAINS THE REAL SYMMETRIC INPUT MATRIX. ONLY THE C LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. C C ON OUTPUT C C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS- C FORMATIONS USED IN THE REDUCTION IN ITS STRICT LOWER C TRIANGLE. THE FULL UPPER TRIANGLE OF A IS UNALTERED. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C DO 100 I = 1, N D(I) = A(N,I) A(N,I) = A(I,I) 100 CONTINUE C .......... FOR I=N STEP -1 UNTIL 1 DO -- .......... DO 300 II = 1, N I = N + 1 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 1) GO TO 130 C .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... DO 120 K = 1, L 120 SCALE = SCALE + DABS(D(K)) C IF (SCALE .NE. 0.0D0) GO TO 140 C DO 125 J = 1, L D(J) = A(L,J) A(L,J) = A(I,J) A(I,J) = 0.0D0 125 CONTINUE C 130 E(I) = 0.0D0 E2(I) = 0.0D0 GO TO 300 C 140 DO 150 K = 1, L D(K) = D(K) / SCALE H = H + D(K) * D(K) 150 CONTINUE C E2(I) = SCALE * SCALE * H F = D(L) G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G D(L) = F - G IF (L .EQ. 1) GO TO 285 C .......... FORM A*U .......... DO 170 J = 1, L 170 E(J) = 0.0D0 C DO 240 J = 1, L F = D(J) G = E(J) + A(J,J) * F JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L G = G + A(K,J) * D(K) E(K) = E(K) + A(K,J) * F 200 CONTINUE C 220 E(J) = G 240 CONTINUE C .......... FORM P .......... F = 0.0D0 C DO 245 J = 1, L E(J) = E(J) / H F = F + E(J) * D(J) 245 CONTINUE C H = F / (H + H) C .......... FORM Q .......... DO 250 J = 1, L 250 E(J) = E(J) - H * D(J) C .......... FORM REDUCED A .......... DO 280 J = 1, L F = D(J) G = E(J) C DO 260 K = J, L 260 A(K,J) = A(K,J) - F * E(K) - G * D(K) C 280 CONTINUE C 285 DO 290 J = 1, L F = D(J) D(J) = A(L,J) A(L,J) = A(I,J) A(I,J) = F * SCALE 290 CONTINUE C 300 CONTINUE C RETURN END c--------------------------------------------------------------- c SUBROUTINE TQLRAT(N,D,E2,IERR) C INTEGER I,J,L,M,N,II,L1,MML,IERR DOUBLE PRECISION D(N),E2(N) DOUBLE PRECISION B,C,F,G,H,P,R,S,T,EPSLON,PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT, C ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH. C C THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC C TRIDIAGONAL MATRIX BY THE RATIONAL QL METHOD. C C ON INPUT C C N IS THE ORDER OF THE MATRIX. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. C C E2 CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF THE C INPUT MATRIX IN ITS LAST N-1 POSITIONS. E2(1) IS ARBITRARY. C C ON OUTPUT C C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND C ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE C THE SMALLEST EIGENVALUES. C C E2 HAS BEEN DESTROYED. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E2(I-1) = E2(I) C F = 0.0D0 T = 0.0D0 E2(N) = 0.0D0 C DO 290 L = 1, N J = 0 H = DABS(D(L)) + DSQRT(E2(L)) IF (T .GT. H) GO TO 105 T = H B = EPSLON(T) C = B * B C .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT .......... 105 DO 110 M = L, N IF (E2(M) .LE. C) GO TO 120 C .......... E2(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 210 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C .......... FORM SHIFT .......... L1 = L + 1 S = DSQRT(E2(L)) G = D(L) P = (D(L1) - G) / (2.0D0 * S) R = PYTHAG(P,1.0D0) D(L) = S / (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 .......... RATIONAL QL TRANSFORMATION .......... G = D(M) IF (G .EQ. 0.0D0) G = B H = G S = 0.0D0 MML = M - L C .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... DO 200 II = 1, MML I = M - II P = G * H R = P + E2(I) E2(I+1) = S * R S = E2(I) / R D(I+1) = H + S * (H + D(I)) G = D(I) - E2(I) / G IF (G .EQ. 0.0D0) G = B H = G * P / R 200 CONTINUE C E2(L) = S * G D(L) = H C .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST .......... IF (H .EQ. 0.0D0) GO TO 210 IF (DABS(E2(L)) .LE. DABS(C/H)) GO TO 210 E2(L) = H * E2(L) IF (E2(L) .NE. 0.0D0) GO TO 130 210 P = D(L) + F C .......... ORDER EIGENVALUES .......... IF (L .EQ. 1) GO TO 250 C .......... FOR I=L STEP -1 UNTIL 2 DO -- .......... DO 230 II = 2, L I = L + 2 - II IF (P .GE. D(I-1)) GO TO 270 D(I) = D(I-1) 230 CONTINUE C 250 I = 1 270 D(I) = P 290 CONTINUE C GO TO 1001 C .......... SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS .......... 1000 IERR = L 1001 RETURN END c--------------------------------------------------------------- c DOUBLE PRECISION FUNCTION EPSLON (X) DOUBLE PRECISION X C C ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X. C DOUBLE PRECISION A,B,C,EPS C C THIS PROGRAM SHOULD FUNCTION PROPERLY ON ALL SYSTEMS C SATISFYING THE FOLLOWING TWO ASSUMPTIONS, C 1. THE BASE USED IN REPRESENTING FLOATING POINT C NUMBERS IS NOT A POWER OF THREE. C 2. THE QUANTITY A IN STATEMENT 10 IS REPRESENTED TO C THE ACCURACY USED IN FLOATING POINT VARIABLES C THAT ARE STORED IN MEMORY. C THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO C FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING C ASSUMPTION 2. C UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT, C A IS NOT EXACTLY EQUAL TO FOUR-THIRDS, C B HAS A ZERO FOR ITS LAST BIT OR DIGIT, C C IS NOT EXACTLY EQUAL TO ONE, C EPS MEASURES THE SEPARATION OF 1.0 FROM C THE NEXT LARGER FLOATING POINT NUMBER. C THE DEVELOPERS OF EISPACK WOULD APPRECIATE BEING INFORMED C ABOUT ANY SYSTEMS WHERE THESE ASSUMPTIONS DO NOT HOLD. C C THIS VERSION DATED 4/6/83. C A = 4.0D0/3.0D0 10 B = A - 1.0D0 C = B + B + B EPS = DABS(C-1.0D0) IF (EPS .EQ. 0.0D0) GO TO 10 EPSLON = EPS*DABS(X) RETURN END DOUBLE PRECISION FUNCTION PYTHAG(A,B) DOUBLE PRECISION A,B C C FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW C DOUBLE PRECISION P,R,S,T,U P = DMAX1(DABS(A),DABS(B)) IF (P .EQ. 0.0D0) GO TO 20 R = (DMIN1(DABS(A),DABS(B))/P)**2 10 CONTINUE T = 4.0D0 + R IF (T .EQ. 4.0D0) GO TO 20 S = R/T U = 1.0D0 + 2.0D0*S P = U*P R = (S/U)**2 * R GO TO 10 20 PYTHAG = P RETURN END