LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
cheevd.f
Go to the documentation of this file.
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
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:86
subroutine chetrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
CHETRD
Definition: chetrd.f:192
subroutine cheevd(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
CHEEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices
Definition: cheevd.f:205
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: clascl.f:143
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
subroutine cunmtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMTR
Definition: cunmtr.f:172
subroutine cstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
CSTEDC
Definition: cstedc.f:212
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79