1 SUBROUTINE CHIEV(A,LDA,N,E,V,LDV,WORK,JOB,INFO) 2C***BEGIN PROLOGUE CHIEV 3C***DATE WRITTEN 800808 (YYMMDD) 4C***REVISION DATE 820801 (YYMMDD) 5C***REVISION HISTORY (YYMMDD) 6C 000330 Modified array declarations. (JEC) 7C***CATEGORY NO. D4A3 8C***KEYWORDS COMPLEX,COMPLEX HERMITIAN,EIGENANALYSIS,EIGENVALUE, 9C EIGENVECTOR,MATRIX 10C***AUTHOR KAHANER, D. K., (NBS) 11C MOLER, C. B., (U. OF NEW MEXICO) 12C STEWART, G. W., (U. OF MARYLAND) 13C***PURPOSE To compute the eigenvalues and, optionally, the eigen- 14C vectors of a COMPLEX HERMITIAN matrix. 15C***DESCRIPTION 16C 17C David Kahaner, Cleve Moler, G. W. Stewart, 18C N.B.S. U.N.M. N.B.S./U.MD. 19C 20C Abstract 21C CHIEV computes the eigenvalues and, optionally, 22C the eigenvectors of a complex Hermitian matrix. 23C 24C Call Sequence Parameters- 25C (the values of parameters marked with * (star) will be changed 26C by CHIEV.) 27C 28C A* COMPLEX(LDA,N) 29C complex Hermitian input matrix. 30C Only the upper triangle of A need be 31C filled in. Elements on diagonal must be real. 32C 33C LDA INTEGER 34C set by the user to 35C the leading dimension of the complex array A. 36C 37C N INTEGER 38C set by the user to 39C the order of the matrices A and V, and 40C the number of elements in E. 41C 42C E* REAL(N) 43C on return from CHIEV E contains the eigenvalues of A. 44C See also INFO below. 45C 46C V* COMPLEX(LDV,N) 47C on return from CHIEV if the user has set JOB 48C = 0 V is not referenced. 49C = nonzero the N eigenvectors of A are stored in the 50C first N columns of V. See also INFO below. 51C 52C LDV INTEGER 53C set by the user to 54C the leading dimension of the array V if JOB is also 55C set nonzero. In that case N must be .LE. LDV. 56C If JOB is set to zero LDV is not referenced. 57C 58C WORK* REAL(4N) 59C temporary storage vector. Contents changed by CHIEV. 60C 61C JOB INTEGER 62C set by the user to 63C = 0 eigenvalues only to be calculated by CHIEV. 64C Neither V nor LDV are referenced. 65C = nonzero eigenvalues and vectors to be calculated. 66C In this case A and V must be distinct arrays 67C also if LDA .GT. LDV CHIEV changes all the 68C elements of A thru column N. If LDA < LDV 69C CHIEV changes all the elements of V through 70C column N. If LDA = LDV only A(I,J) and V(I, 71C J) for I,J = 1,...,N are changed by CHIEV. 72C 73C INFO* INTEGER 74C on return from CHIEV the value of INFO is 75C = 0 normal return, calculation successful. 76C = K if the eigenvalue iteration fails to converge, 77C eigenvalues (and eigenvectors if requested) 78C 1 through K-1 are correct. 79C 80C Error Messages 81C No. 1 recoverable N is greater than LDA 82C No. 2 recoverable N is less than one. 83C No. 3 recoverable JOB is nonzero and N is greater than LDV 84C No. 4 warning LDA > LDV, elements of A other than the 85C N by N input elements have been changed 86C No. 5 warning LDA < LDV, elements of V other than the 87C N by N output elements have been changed 88C No. 6 recoverable nonreal element on diagonal of A. 89C 90C 91C Subroutines Used 92C 93C EISPACK- HTRIBK, HTRIDI, IMTQL2, TQLRAT 94C BLAS- SCOPY, SCOPYM 95C SLATEC- XERROR 96C***REFERENCES (NONE) 97C***ROUTINES CALLED HTRIBK,HTRIDI,IMTQL2,SCOPY,SCOPYM,TQLRAT,XERROR 98C***END PROLOGUE CHIEV 99 INTEGER I,INFO,J,JOB,K,L,LDA,LDV,M,MDIM,MIN0,N 100 REAL A(*),E(*),WORK(*),V(*) 101C***FIRST EXECUTABLE STATEMENT CHIEV 102 IF(N .GT. LDA) CALL XERROR( 'CHIEV-N .GT. LDA.',17,1,1) 103 IF(N .GT. LDA) RETURN 104 IF(N .LT. 1) CALL XERROR( 'CHIEV-N .LT. 1',14,2,1) 105 IF(N .LT. 1) RETURN 106 IF(N .EQ. 1 .AND. JOB .EQ. 0) GO TO 35 107 MDIM = 2 * LDA 108 IF(JOB .EQ. 0) GO TO 5 109 IF(N .GT. LDV) CALL XERROR( 'CHIEV-JOB .NE. 0, AND N .GT. LDV.', 110 1 33,3,1) 111 IF(N .GT. LDV) RETURN 112 IF(N .EQ. 1) GO TO 35 113C 114C REARRANGE A IF NECESSARY WHEN LDA.GT.LDV AND JOB .NE.0 115C 116 MDIM = MIN0(MDIM,2 * LDV) 117 IF(LDA.LT.LDV) CALL XERROR( 'CHIEV-LDA.LT.LDV, ELEMENTS OF V OTHE 118 1R THAN THE N BY N OUTPUT ELEMENTS HAVE BEEN CHANGED.',89,5,0) 119 IF(LDA.LE.LDV) GO TO 5 120 CALL XERROR( 'CHIEV-LDA.GT.LDV, ELEMENTS OF A OTHER THAN THE N BY 121 1 N INPUT ELEMENTS HAVE BEEN CHANGED.',87,4,0) 122 L = N - 1 123 DO 4 J=1,L 124 M = 1+J*2*LDV 125 K = 1+J*2*LDA 126 CALL SCOPY(2*N,A(K),1,A(M),1) 127 4 CONTINUE 128 5 CONTINUE 129C 130C FILL IN LOWER TRIANGLE OF A, COLUMN BY COLUMN. 131C 132 DO 6 J = 1,N 133 K = (J-1)*(MDIM+2)+1 134 IF(A(K+1).NE.0.0) CALL XERROR 135 1 ( 'CHIEV-NONREAL ELEMENT ON DIAGONAL OF A',38,6,1) 136 IF(A(K+1) .NE.0.0) RETURN 137 CALL SCOPY(N-J+1,A(K),MDIM,A(K),2) 138 CALL SCOPYM(N-J+1,A(K+1),MDIM,A(K+1),2) 139 6 CONTINUE 140C 141C SEPARATE REAL AND IMAGINARY PARTS 142C 143 DO 10 J = 1, N 144 K = (J-1) * MDIM +1 145 L = K + N 146 CALL SCOPY(N,A(K+1),2,WORK(1),1) 147 CALL SCOPY(N,A(K),2,A(K),1) 148 CALL SCOPY(N,WORK(1),1,A(L),1) 149 10 CONTINUE 150C 151C REDUCE A TO TRIDIAGONAL MATRIX. 152C 153 CALL HTRIDI(MDIM,N,A(1),A(N+1),E,WORK(1),WORK(N+1), 154 1 WORK(2*N+1)) 155 IF(JOB .NE. 0) GOTO 15 156C 157C EIGENVALUES ONLY. 158C 159 CALL TQLRAT(N,E,WORK(N+1),INFO) 160 RETURN 161C 162C EIGENVALUES AND EIGENVECTORS. 163C 164 15 DO 17 J = 1,N 165 K = (J-1) * MDIM + 1 166 M = K + N - 1 167 DO 16 I = K,M 168 16 V(I) = 0. 169 I = K + J - 1 170 V(I) = 1. 171 17 CONTINUE 172 CALL IMTQL2(MDIM,N,E,WORK(1),V,INFO) 173 IF(INFO .NE. 0) RETURN 174 CALL HTRIBK(MDIM,N,A(1),A(N+1),WORK(2*N+1),N,V(1),V(N+1)) 175C 176C CONVERT EIGENVECTORS TO COMPLEX STORAGE. 177C 178 DO 20 J = 1,N 179 K = (J-1) * MDIM + 1 180 I = (J-1) * 2 * LDV + 1 181 L = K + N 182 CALL SCOPY(N,V(K),1,WORK(1),1) 183 CALL SCOPY(N,V(L),1,V(I+1),2) 184 CALL SCOPY(N,WORK(1),1,V(I),2) 185 20 CONTINUE 186 RETURN 187C 188C TAKE CARE OF N=1 CASE. 189C 190 35 IF(A(2) .NE. 0.) CALL XERROR 191 1 ( 'CHIEV- NONREAL ELEMENT ON DIAGONAL OF A',39,6,1) 192 IF(A(2) .NE. 0.) RETURN 193 E(1) = A(1) 194 INFO = 0 195 IF(JOB .EQ. 0) RETURN 196 V(1) = A(1) 197 V(2) = 0. 198 RETURN 199 END 200