1*> \brief <b> CHEEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices</b>
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CHEEVD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cheevd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cheevd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cheevd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
22*                          LRWORK, IWORK, LIWORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          JOBZ, UPLO
26*       INTEGER            INFO, LDA, LIWORK, LRWORK, LWORK, N
27*       ..
28*       .. Array Arguments ..
29*       INTEGER            IWORK( * )
30*       REAL               RWORK( * ), W( * )
31*       COMPLEX            A( LDA, * ), WORK( * )
32*       ..
33*
34*
35*> \par Purpose:
36*  =============
37*>
38*> \verbatim
39*>
40*> CHEEVD computes all eigenvalues and, optionally, eigenvectors of a
41*> complex Hermitian matrix A.  If eigenvectors are desired, it uses a
42*> divide and conquer algorithm.
43*>
44*> The divide and conquer algorithm makes very mild assumptions about
45*> floating point arithmetic. It will work on machines with a guard
46*> digit in add/subtract, or on those binary machines without guard
47*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
48*> Cray-2. It could conceivably fail on hexadecimal or decimal machines
49*> without guard digits, but we know of none.
50*> \endverbatim
51*
52*  Arguments:
53*  ==========
54*
55*> \param[in] JOBZ
56*> \verbatim
57*>          JOBZ is CHARACTER*1
58*>          = 'N':  Compute eigenvalues only;
59*>          = 'V':  Compute eigenvalues and eigenvectors.
60*> \endverbatim
61*>
62*> \param[in] UPLO
63*> \verbatim
64*>          UPLO is CHARACTER*1
65*>          = 'U':  Upper triangle of A is stored;
66*>          = 'L':  Lower triangle of A is stored.
67*> \endverbatim
68*>
69*> \param[in] N
70*> \verbatim
71*>          N is INTEGER
72*>          The order of the matrix A.  N >= 0.
73*> \endverbatim
74*>
75*> \param[in,out] A
76*> \verbatim
77*>          A is COMPLEX array, dimension (LDA, N)
78*>          On entry, the Hermitian matrix A.  If UPLO = 'U', the
79*>          leading N-by-N upper triangular part of A contains the
80*>          upper triangular part of the matrix A.  If UPLO = 'L',
81*>          the leading N-by-N lower triangular part of A contains
82*>          the lower triangular part of the matrix A.
83*>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
84*>          orthonormal eigenvectors of the matrix A.
85*>          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
86*>          or the upper triangle (if UPLO='U') of A, including the
87*>          diagonal, is destroyed.
88*> \endverbatim
89*>
90*> \param[in] LDA
91*> \verbatim
92*>          LDA is INTEGER
93*>          The leading dimension of the array A.  LDA >= max(1,N).
94*> \endverbatim
95*>
96*> \param[out] W
97*> \verbatim
98*>          W is REAL array, dimension (N)
99*>          If INFO = 0, the eigenvalues in ascending order.
100*> \endverbatim
101*>
102*> \param[out] WORK
103*> \verbatim
104*>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
105*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
106*> \endverbatim
107*>
108*> \param[in] LWORK
109*> \verbatim
110*>          LWORK is INTEGER
111*>          The length of the array WORK.
112*>          If N <= 1,                LWORK must be at least 1.
113*>          If JOBZ  = 'N' and N > 1, LWORK must be at least N + 1.
114*>          If JOBZ  = 'V' and N > 1, LWORK must be at least 2*N + N**2.
115*>
116*>          If LWORK = -1, then a workspace query is assumed; the routine
117*>          only calculates the optimal sizes of the WORK, RWORK and
118*>          IWORK arrays, returns these values as the first entries of
119*>          the WORK, RWORK and IWORK arrays, and no error message
120*>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
121*> \endverbatim
122*>
123*> \param[out] RWORK
124*> \verbatim
125*>          RWORK is REAL array,
126*>                                         dimension (LRWORK)
127*>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
128*> \endverbatim
129*>
130*> \param[in] LRWORK
131*> \verbatim
132*>          LRWORK is INTEGER
133*>          The dimension of the array RWORK.
134*>          If N <= 1,                LRWORK must be at least 1.
135*>          If JOBZ  = 'N' and N > 1, LRWORK must be at least N.
136*>          If JOBZ  = 'V' and N > 1, LRWORK must be at least
137*>                         1 + 5*N + 2*N**2.
138*>
139*>          If LRWORK = -1, then a workspace query is assumed; the
140*>          routine only calculates the optimal sizes of the WORK, RWORK
141*>          and IWORK arrays, returns these values as the first entries
142*>          of the WORK, RWORK and IWORK arrays, and no error message
143*>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
144*> \endverbatim
145*>
146*> \param[out] IWORK
147*> \verbatim
148*>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
149*>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
150*> \endverbatim
151*>
152*> \param[in] LIWORK
153*> \verbatim
154*>          LIWORK is INTEGER
155*>          The dimension of the array IWORK.
156*>          If N <= 1,                LIWORK must be at least 1.
157*>          If JOBZ  = 'N' and N > 1, LIWORK must be at least 1.
158*>          If JOBZ  = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
159*>
160*>          If LIWORK = -1, then a workspace query is assumed; the
161*>          routine only calculates the optimal sizes of the WORK, RWORK
162*>          and IWORK arrays, returns these values as the first entries
163*>          of the WORK, RWORK and IWORK arrays, and no error message
164*>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
165*> \endverbatim
166*>
167*> \param[out] INFO
168*> \verbatim
169*>          INFO is INTEGER
170*>          = 0:  successful exit
171*>          < 0:  if INFO = -i, the i-th argument had an illegal value
172*>          > 0:  if INFO = i and JOBZ = 'N', then the algorithm failed
173*>                to converge; i off-diagonal elements of an intermediate
174*>                tridiagonal form did not converge to zero;
175*>                if INFO = i and JOBZ = 'V', then the algorithm failed
176*>                to compute an eigenvalue while working on the submatrix
177*>                lying in rows and columns INFO/(N+1) through
178*>                mod(INFO,N+1).
179*> \endverbatim
180*
181*  Authors:
182*  ========
183*
184*> \author Univ. of Tennessee
185*> \author Univ. of California Berkeley
186*> \author Univ. of Colorado Denver
187*> \author NAG Ltd.
188*
189*> \ingroup complexHEeigen
190*
191*> \par Further Details:
192*  =====================
193*>
194*>  Modified description of INFO. Sven, 16 Feb 05.
195*
196*> \par Contributors:
197*  ==================
198*>
199*> Jeff Rutter, Computer Science Division, University of California
200*> at Berkeley, USA
201*>
202*  =====================================================================
203      SUBROUTINE CHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
204     $                   LRWORK, IWORK, LIWORK, INFO )
205*
206*  -- LAPACK driver routine --
207*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
208*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
209*
210*     .. Scalar Arguments ..
211      CHARACTER          JOBZ, UPLO
212      INTEGER            INFO, LDA, LIWORK, LRWORK, LWORK, N
213*     ..
214*     .. Array Arguments ..
215      INTEGER            IWORK( * )
216      REAL               RWORK( * ), W( * )
217      COMPLEX            A( LDA, * ), WORK( * )
218*     ..
219*
220*  =====================================================================
221*
222*     .. Parameters ..
223      REAL               ZERO, ONE
224      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
225      COMPLEX            CONE
226      PARAMETER          ( CONE = ( 1.0E0, 0.0E0 ) )
227*     ..
228*     .. Local Scalars ..
229      LOGICAL            LOWER, LQUERY, WANTZ
230      INTEGER            IINFO, IMAX, INDE, INDRWK, INDTAU, INDWK2,
231     $                   INDWRK, ISCALE, LIOPT, LIWMIN, LLRWK, LLWORK,
232     $                   LLWRK2, LOPT, LROPT, LRWMIN, LWMIN
233      REAL               ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
234     $                   SMLNUM
235*     ..
236*     .. External Functions ..
237      LOGICAL            LSAME
238      INTEGER            ILAENV
239      REAL               CLANHE, SLAMCH
240      EXTERNAL           ILAENV, LSAME, CLANHE, SLAMCH
241*     ..
242*     .. External Subroutines ..
243      EXTERNAL           CHETRD, CLACPY, CLASCL, CSTEDC, CUNMTR, SSCAL,
244     $                   SSTERF, XERBLA
245*     ..
246*     .. Intrinsic Functions ..
247      INTRINSIC          MAX, SQRT
248*     ..
249*     .. Executable Statements ..
250*
251*     Test the input parameters.
252*
253      WANTZ = LSAME( JOBZ, 'V' )
254      LOWER = LSAME( UPLO, 'L' )
255      LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
256*
257      INFO = 0
258      IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
259         INFO = -1
260      ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
261         INFO = -2
262      ELSE IF( N.LT.0 ) THEN
263         INFO = -3
264      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
265         INFO = -5
266      END IF
267*
268      IF( INFO.EQ.0 ) THEN
269         IF( N.LE.1 ) THEN
270            LWMIN = 1
271            LRWMIN = 1
272            LIWMIN = 1
273            LOPT = LWMIN
274            LROPT = LRWMIN
275            LIOPT = LIWMIN
276         ELSE
277            IF( WANTZ ) THEN
278               LWMIN = 2*N + N*N
279               LRWMIN = 1 + 5*N + 2*N**2
280               LIWMIN = 3 + 5*N
281            ELSE
282               LWMIN = N + 1
283               LRWMIN = N
284               LIWMIN = 1
285            END IF
286            LOPT = MAX( LWMIN, N +
287     $                  ILAENV( 1, 'CHETRD', UPLO, N, -1, -1, -1 ) )
288            LROPT = LRWMIN
289            LIOPT = LIWMIN
290         END IF
291         WORK( 1 ) = LOPT
292         RWORK( 1 ) = LROPT
293         IWORK( 1 ) = LIOPT
294*
295         IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
296            INFO = -8
297         ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
298            INFO = -10
299         ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
300            INFO = -12
301         END IF
302      END IF
303*
304      IF( INFO.NE.0 ) THEN
305         CALL XERBLA( 'CHEEVD', -INFO )
306         RETURN
307      ELSE IF( LQUERY ) THEN
308         RETURN
309      END IF
310*
311*     Quick return if possible
312*
313      IF( N.EQ.0 )
314     $   RETURN
315*
316      IF( N.EQ.1 ) THEN
317         W( 1 ) = REAL( A( 1, 1 ) )
318         IF( WANTZ )
319     $      A( 1, 1 ) = CONE
320         RETURN
321      END IF
322*
323*     Get machine constants.
324*
325      SAFMIN = SLAMCH( 'Safe minimum' )
326      EPS = SLAMCH( 'Precision' )
327      SMLNUM = SAFMIN / EPS
328      BIGNUM = ONE / SMLNUM
329      RMIN = SQRT( SMLNUM )
330      RMAX = SQRT( BIGNUM )
331*
332*     Scale matrix to allowable range, if necessary.
333*
334      ANRM = CLANHE( 'M', UPLO, N, A, LDA, RWORK )
335      ISCALE = 0
336      IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
337         ISCALE = 1
338         SIGMA = RMIN / ANRM
339      ELSE IF( ANRM.GT.RMAX ) THEN
340         ISCALE = 1
341         SIGMA = RMAX / ANRM
342      END IF
343      IF( ISCALE.EQ.1 )
344     $   CALL CLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
345*
346*     Call CHETRD to reduce Hermitian matrix to tridiagonal form.
347*
348      INDE = 1
349      INDTAU = 1
350      INDWRK = INDTAU + N
351      INDRWK = INDE + N
352      INDWK2 = INDWRK + N*N
353      LLWORK = LWORK - INDWRK + 1
354      LLWRK2 = LWORK - INDWK2 + 1
355      LLRWK = LRWORK - INDRWK + 1
356      CALL CHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ),
357     $             WORK( INDWRK ), LLWORK, IINFO )
358*
359*     For eigenvalues only, call SSTERF.  For eigenvectors, first call
360*     CSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
361*     tridiagonal matrix, then call CUNMTR to multiply it to the
362*     Householder transformations represented as Householder vectors in
363*     A.
364*
365      IF( .NOT.WANTZ ) THEN
366         CALL SSTERF( N, W, RWORK( INDE ), INFO )
367      ELSE
368         CALL CSTEDC( 'I', N, W, RWORK( INDE ), WORK( INDWRK ), N,
369     $                WORK( INDWK2 ), LLWRK2, RWORK( INDRWK ), LLRWK,
370     $                IWORK, LIWORK, INFO )
371         CALL CUNMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ),
372     $                WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO )
373         CALL CLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA )
374      END IF
375*
376*     If matrix was scaled, then rescale eigenvalues appropriately.
377*
378      IF( ISCALE.EQ.1 ) THEN
379         IF( INFO.EQ.0 ) THEN
380            IMAX = N
381         ELSE
382            IMAX = INFO - 1
383         END IF
384         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
385      END IF
386*
387      WORK( 1 ) = LOPT
388      RWORK( 1 ) = LROPT
389      IWORK( 1 ) = LIOPT
390*
391      RETURN
392*
393*     End of CHEEVD
394*
395      END
396