subroutine DHERQL(AR, AI, LDA, N, EVAL, VR, VI, WORK, IERR) c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. c ALL RIGHTS RESERVED. c Based on Government Sponsored Research NAS7-03001. C>> 1994-10-20 DHERQL Krogh Changes to use M77CON C>> 1994-04-20 DHERQL CLL Edited comment to make DP & SP files similar. c>> 1992-04-23 CLL Declaring all variables. c>> 1992-04-08 DHERQL Removed unused label 220. c>> 1991-10-23 DHERQL Krogh Initial version, converted from EISPACK. c c This subroutine reduces a complex hermitian matrix to a real c symmetric tridiagonal matrix using unitary similarity c transformations, computes the eigenvalues and eigenvectors of the c tridiagonal matrix, and then computes the eigenvectors of the c original matrix. c c This subroutine is a concatenation of a slight modification of the c EISPACK subroutine HTRIDI, a call to DIMQL which is a slight c modification of the EISPACK subroutine IMTQL2, and a slight c modification of the EISPACK subroutine HTRIBK. These subroutines c in turn are a translation of a complex analogue of the ALGOL c procedure TRED1, Num. Math. 11, 181-195(1968) by Martin, Reinsch, c and Wilkinson. Handbook for Auto. Comp., Vol.ii-Linear Algebra, c 212-226(1971), and the complex analogue of the ALGOL procedure c TRBAK1, Num. Math. 11, 181-195(1968) by Martin, Reinsch, and c Wilkinson. Handbook for Auto. Comp., vol.ii-Linear Algebra, c 212-226(1971). c EISPACK routines are dated August 1983. c c On input c c AR and AI contain the real and imaginary parts, c respectively, of the complex hermitian input matrix. c only the lower triangle of the matrix need be supplied. c c LDA 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 On output c c AR and AI contain information about the unitary transformations c used in the reduction in their full lower triangles. Their c strict upper triangles and the diagonal of AR are unaltered. c c EVAL contains the eigenvalues in ascending order. If an c error exit is made, the eigenvalues are correct but c unordered for indices 1,2,...,IERR-1, and no eigenvectors are c computed. c c VR and VI contain the the real and imaginary parts of the c eigenvectors with column k corresponding to the eigenvalue in c EVAL(k). Note that the last component of each returned vector c is real and that vector euclidean norms are preserved. c c WORK used for working storage. c c IERR is set to: c zero for normal return, c J if the J-th eigenvalue has not been determined after 30 c iterations. c c ------------------------------------------------------------------ c--D replaces "?": ?HERQL, ?IMQL c Both versions use IERM1 c ------------------------------------------------------------------ integer I, IERR, J,K,L,N,LDA,N1,N2 double precision AR(LDA,N),AI(LDA,N),EVAL(N),WORK(*) double precision VR(LDA,N),VI(LDA,N) double precision F,G,H,FI,GI,HH,SI,SCALE double precision S c ------------------------------------------------------------------ N1 = N N2 = N1 + N1 WORK(N2) = 1.0D0 WORK(N2 + N1) = 0.0D0 c do 100 I = 1, N1 EVAL(I) = AR(I,I) 100 continue if (N1 .le. 1) then if (N1 .eq. 1) then AR(1,1) = 1.0d0 AI(1,1) = 0.0d0 else call IERM1('DHERQL',-1, 0, 'Require N > 0.', 'N', N1, '.') IERR = -1 end if return end if do 300 I = N1, 1, -1 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 SCALE = SCALE + abs(AR(I,K)) + abs(AI(I,K)) 120 continue c if (SCALE .ne. 0.0D0) go to 140 WORK(N1 + L) = 1.0D0 WORK(N2 + L) = 0.0D0 130 WORK(I) = 0.0D0 go to 290 c 140 do 150 K = 1, L AR(I,K) = AR(I,K) / SCALE AI(I,K) = AI(I,K) / SCALE H = H + AR(I,K) * AR(I,K) + AI(I,K) * AI(I,K) 150 continue c G = sqrt(H) WORK(I) = SCALE * G F = abs(AR(I,L)) SI = abs(AI(I,L)) if (F .gt. SI) then F = F * sqrt(1.0d0 + (SI / F) **2) else if (SI .ne. 0.0d0) then F = SI * sqrt(1.0d0 + (F / SI) **2) else WORK(N1 + L) = -WORK(N1 + I) SI = WORK(N2 + I) AR(I,L) = G go to 170 end if c .......... Form next diagonal element of matrix t .......... WORK(N1+L) = (AI(I,L) * WORK(N2+I) - AR(I,L) * WORK(N1+I)) / F SI = (AR(I,L) * WORK(N2 + I) + AI(I,L) * WORK(N1 + I)) / F H = H + F * G G = 1.0D0 + G / F AR(I,L) = G * AR(I,L) AI(I,L) = G * AI(I,L) if (L .eq. 1) go to 270 F = 0.0D0 c 170 do 240 J = 1, L G = 0.0D0 GI = 0.0D0 c .......... form element of a*u .......... do 180 K = 1, J G = G + AR(J,K) * AR(I,K) + AI(J,K) * AI(I,K) GI = GI - AR(J,K) * AI(I,K) + AI(J,K) * AR(I,K) 180 continue c do 200 K = J+1, L G = G + AR(K,J) * AR(I,K) - AI(K,J) * AI(I,K) GI = GI - AR(K,J) * AI(I,K) - AI(K,J) * AR(I,K) 200 continue c .......... form element of p .......... WORK(J) = G / H WORK(N2 + J) = GI / H F = F + WORK(J) * AR(I,J) - WORK(N2 + J) * AI(I,J) 240 continue c HH = F / (H + H) c .......... form reduced a .......... do 260 J = 1, L F = AR(I,J) G = WORK(J) - HH * F WORK(J) = G FI = -AI(I,J) GI = WORK(N2 + J) - HH * FI WORK(N2 + J) = -GI c do 250 K = 1, J AR(J,K) = AR(J,K) - F * WORK(K) - G * AR(I,K) x + FI * WORK(N2 + K) + GI * AI(I,K) AI(J,K) = AI(J,K) - F * WORK(N2 + K) - G * AI(I,K) x - FI * WORK(K) - GI * AR(I,K) 250 continue 260 continue c 270 do 280 K = 1, L AR(I,K) = SCALE * AR(I,K) AI(I,K) = SCALE * AI(I,K) 280 continue c WORK(N2 + L) = -SI 290 HH = EVAL(I) EVAL(I) = AR(I,I) AR(I,I) = HH AI(I,I) = SCALE * sqrt(H) 300 continue c c ------------------------------------------------------------------ c do 420 I = 1, N1 do 410 J = 1, N1 VR(I,J) = 0.0d0 410 continue VR(I,I) = 1.0d0 420 continue call DIMQL(VR, LDA, N1, EVAL, WORK, IERR) if (IERR .ne. 0) return c c ------------------------------------------------------------------ c c .......... transform the eigenvectors of the real symmetric c tridiagonal matrix to those of the hermitian c tridiagonal matrix. .......... do 550 K = 1, N1 c do 540 J = 1, N1 VI(K,J) = -VR(K,J) * WORK(N2 + K) VR(K,J) = VR(K,J) * WORK(N1 + K) 540 continue 550 continue c c .......... Recover and apply the Householder matrices .......... do 640 I = 2, N1 L = I - 1 H = AI(I,I) if (H .ne. 0.0D0) then do 630 J = 1, N1 S = 0.0D0 SI = 0.0D0 do 610 K = 1, L S = S + AR(I,K) * VR(K,J) - AI(I,K) * VI(K,J) SI = SI + AR(I,K) * VI(K,J) + AI(I,K) * VR(K,J) 610 continue c .......... double divisions avoid possible underflow .......... S = (S / H) / H SI = (SI / H) / H c do 620 K = 1, L VR(K,J) = VR(K,J) - S * AR(I,K) - SI * AI(I,K) VI(K,J) = VI(K,J) - SI * AR(I,K) + S * AI(I,K) 620 continue 630 continue end if 640 continue c return end