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