LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
cheevx.f
Go to the documentation of this file.
1 *> \brief <b> CHEEVX 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 CHEEVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cheevx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cheevx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cheevx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CHEEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
22 * ABSTOL, M, W, Z, LDZ, WORK, LWORK, RWORK,
23 * IWORK, IFAIL, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBZ, RANGE, UPLO
27 * INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N
28 * REAL ABSTOL, VL, VU
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IFAIL( * ), IWORK( * )
32 * REAL RWORK( * ), W( * )
33 * COMPLEX A( LDA, * ), WORK( * ), Z( LDZ, * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> CHEEVX computes selected eigenvalues and, optionally, eigenvectors
43 *> of a complex Hermitian matrix A. Eigenvalues and eigenvectors can
44 *> be selected by specifying either a range of values or a range of
45 *> indices for the desired eigenvalues.
46 *> \endverbatim
47 *
48 * Arguments:
49 * ==========
50 *
51 *> \param[in] JOBZ
52 *> \verbatim
53 *> JOBZ is CHARACTER*1
54 *> = 'N': Compute eigenvalues only;
55 *> = 'V': Compute eigenvalues and eigenvectors.
56 *> \endverbatim
57 *>
58 *> \param[in] RANGE
59 *> \verbatim
60 *> RANGE is CHARACTER*1
61 *> = 'A': all eigenvalues will be found.
62 *> = 'V': all eigenvalues in the half-open interval (VL,VU]
63 *> will be found.
64 *> = 'I': the IL-th through IU-th eigenvalues will be found.
65 *> \endverbatim
66 *>
67 *> \param[in] UPLO
68 *> \verbatim
69 *> UPLO is CHARACTER*1
70 *> = 'U': Upper triangle of A is stored;
71 *> = 'L': Lower triangle of A is stored.
72 *> \endverbatim
73 *>
74 *> \param[in] N
75 *> \verbatim
76 *> N is INTEGER
77 *> The order of the matrix A. N >= 0.
78 *> \endverbatim
79 *>
80 *> \param[in,out] A
81 *> \verbatim
82 *> A is COMPLEX array, dimension (LDA, N)
83 *> On entry, the Hermitian matrix A. If UPLO = 'U', the
84 *> leading N-by-N upper triangular part of A contains the
85 *> upper triangular part of the matrix A. If UPLO = 'L',
86 *> the leading N-by-N lower triangular part of A contains
87 *> the lower triangular part of the matrix A.
88 *> On exit, the lower triangle (if UPLO='L') or the upper
89 *> triangle (if UPLO='U') of A, including the diagonal, is
90 *> destroyed.
91 *> \endverbatim
92 *>
93 *> \param[in] LDA
94 *> \verbatim
95 *> LDA is INTEGER
96 *> The leading dimension of the array A. LDA >= max(1,N).
97 *> \endverbatim
98 *>
99 *> \param[in] VL
100 *> \verbatim
101 *> VL is REAL
102 *> \endverbatim
103 *>
104 *> \param[in] VU
105 *> \verbatim
106 *> VU is REAL
107 *> If RANGE='V', the lower and upper bounds of the interval to
108 *> be searched for eigenvalues. VL < VU.
109 *> Not referenced if RANGE = 'A' or 'I'.
110 *> \endverbatim
111 *>
112 *> \param[in] IL
113 *> \verbatim
114 *> IL is INTEGER
115 *> \endverbatim
116 *>
117 *> \param[in] IU
118 *> \verbatim
119 *> IU is INTEGER
120 *> If RANGE='I', the indices (in ascending order) of the
121 *> smallest and largest eigenvalues to be returned.
122 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
123 *> Not referenced if RANGE = 'A' or 'V'.
124 *> \endverbatim
125 *>
126 *> \param[in] ABSTOL
127 *> \verbatim
128 *> ABSTOL is REAL
129 *> The absolute error tolerance for the eigenvalues.
130 *> An approximate eigenvalue is accepted as converged
131 *> when it is determined to lie in an interval [a,b]
132 *> of width less than or equal to
133 *>
134 *> ABSTOL + EPS * max( |a|,|b| ) ,
135 *>
136 *> where EPS is the machine precision. If ABSTOL is less than
137 *> or equal to zero, then EPS*|T| will be used in its place,
138 *> where |T| is the 1-norm of the tridiagonal matrix obtained
139 *> by reducing A to tridiagonal form.
140 *>
141 *> Eigenvalues will be computed most accurately when ABSTOL is
142 *> set to twice the underflow threshold 2*SLAMCH('S'), not zero.
143 *> If this routine returns with INFO>0, indicating that some
144 *> eigenvectors did not converge, try setting ABSTOL to
145 *> 2*SLAMCH('S').
146 *>
147 *> See "Computing Small Singular Values of Bidiagonal Matrices
148 *> with Guaranteed High Relative Accuracy," by Demmel and
149 *> Kahan, LAPACK Working Note #3.
150 *> \endverbatim
151 *>
152 *> \param[out] M
153 *> \verbatim
154 *> M is INTEGER
155 *> The total number of eigenvalues found. 0 <= M <= N.
156 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
157 *> \endverbatim
158 *>
159 *> \param[out] W
160 *> \verbatim
161 *> W is REAL array, dimension (N)
162 *> On normal exit, the first M elements contain the selected
163 *> eigenvalues in ascending order.
164 *> \endverbatim
165 *>
166 *> \param[out] Z
167 *> \verbatim
168 *> Z is COMPLEX array, dimension (LDZ, max(1,M))
169 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
170 *> contain the orthonormal eigenvectors of the matrix A
171 *> corresponding to the selected eigenvalues, with the i-th
172 *> column of Z holding the eigenvector associated with W(i).
173 *> If an eigenvector fails to converge, then that column of Z
174 *> contains the latest approximation to the eigenvector, and the
175 *> index of the eigenvector is returned in IFAIL.
176 *> If JOBZ = 'N', then Z is not referenced.
177 *> Note: the user must ensure that at least max(1,M) columns are
178 *> supplied in the array Z; if RANGE = 'V', the exact value of M
179 *> is not known in advance and an upper bound must be used.
180 *> \endverbatim
181 *>
182 *> \param[in] LDZ
183 *> \verbatim
184 *> LDZ is INTEGER
185 *> The leading dimension of the array Z. LDZ >= 1, and if
186 *> JOBZ = 'V', LDZ >= max(1,N).
187 *> \endverbatim
188 *>
189 *> \param[out] WORK
190 *> \verbatim
191 *> WORK is COMPLEX array, dimension (MAX(1,LWORK))
192 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
193 *> \endverbatim
194 *>
195 *> \param[in] LWORK
196 *> \verbatim
197 *> LWORK is INTEGER
198 *> The length of the array WORK. LWORK >= 1, when N <= 1;
199 *> otherwise 2*N.
200 *> For optimal efficiency, LWORK >= (NB+1)*N,
201 *> where NB is the max of the blocksize for CHETRD and for
202 *> CUNMTR as returned by ILAENV.
203 *>
204 *> If LWORK = -1, then a workspace query is assumed; the routine
205 *> only calculates the optimal size of the WORK array, returns
206 *> this value as the first entry of the WORK array, and no error
207 *> message related to LWORK is issued by XERBLA.
208 *> \endverbatim
209 *>
210 *> \param[out] RWORK
211 *> \verbatim
212 *> RWORK is REAL array, dimension (7*N)
213 *> \endverbatim
214 *>
215 *> \param[out] IWORK
216 *> \verbatim
217 *> IWORK is INTEGER array, dimension (5*N)
218 *> \endverbatim
219 *>
220 *> \param[out] IFAIL
221 *> \verbatim
222 *> IFAIL is INTEGER array, dimension (N)
223 *> If JOBZ = 'V', then if INFO = 0, the first M elements of
224 *> IFAIL are zero. If INFO > 0, then IFAIL contains the
225 *> indices of the eigenvectors that failed to converge.
226 *> If JOBZ = 'N', then IFAIL is not referenced.
227 *> \endverbatim
228 *>
229 *> \param[out] INFO
230 *> \verbatim
231 *> INFO is INTEGER
232 *> = 0: successful exit
233 *> < 0: if INFO = -i, the i-th argument had an illegal value
234 *> > 0: if INFO = i, then i eigenvectors failed to converge.
235 *> Their indices are stored in array IFAIL.
236 *> \endverbatim
237 *
238 * Authors:
239 * ========
240 *
241 *> \author Univ. of Tennessee
242 *> \author Univ. of California Berkeley
243 *> \author Univ. of Colorado Denver
244 *> \author NAG Ltd.
245 *
246 *> \date November 2011
247 *
248 *> \ingroup complexHEeigen
249 *
250 * =====================================================================
251  SUBROUTINE cheevx( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
252  $ abstol, m, w, z, ldz, work, lwork, rwork,
253  $ iwork, ifail, info )
254 *
255 * -- LAPACK driver routine (version 3.4.0) --
256 * -- LAPACK is a software package provided by Univ. of Tennessee, --
257 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
258 * November 2011
259 *
260 * .. Scalar Arguments ..
261  CHARACTER jobz, range, uplo
262  INTEGER il, info, iu, lda, ldz, lwork, m, n
263  REAL abstol, vl, vu
264 * ..
265 * .. Array Arguments ..
266  INTEGER ifail( * ), iwork( * )
267  REAL rwork( * ), w( * )
268  COMPLEX a( lda, * ), work( * ), z( ldz, * )
269 * ..
270 *
271 * =====================================================================
272 *
273 * .. Parameters ..
274  REAL zero, one
275  parameter( zero = 0.0e+0, one = 1.0e+0 )
276  COMPLEX cone
277  parameter( cone = ( 1.0e+0, 0.0e+0 ) )
278 * ..
279 * .. Local Scalars ..
280  LOGICAL alleig, indeig, lower, lquery, test, valeig,
281  $ wantz
282  CHARACTER order
283  INTEGER i, iinfo, imax, indd, inde, indee, indibl,
284  $ indisp, indiwk, indrwk, indtau, indwrk, iscale,
285  $ itmp1, j, jj, llwork, lwkmin, lwkopt, nb,
286  $ nsplit
287  REAL abstll, anrm, bignum, eps, rmax, rmin, safmin,
288  $ sigma, smlnum, tmp1, vll, vuu
289 * ..
290 * .. External Functions ..
291  LOGICAL lsame
292  INTEGER ilaenv
293  REAL slamch, clanhe
294  EXTERNAL lsame, ilaenv, slamch, clanhe
295 * ..
296 * .. External Subroutines ..
297  EXTERNAL scopy, sscal, sstebz, ssterf, xerbla, csscal,
299  $ cunmtr
300 * ..
301 * .. Intrinsic Functions ..
302  INTRINSIC REAL, max, min, sqrt
303 * ..
304 * .. Executable Statements ..
305 *
306 * Test the input parameters.
307 *
308  lower = lsame( uplo, 'L' )
309  wantz = lsame( jobz, 'V' )
310  alleig = lsame( range, 'A' )
311  valeig = lsame( range, 'V' )
312  indeig = lsame( range, 'I' )
313  lquery = ( lwork.EQ.-1 )
314 *
315  info = 0
316  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
317  info = -1
318  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
319  info = -2
320  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
321  info = -3
322  ELSE IF( n.LT.0 ) THEN
323  info = -4
324  ELSE IF( lda.LT.max( 1, n ) ) THEN
325  info = -6
326  ELSE
327  IF( valeig ) THEN
328  IF( n.GT.0 .AND. vu.LE.vl )
329  $ info = -8
330  ELSE IF( indeig ) THEN
331  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
332  info = -9
333  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
334  info = -10
335  END IF
336  END IF
337  END IF
338  IF( info.EQ.0 ) THEN
339  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
340  info = -15
341  END IF
342  END IF
343 *
344  IF( info.EQ.0 ) THEN
345  IF( n.LE.1 ) THEN
346  lwkmin = 1
347  work( 1 ) = lwkmin
348  ELSE
349  lwkmin = 2*n
350  nb = ilaenv( 1, 'CHETRD', uplo, n, -1, -1, -1 )
351  nb = max( nb, ilaenv( 1, 'CUNMTR', uplo, n, -1, -1, -1 ) )
352  lwkopt = max( 1, ( nb + 1 )*n )
353  work( 1 ) = lwkopt
354  END IF
355 *
356  IF( lwork.LT.lwkmin .AND. .NOT.lquery )
357  $ info = -17
358  END IF
359 *
360  IF( info.NE.0 ) THEN
361  CALL xerbla( 'CHEEVX', -info )
362  RETURN
363  ELSE IF( lquery ) THEN
364  RETURN
365  END IF
366 *
367 * Quick return if possible
368 *
369  m = 0
370  IF( n.EQ.0 ) THEN
371  RETURN
372  END IF
373 *
374  IF( n.EQ.1 ) THEN
375  IF( alleig .OR. indeig ) THEN
376  m = 1
377  w( 1 ) = a( 1, 1 )
378  ELSE IF( valeig ) THEN
379  IF( vl.LT.REAL( A( 1, 1 ) ) .AND. vu.GE.REAL( A( 1, 1 ) ) )
380  $ THEN
381  m = 1
382  w( 1 ) = a( 1, 1 )
383  END IF
384  END IF
385  IF( wantz )
386  $ z( 1, 1 ) = cone
387  RETURN
388  END IF
389 *
390 * Get machine constants.
391 *
392  safmin = slamch( 'Safe minimum' )
393  eps = slamch( 'Precision' )
394  smlnum = safmin / eps
395  bignum = one / smlnum
396  rmin = sqrt( smlnum )
397  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
398 *
399 * Scale matrix to allowable range, if necessary.
400 *
401  iscale = 0
402  abstll = abstol
403  IF( valeig ) THEN
404  vll = vl
405  vuu = vu
406  END IF
407  anrm = clanhe( 'M', uplo, n, a, lda, rwork )
408  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
409  iscale = 1
410  sigma = rmin / anrm
411  ELSE IF( anrm.GT.rmax ) THEN
412  iscale = 1
413  sigma = rmax / anrm
414  END IF
415  IF( iscale.EQ.1 ) THEN
416  IF( lower ) THEN
417  DO 10 j = 1, n
418  CALL csscal( n-j+1, sigma, a( j, j ), 1 )
419  10 CONTINUE
420  ELSE
421  DO 20 j = 1, n
422  CALL csscal( j, sigma, a( 1, j ), 1 )
423  20 CONTINUE
424  END IF
425  IF( abstol.GT.0 )
426  $ abstll = abstol*sigma
427  IF( valeig ) THEN
428  vll = vl*sigma
429  vuu = vu*sigma
430  END IF
431  END IF
432 *
433 * Call CHETRD to reduce Hermitian matrix to tridiagonal form.
434 *
435  indd = 1
436  inde = indd + n
437  indrwk = inde + n
438  indtau = 1
439  indwrk = indtau + n
440  llwork = lwork - indwrk + 1
441  CALL chetrd( uplo, n, a, lda, rwork( indd ), rwork( inde ),
442  $ work( indtau ), work( indwrk ), llwork, iinfo )
443 *
444 * If all eigenvalues are desired and ABSTOL is less than or equal to
445 * zero, then call SSTERF or CUNGTR and CSTEQR. If this fails for
446 * some eigenvalue, then try SSTEBZ.
447 *
448  test = .false.
449  IF( indeig ) THEN
450  IF( il.EQ.1 .AND. iu.EQ.n ) THEN
451  test = .true.
452  END IF
453  END IF
454  IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
455  CALL scopy( n, rwork( indd ), 1, w, 1 )
456  indee = indrwk + 2*n
457  IF( .NOT.wantz ) THEN
458  CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
459  CALL ssterf( n, w, rwork( indee ), info )
460  ELSE
461  CALL clacpy( 'A', n, n, a, lda, z, ldz )
462  CALL cungtr( uplo, n, z, ldz, work( indtau ),
463  $ work( indwrk ), llwork, iinfo )
464  CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
465  CALL csteqr( jobz, n, w, rwork( indee ), z, ldz,
466  $ rwork( indrwk ), info )
467  IF( info.EQ.0 ) THEN
468  DO 30 i = 1, n
469  ifail( i ) = 0
470  30 CONTINUE
471  END IF
472  END IF
473  IF( info.EQ.0 ) THEN
474  m = n
475  go to 40
476  END IF
477  info = 0
478  END IF
479 *
480 * Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN.
481 *
482  IF( wantz ) THEN
483  order = 'B'
484  ELSE
485  order = 'E'
486  END IF
487  indibl = 1
488  indisp = indibl + n
489  indiwk = indisp + n
490  CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
491  $ rwork( indd ), rwork( inde ), m, nsplit, w,
492  $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
493  $ iwork( indiwk ), info )
494 *
495  IF( wantz ) THEN
496  CALL cstein( n, rwork( indd ), rwork( inde ), m, w,
497  $ iwork( indibl ), iwork( indisp ), z, ldz,
498  $ rwork( indrwk ), iwork( indiwk ), ifail, info )
499 *
500 * Apply unitary matrix used in reduction to tridiagonal
501 * form to eigenvectors returned by CSTEIN.
502 *
503  CALL cunmtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
504  $ ldz, work( indwrk ), llwork, iinfo )
505  END IF
506 *
507 * If matrix was scaled, then rescale eigenvalues appropriately.
508 *
509  40 CONTINUE
510  IF( iscale.EQ.1 ) THEN
511  IF( info.EQ.0 ) THEN
512  imax = m
513  ELSE
514  imax = info - 1
515  END IF
516  CALL sscal( imax, one / sigma, w, 1 )
517  END IF
518 *
519 * If eigenvalues are not in order, then sort them, along with
520 * eigenvectors.
521 *
522  IF( wantz ) THEN
523  DO 60 j = 1, m - 1
524  i = 0
525  tmp1 = w( j )
526  DO 50 jj = j + 1, m
527  IF( w( jj ).LT.tmp1 ) THEN
528  i = jj
529  tmp1 = w( jj )
530  END IF
531  50 CONTINUE
532 *
533  IF( i.NE.0 ) THEN
534  itmp1 = iwork( indibl+i-1 )
535  w( i ) = w( j )
536  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
537  w( j ) = tmp1
538  iwork( indibl+j-1 ) = itmp1
539  CALL cswap( n, z( 1, i ), 1, z( 1, j ), 1 )
540  IF( info.NE.0 ) THEN
541  itmp1 = ifail( i )
542  ifail( i ) = ifail( j )
543  ifail( j ) = itmp1
544  END IF
545  END IF
546  60 CONTINUE
547  END IF
548 *
549 * Set WORK(1) to optimal complex workspace size.
550 *
551  work( 1 ) = lwkopt
552 *
553  RETURN
554 *
555 * End of CHEEVX
556 *
557  END
real function clanhe(NORM, UPLO, N, A, LDA, WORK)
CLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian matrix.
Definition: clanhe.f:125
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:61
subroutine cunmtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMTR
Definition: cunmtr.f:172
subroutine cungtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
CUNGTR
Definition: cungtr.f:124
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:51
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:52
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:54
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:104
subroutine csteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
CSTEQR
Definition: csteqr.f:133
subroutine cstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
CSTEIN
Definition: cstein.f:182
subroutine sstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
SSTEBZ
Definition: sstebz.f:262
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
Definition: xerbla-fortran:9
subroutine cheevx(JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, LWORK, RWORK, IWORK, IFAIL, INFO)
CHEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices ...
Definition: cheevx.f:251
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:81
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:53
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:87
subroutine chetrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
CHETRD
Definition: chetrd.f:193
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:54