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*> \date December 2016
190*
191*> \ingroup complexHEeigen
192*
193*> \par Further Details:
194*  =====================
195*>
196*>  Modified description of INFO. Sven, 16 Feb 05.
197*
198*> \par Contributors:
199*  ==================
200*>
201*> Jeff Rutter, Computer Science Division, University of California
202*> at Berkeley, USA
203*>
204*  =====================================================================
205      SUBROUTINE CHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
206     $                   LRWORK, IWORK, LIWORK, INFO )
207*
208*  -- LAPACK driver routine (version 3.7.0) --
209*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
210*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
211*     December 2016
212*
213*     .. Scalar Arguments ..
214      CHARACTER          JOBZ, UPLO
215      INTEGER            INFO, LDA, LIWORK, LRWORK, LWORK, N
216*     ..
217*     .. Array Arguments ..
218      INTEGER            IWORK( * )
219      REAL               RWORK( * ), W( * )
220      COMPLEX            A( LDA, * ), WORK( * )
221*     ..
222*
223*  =====================================================================
224*
225*     .. Parameters ..
226      REAL               ZERO, ONE
227      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
228      COMPLEX            CONE
229      PARAMETER          ( CONE = ( 1.0E0, 0.0E0 ) )
230*     ..
231*     .. Local Scalars ..
232      LOGICAL            LOWER, LQUERY, WANTZ
233      INTEGER            IINFO, IMAX, INDE, INDRWK, INDTAU, INDWK2,
234     $                   INDWRK, ISCALE, LIOPT, LIWMIN, LLRWK, LLWORK,
235     $                   LLWRK2, LOPT, LROPT, LRWMIN, LWMIN
236      REAL               ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
237     $                   SMLNUM
238*     ..
239*     .. External Functions ..
240      LOGICAL            LSAME
241      INTEGER            ILAENV
242      REAL               CLANHE, SLAMCH
243      EXTERNAL           ILAENV, LSAME, CLANHE, SLAMCH
244*     ..
245*     .. External Subroutines ..
246      EXTERNAL           CHETRD, CLACPY, CLASCL, CSTEDC, CUNMTR, SSCAL,
247     $                   SSTERF, XERBLA
248*     ..
249*     .. Intrinsic Functions ..
250      INTRINSIC          MAX, SQRT
251*     ..
252*     .. Executable Statements ..
253*
254*     Test the input parameters.
255*
256      WANTZ = LSAME( JOBZ, 'V' )
257      LOWER = LSAME( UPLO, 'L' )
258      LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
259*
260      INFO = 0
261      IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
262         INFO = -1
263      ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
264         INFO = -2
265      ELSE IF( N.LT.0 ) THEN
266         INFO = -3
267      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
268         INFO = -5
269      END IF
270*
271      IF( INFO.EQ.0 ) THEN
272         IF( N.LE.1 ) THEN
273            LWMIN = 1
274            LRWMIN = 1
275            LIWMIN = 1
276            LOPT = LWMIN
277            LROPT = LRWMIN
278            LIOPT = LIWMIN
279         ELSE
280            IF( WANTZ ) THEN
281               LWMIN = 2*N + N*N
282               LRWMIN = 1 + 5*N + 2*N**2
283               LIWMIN = 3 + 5*N
284            ELSE
285               LWMIN = N + 1
286               LRWMIN = N
287               LIWMIN = 1
288            END IF
289            LOPT = MAX( LWMIN, N +
290     $                  ILAENV( 1, 'CHETRD', UPLO, N, -1, -1, -1 ) )
291            LROPT = LRWMIN
292            LIOPT = LIWMIN
293         END IF
294         WORK( 1 ) = LOPT
295         RWORK( 1 ) = LROPT
296         IWORK( 1 ) = LIOPT
297*
298         IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
299            INFO = -8
300         ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
301            INFO = -10
302         ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
303            INFO = -12
304         END IF
305      END IF
306*
307      IF( INFO.NE.0 ) THEN
308         CALL XERBLA( 'CHEEVD', -INFO )
309         RETURN
310      ELSE IF( LQUERY ) THEN
311         RETURN
312      END IF
313*
314*     Quick return if possible
315*
316      IF( N.EQ.0 )
317     $   RETURN
318*
319      IF( N.EQ.1 ) THEN
320         W( 1 ) = A( 1, 1 )
321         IF( WANTZ )
322     $      A( 1, 1 ) = CONE
323         RETURN
324      END IF
325*
326*     Get machine constants.
327*
328      SAFMIN = SLAMCH( 'Safe minimum' )
329      EPS = SLAMCH( 'Precision' )
330      SMLNUM = SAFMIN / EPS
331      BIGNUM = ONE / SMLNUM
332      RMIN = SQRT( SMLNUM )
333      RMAX = SQRT( BIGNUM )
334*
335*     Scale matrix to allowable range, if necessary.
336*
337      ANRM = CLANHE( 'M', UPLO, N, A, LDA, RWORK )
338      ISCALE = 0
339      IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
340         ISCALE = 1
341         SIGMA = RMIN / ANRM
342      ELSE IF( ANRM.GT.RMAX ) THEN
343         ISCALE = 1
344         SIGMA = RMAX / ANRM
345      END IF
346      IF( ISCALE.EQ.1 )
347     $   CALL CLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
348*
349*     Call CHETRD to reduce Hermitian matrix to tridiagonal form.
350*
351      INDE = 1
352      INDTAU = 1
353      INDWRK = INDTAU + N
354      INDRWK = INDE + N
355      INDWK2 = INDWRK + N*N
356      LLWORK = LWORK - INDWRK + 1
357      LLWRK2 = LWORK - INDWK2 + 1
358      LLRWK = LRWORK - INDRWK + 1
359      CALL CHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ),
360     $             WORK( INDWRK ), LLWORK, IINFO )
361*
362*     For eigenvalues only, call SSTERF.  For eigenvectors, first call
363*     CSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
364*     tridiagonal matrix, then call CUNMTR to multiply it to the
365*     Householder transformations represented as Householder vectors in
366*     A.
367*
368      IF( .NOT.WANTZ ) THEN
369         CALL SSTERF( N, W, RWORK( INDE ), INFO )
370      ELSE
371         CALL CSTEDC( 'I', N, W, RWORK( INDE ), WORK( INDWRK ), N,
372     $                WORK( INDWK2 ), LLWRK2, RWORK( INDRWK ), LLRWK,
373     $                IWORK, LIWORK, INFO )
374         CALL CUNMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ),
375     $                WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO )
376         CALL CLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA )
377      END IF
378*
379*     If matrix was scaled, then rescale eigenvalues appropriately.
380*
381      IF( ISCALE.EQ.1 ) THEN
382         IF( INFO.EQ.0 ) THEN
383            IMAX = N
384         ELSE
385            IMAX = INFO - 1
386         END IF
387         CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
388      END IF
389*
390      WORK( 1 ) = LOPT
391      RWORK( 1 ) = LROPT
392      IWORK( 1 ) = LIOPT
393*
394      RETURN
395*
396*     End of CHEEVD
397*
398      END
399