1      SUBROUTINE CHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
2     $                   LRWORK, IWORK, LIWORK, INFO )
3*
4*  -- LAPACK driver routine (version 3.0) --
5*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6*     Courant Institute, Argonne National Lab, and Rice University
7*     June 30, 1999
8*
9*     .. Scalar Arguments ..
10      CHARACTER          JOBZ, UPLO
11      INTEGER            INFO, LDA, LIWORK, LRWORK, LWORK, N
12*     ..
13*     .. Array Arguments ..
14      INTEGER            IWORK( * )
15      REAL               RWORK( * ), W( * )
16      COMPLEX            A( LDA, * ), WORK( * )
17*     ..
18*
19*  Purpose
20*  =======
21*
22*  CHEEVD computes all eigenvalues and, optionally, eigenvectors of a
23*  complex Hermitian matrix A.  If eigenvectors are desired, it uses a
24*  divide and conquer algorithm.
25*
26*  The divide and conquer algorithm makes very mild assumptions about
27*  floating point arithmetic. It will work on machines with a guard
28*  digit in add/subtract, or on those binary machines without guard
29*  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
30*  Cray-2. It could conceivably fail on hexadecimal or decimal machines
31*  without guard digits, but we know of none.
32*
33*  Arguments
34*  =========
35*
36*  JOBZ    (input) CHARACTER*1
37*          = 'N':  Compute eigenvalues only;
38*          = 'V':  Compute eigenvalues and eigenvectors.
39*
40*  UPLO    (input) CHARACTER*1
41*          = 'U':  Upper triangle of A is stored;
42*          = 'L':  Lower triangle of A is stored.
43*
44*  N       (input) INTEGER
45*          The order of the matrix A.  N >= 0.
46*
47*  A       (input/output) COMPLEX array, dimension (LDA, N)
48*          On entry, the Hermitian matrix A.  If UPLO = 'U', the
49*          leading N-by-N upper triangular part of A contains the
50*          upper triangular part of the matrix A.  If UPLO = 'L',
51*          the leading N-by-N lower triangular part of A contains
52*          the lower triangular part of the matrix A.
53*          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
54*          orthonormal eigenvectors of the matrix A.
55*          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
56*          or the upper triangle (if UPLO='U') of A, including the
57*          diagonal, is destroyed.
58*
59*  LDA     (input) INTEGER
60*          The leading dimension of the array A.  LDA >= max(1,N).
61*
62*  W       (output) REAL array, dimension (N)
63*          If INFO = 0, the eigenvalues in ascending order.
64*
65*  WORK    (workspace/output) COMPLEX array, dimension (LWORK)
66*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
67*
68*  LWORK   (input) INTEGER
69*          The length of the array WORK.
70*          If N <= 1,                LWORK must be at least 1.
71*          If JOBZ  = 'N' and N > 1, LWORK must be at least N + 1.
72*          If JOBZ  = 'V' and N > 1, LWORK must be at least 2*N + N**2.
73*
74*          If LWORK = -1, then a workspace query is assumed; the routine
75*          only calculates the optimal size of the WORK array, returns
76*          this value as the first entry of the WORK array, and no error
77*          message related to LWORK is issued by XERBLA.
78*
79*  RWORK   (workspace/output) REAL array,
80*                                         dimension (LRWORK)
81*          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
82*
83*  LRWORK  (input) INTEGER
84*          The dimension of the array RWORK.
85*          If N <= 1,                LRWORK must be at least 1.
86*          If JOBZ  = 'N' and N > 1, LRWORK must be at least N.
87*          If JOBZ  = 'V' and N > 1, LRWORK must be at least
88*                         1 + 5*N + 2*N**2.
89*
90*          If LRWORK = -1, then a workspace query is assumed; the
91*          routine only calculates the optimal size of the RWORK array,
92*          returns this value as the first entry of the RWORK array, and
93*          no error message related to LRWORK is issued by XERBLA.
94*
95*  IWORK   (workspace/output) INTEGER array, dimension (LIWORK)
96*          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
97*
98*  LIWORK  (input) INTEGER
99*          The dimension of the array IWORK.
100*          If N <= 1,                LIWORK must be at least 1.
101*          If JOBZ  = 'N' and N > 1, LIWORK must be at least 1.
102*          If JOBZ  = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
103*
104*          If LIWORK = -1, then a workspace query is assumed; the
105*          routine only calculates the optimal size of the IWORK array,
106*          returns this value as the first entry of the IWORK array, and
107*          no error message related to LIWORK is issued by XERBLA.
108*
109*  INFO    (output) INTEGER
110*          = 0:  successful exit
111*          < 0:  if INFO = -i, the i-th argument had an illegal value
112*          > 0:  if INFO = i, the algorithm failed to converge; i
113*                off-diagonal elements of an intermediate tridiagonal
114*                form did not converge to zero.
115*
116*  Further Details
117*  ===============
118*
119*  Based on contributions by
120*     Jeff Rutter, Computer Science Division, University of California
121*     at Berkeley, USA
122*
123*  =====================================================================
124*
125*     .. Parameters ..
126      REAL               ZERO, ONE
127      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
128      COMPLEX            CONE
129      PARAMETER          ( CONE = ( 1.0E0, 0.0E0 ) )
130*     ..
131*     .. Local Scalars ..
132      LOGICAL            LOWER, LQUERY, WANTZ
133      INTEGER            IINFO, IMAX, INDE, INDRWK, INDTAU, INDWK2,
134     $                   INDWRK, ISCALE, LIOPT, LIWMIN, LLRWK, LLWORK,
135     $                   LLWRK2, LOPT, LROPT, LRWMIN, LWMIN
136      REAL               ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
137     $                   SMLNUM
138*     ..
139*     .. External Functions ..
140      LOGICAL            LSAME
141      REAL               CLANHE, SLAMCH
142      EXTERNAL           LSAME, CLANHE, SLAMCH
143*     ..
144*     .. External Subroutines ..
145      EXTERNAL           CHETRD, CLACPY, CLASCL, CSTEDC, CUNMTR, SSCAL,
146     $                   SSTERF, XERBLA
147*     ..
148*     .. Intrinsic Functions ..
149      INTRINSIC          INT, MAX, REAL, SQRT
150*     ..
151*     .. Executable Statements ..
152*
153*     Test the input parameters.
154*
155      WANTZ = LSAME( JOBZ, 'V' )
156      LOWER = LSAME( UPLO, 'L' )
157      LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
158*
159      INFO = 0
160      IF( N.LE.1 ) THEN
161         LWMIN = 1
162         LRWMIN = 1
163         LIWMIN = 1
164         LOPT = LWMIN
165         LROPT = LRWMIN
166         LIOPT = LIWMIN
167      ELSE
168         IF( WANTZ ) THEN
169            LWMIN = 2*N + N*N
170            LRWMIN = 1 + 5*N + 2*N**2
171            LIWMIN = 3 + 5*N
172         ELSE
173            LWMIN = N + 1
174            LRWMIN = N
175            LIWMIN = 1
176         END IF
177         LOPT = LWMIN
178         LROPT = LRWMIN
179         LIOPT = LIWMIN
180      END IF
181      IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
182         INFO = -1
183      ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
184         INFO = -2
185      ELSE IF( N.LT.0 ) THEN
186         INFO = -3
187      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
188         INFO = -5
189      ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
190         INFO = -8
191      ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
192         INFO = -10
193      ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
194         INFO = -12
195      END IF
196*
197      IF( INFO.EQ.0 ) THEN
198         WORK( 1 ) = LOPT
199         RWORK( 1 ) = LROPT
200         IWORK( 1 ) = LIOPT
201      END IF
202*
203      IF( INFO.NE.0 ) THEN
204         CALL XERBLA( 'CHEEVD', -INFO )
205         RETURN
206      ELSE IF( LQUERY ) THEN
207         RETURN
208      END IF
209*
210*     Quick return if possible
211*
212      IF( N.EQ.0 )
213     $   RETURN
214*
215      IF( N.EQ.1 ) THEN
216         W( 1 ) = A( 1, 1 )
217         IF( WANTZ )
218     $      A( 1, 1 ) = CONE
219         RETURN
220      END IF
221*
222*     Get machine constants.
223*
224      SAFMIN = SLAMCH( 'Safe minimum' )
225      EPS = SLAMCH( 'Precision' )
226      SMLNUM = SAFMIN / EPS
227      BIGNUM = ONE / SMLNUM
228      RMIN = SQRT( SMLNUM )
229      RMAX = SQRT( BIGNUM )
230*
231*     Scale matrix to allowable range, if necessary.
232*
233      ANRM = CLANHE( 'M', UPLO, N, A, LDA, RWORK )
234      ISCALE = 0
235      IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
236         ISCALE = 1
237         SIGMA = RMIN / ANRM
238      ELSE IF( ANRM.GT.RMAX ) THEN
239         ISCALE = 1
240         SIGMA = RMAX / ANRM
241      END IF
242      IF( ISCALE.EQ.1 )
243     $   CALL CLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
244*
245*     Call CHETRD to reduce Hermitian matrix to tridiagonal form.
246*
247      INDE = 1
248      INDTAU = 1
249      INDWRK = INDTAU + N
250      INDRWK = INDE + N
251      INDWK2 = INDWRK + N*N
252      LLWORK = LWORK - INDWRK + 1
253      LLWRK2 = LWORK - INDWK2 + 1
254      LLRWK = LRWORK - INDRWK + 1
255      CALL CHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ),
256     $             WORK( INDWRK ), LLWORK, IINFO )
257      LOPT = MAX( REAL( LOPT ), REAL( N )+REAL( WORK( INDWRK ) ) )
258*
259*     For eigenvalues only, call SSTERF.  For eigenvectors, first call
260*     CSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
261*     tridiagonal matrix, then call CUNMTR to multiply it to the
262*     Householder transformations represented as Householder vectors in
263*     A.
264*
265      IF( .NOT.WANTZ ) THEN
266         CALL SSTERF( N, W, RWORK( INDE ), INFO )
267      ELSE
268         CALL CSTEDC( 'I', N, W, RWORK( INDE ), WORK( INDWRK ), N,
269     $                WORK( INDWK2 ), LLWRK2, RWORK( INDRWK ), LLRWK,
270     $                IWORK, LIWORK, INFO )
271         CALL CUNMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ),
272     $                WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO )
273         CALL CLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA )
274         LOPT = MAX( LOPT, N+N**2+INT( WORK( INDWK2 ) ) )
275      END IF
276*
277*     If matrix was scaled, then rescale eigenvalues appropriately.
278*
279      IF( ISCALE.EQ.1 ) THEN
280         IF( INFO.EQ.0 ) THEN
281            IMAX = N
282         ELSE
283            IMAX = INFO - 1
284         END IF
285         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
286      END IF
287*
288      WORK( 1 ) = LOPT
289      RWORK( 1 ) = LROPT
290      IWORK( 1 ) = LIOPT
291*
292      RETURN
293*
294*     End of CHEEVD
295*
296      END
297