LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
zhbevx.f
Go to the documentation of this file.
1 *> \brief <b> ZHBEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER 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 ZHBEVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhbevx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhbevx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhbevx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ, VL,
22 * VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, RWORK,
23 * IWORK, IFAIL, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBZ, RANGE, UPLO
27 * INTEGER IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N
28 * DOUBLE PRECISION ABSTOL, VL, VU
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IFAIL( * ), IWORK( * )
32 * DOUBLE PRECISION RWORK( * ), W( * )
33 * COMPLEX*16 AB( LDAB, * ), Q( LDQ, * ), WORK( * ),
34 * $ Z( LDZ, * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> ZHBEVX computes selected eigenvalues and, optionally, eigenvectors
44 *> of a complex Hermitian band matrix A. Eigenvalues and eigenvectors
45 *> can be selected by specifying either a range of values or a range of
46 *> indices for the desired eigenvalues.
47 *> \endverbatim
48 *
49 * Arguments:
50 * ==========
51 *
52 *> \param[in] JOBZ
53 *> \verbatim
54 *> JOBZ is CHARACTER*1
55 *> = 'N': Compute eigenvalues only;
56 *> = 'V': Compute eigenvalues and eigenvectors.
57 *> \endverbatim
58 *>
59 *> \param[in] RANGE
60 *> \verbatim
61 *> RANGE is CHARACTER*1
62 *> = 'A': all eigenvalues will be found;
63 *> = 'V': all eigenvalues in the half-open interval (VL,VU]
64 *> will be found;
65 *> = 'I': the IL-th through IU-th eigenvalues will be found.
66 *> \endverbatim
67 *>
68 *> \param[in] UPLO
69 *> \verbatim
70 *> UPLO is CHARACTER*1
71 *> = 'U': Upper triangle of A is stored;
72 *> = 'L': Lower triangle of A is stored.
73 *> \endverbatim
74 *>
75 *> \param[in] N
76 *> \verbatim
77 *> N is INTEGER
78 *> The order of the matrix A. N >= 0.
79 *> \endverbatim
80 *>
81 *> \param[in] KD
82 *> \verbatim
83 *> KD is INTEGER
84 *> The number of superdiagonals of the matrix A if UPLO = 'U',
85 *> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
86 *> \endverbatim
87 *>
88 *> \param[in,out] AB
89 *> \verbatim
90 *> AB is COMPLEX*16 array, dimension (LDAB, N)
91 *> On entry, the upper or lower triangle of the Hermitian band
92 *> matrix A, stored in the first KD+1 rows of the array. The
93 *> j-th column of A is stored in the j-th column of the array AB
94 *> as follows:
95 *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
96 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
97 *>
98 *> On exit, AB is overwritten by values generated during the
99 *> reduction to tridiagonal form.
100 *> \endverbatim
101 *>
102 *> \param[in] LDAB
103 *> \verbatim
104 *> LDAB is INTEGER
105 *> The leading dimension of the array AB. LDAB >= KD + 1.
106 *> \endverbatim
107 *>
108 *> \param[out] Q
109 *> \verbatim
110 *> Q is COMPLEX*16 array, dimension (LDQ, N)
111 *> If JOBZ = 'V', the N-by-N unitary matrix used in the
112 *> reduction to tridiagonal form.
113 *> If JOBZ = 'N', the array Q is not referenced.
114 *> \endverbatim
115 *>
116 *> \param[in] LDQ
117 *> \verbatim
118 *> LDQ is INTEGER
119 *> The leading dimension of the array Q. If JOBZ = 'V', then
120 *> LDQ >= max(1,N).
121 *> \endverbatim
122 *>
123 *> \param[in] VL
124 *> \verbatim
125 *> VL is DOUBLE PRECISION
126 *> If RANGE='V', the lower bound of the interval to
127 *> be searched for eigenvalues. VL < VU.
128 *> Not referenced if RANGE = 'A' or 'I'.
129 *> \endverbatim
130 *>
131 *> \param[in] VU
132 *> \verbatim
133 *> VU is DOUBLE PRECISION
134 *> If RANGE='V', the upper bound of the interval to
135 *> be searched for eigenvalues. VL < VU.
136 *> Not referenced if RANGE = 'A' or 'I'.
137 *> \endverbatim
138 *>
139 *> \param[in] IL
140 *> \verbatim
141 *> IL is INTEGER
142 *> If RANGE='I', the index of the
143 *> smallest eigenvalue to be returned.
144 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
145 *> Not referenced if RANGE = 'A' or 'V'.
146 *> \endverbatim
147 *>
148 *> \param[in] IU
149 *> \verbatim
150 *> IU is INTEGER
151 *> If RANGE='I', the index of the
152 *> largest eigenvalue to be returned.
153 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
154 *> Not referenced if RANGE = 'A' or 'V'.
155 *> \endverbatim
156 *>
157 *> \param[in] ABSTOL
158 *> \verbatim
159 *> ABSTOL is DOUBLE PRECISION
160 *> The absolute error tolerance for the eigenvalues.
161 *> An approximate eigenvalue is accepted as converged
162 *> when it is determined to lie in an interval [a,b]
163 *> of width less than or equal to
164 *>
165 *> ABSTOL + EPS * max( |a|,|b| ) ,
166 *>
167 *> where EPS is the machine precision. If ABSTOL is less than
168 *> or equal to zero, then EPS*|T| will be used in its place,
169 *> where |T| is the 1-norm of the tridiagonal matrix obtained
170 *> by reducing AB to tridiagonal form.
171 *>
172 *> Eigenvalues will be computed most accurately when ABSTOL is
173 *> set to twice the underflow threshold 2*DLAMCH('S'), not zero.
174 *> If this routine returns with INFO>0, indicating that some
175 *> eigenvectors did not converge, try setting ABSTOL to
176 *> 2*DLAMCH('S').
177 *>
178 *> See "Computing Small Singular Values of Bidiagonal Matrices
179 *> with Guaranteed High Relative Accuracy," by Demmel and
180 *> Kahan, LAPACK Working Note #3.
181 *> \endverbatim
182 *>
183 *> \param[out] M
184 *> \verbatim
185 *> M is INTEGER
186 *> The total number of eigenvalues found. 0 <= M <= N.
187 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
188 *> \endverbatim
189 *>
190 *> \param[out] W
191 *> \verbatim
192 *> W is DOUBLE PRECISION array, dimension (N)
193 *> The first M elements contain the selected eigenvalues in
194 *> ascending order.
195 *> \endverbatim
196 *>
197 *> \param[out] Z
198 *> \verbatim
199 *> Z is COMPLEX*16 array, dimension (LDZ, max(1,M))
200 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
201 *> contain the orthonormal eigenvectors of the matrix A
202 *> corresponding to the selected eigenvalues, with the i-th
203 *> column of Z holding the eigenvector associated with W(i).
204 *> If an eigenvector fails to converge, then that column of Z
205 *> contains the latest approximation to the eigenvector, and the
206 *> index of the eigenvector is returned in IFAIL.
207 *> If JOBZ = 'N', then Z is not referenced.
208 *> Note: the user must ensure that at least max(1,M) columns are
209 *> supplied in the array Z; if RANGE = 'V', the exact value of M
210 *> is not known in advance and an upper bound must be used.
211 *> \endverbatim
212 *>
213 *> \param[in] LDZ
214 *> \verbatim
215 *> LDZ is INTEGER
216 *> The leading dimension of the array Z. LDZ >= 1, and if
217 *> JOBZ = 'V', LDZ >= max(1,N).
218 *> \endverbatim
219 *>
220 *> \param[out] WORK
221 *> \verbatim
222 *> WORK is COMPLEX*16 array, dimension (N)
223 *> \endverbatim
224 *>
225 *> \param[out] RWORK
226 *> \verbatim
227 *> RWORK is DOUBLE PRECISION array, dimension (7*N)
228 *> \endverbatim
229 *>
230 *> \param[out] IWORK
231 *> \verbatim
232 *> IWORK is INTEGER array, dimension (5*N)
233 *> \endverbatim
234 *>
235 *> \param[out] IFAIL
236 *> \verbatim
237 *> IFAIL is INTEGER array, dimension (N)
238 *> If JOBZ = 'V', then if INFO = 0, the first M elements of
239 *> IFAIL are zero. If INFO > 0, then IFAIL contains the
240 *> indices of the eigenvectors that failed to converge.
241 *> If JOBZ = 'N', then IFAIL is not referenced.
242 *> \endverbatim
243 *>
244 *> \param[out] INFO
245 *> \verbatim
246 *> INFO is INTEGER
247 *> = 0: successful exit
248 *> < 0: if INFO = -i, the i-th argument had an illegal value
249 *> > 0: if INFO = i, then i eigenvectors failed to converge.
250 *> Their indices are stored in array IFAIL.
251 *> \endverbatim
252 *
253 * Authors:
254 * ========
255 *
256 *> \author Univ. of Tennessee
257 *> \author Univ. of California Berkeley
258 *> \author Univ. of Colorado Denver
259 *> \author NAG Ltd.
260 *
261 *> \date June 2016
262 *
263 *> \ingroup complex16OTHEReigen
264 *
265 * =====================================================================
266  SUBROUTINE zhbevx( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ, VL,
267  $ vu, il, iu, abstol, m, w, z, ldz, work, rwork,
268  $ iwork, ifail, info )
269 *
270 * -- LAPACK driver routine (version 3.6.1) --
271 * -- LAPACK is a software package provided by Univ. of Tennessee, --
272 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
273 * June 2016
274 *
275 * .. Scalar Arguments ..
276  CHARACTER JOBZ, RANGE, UPLO
277  INTEGER IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N
278  DOUBLE PRECISION ABSTOL, VL, VU
279 * ..
280 * .. Array Arguments ..
281  INTEGER IFAIL( * ), IWORK( * )
282  DOUBLE PRECISION RWORK( * ), W( * )
283  COMPLEX*16 AB( ldab, * ), Q( ldq, * ), WORK( * ),
284  $ z( ldz, * )
285 * ..
286 *
287 * =====================================================================
288 *
289 * .. Parameters ..
290  DOUBLE PRECISION ZERO, ONE
291  parameter ( zero = 0.0d0, one = 1.0d0 )
292  COMPLEX*16 CZERO, CONE
293  parameter ( czero = ( 0.0d0, 0.0d0 ),
294  $ cone = ( 1.0d0, 0.0d0 ) )
295 * ..
296 * .. Local Scalars ..
297  LOGICAL ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ
298  CHARACTER ORDER
299  INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
300  $ indisp, indiwk, indrwk, indwrk, iscale, itmp1,
301  $ j, jj, nsplit
302  DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
303  $ sigma, smlnum, tmp1, vll, vuu
304  COMPLEX*16 CTMP1
305 * ..
306 * .. External Functions ..
307  LOGICAL LSAME
308  DOUBLE PRECISION DLAMCH, ZLANHB
309  EXTERNAL lsame, dlamch, zlanhb
310 * ..
311 * .. External Subroutines ..
312  EXTERNAL dcopy, dscal, dstebz, dsterf, xerbla, zcopy,
314  $ zswap
315 * ..
316 * .. Intrinsic Functions ..
317  INTRINSIC dble, max, min, sqrt
318 * ..
319 * .. Executable Statements ..
320 *
321 * Test the input parameters.
322 *
323  wantz = lsame( jobz, 'V' )
324  alleig = lsame( range, 'A' )
325  valeig = lsame( range, 'V' )
326  indeig = lsame( range, 'I' )
327  lower = lsame( uplo, 'L' )
328 *
329  info = 0
330  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
331  info = -1
332  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
333  info = -2
334  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
335  info = -3
336  ELSE IF( n.LT.0 ) THEN
337  info = -4
338  ELSE IF( kd.LT.0 ) THEN
339  info = -5
340  ELSE IF( ldab.LT.kd+1 ) THEN
341  info = -7
342  ELSE IF( wantz .AND. ldq.LT.max( 1, n ) ) THEN
343  info = -9
344  ELSE
345  IF( valeig ) THEN
346  IF( n.GT.0 .AND. vu.LE.vl )
347  $ info = -11
348  ELSE IF( indeig ) THEN
349  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
350  info = -12
351  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
352  info = -13
353  END IF
354  END IF
355  END IF
356  IF( info.EQ.0 ) THEN
357  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
358  $ info = -18
359  END IF
360 *
361  IF( info.NE.0 ) THEN
362  CALL xerbla( 'ZHBEVX', -info )
363  RETURN
364  END IF
365 *
366 * Quick return if possible
367 *
368  m = 0
369  IF( n.EQ.0 )
370  $ RETURN
371 *
372  IF( n.EQ.1 ) THEN
373  m = 1
374  IF( lower ) THEN
375  ctmp1 = ab( 1, 1 )
376  ELSE
377  ctmp1 = ab( kd+1, 1 )
378  END IF
379  tmp1 = dble( ctmp1 )
380  IF( valeig ) THEN
381  IF( .NOT.( vl.LT.tmp1 .AND. vu.GE.tmp1 ) )
382  $ m = 0
383  END IF
384  IF( m.EQ.1 ) THEN
385  w( 1 ) = ctmp1
386  IF( wantz )
387  $ z( 1, 1 ) = cone
388  END IF
389  RETURN
390  END IF
391 *
392 * Get machine constants.
393 *
394  safmin = dlamch( 'Safe minimum' )
395  eps = dlamch( 'Precision' )
396  smlnum = safmin / eps
397  bignum = one / smlnum
398  rmin = sqrt( smlnum )
399  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
400 *
401 * Scale matrix to allowable range, if necessary.
402 *
403  iscale = 0
404  abstll = abstol
405  IF( valeig ) THEN
406  vll = vl
407  vuu = vu
408  ELSE
409  vll = zero
410  vuu = zero
411  END IF
412  anrm = zlanhb( 'M', uplo, n, kd, ab, ldab, rwork )
413  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
414  iscale = 1
415  sigma = rmin / anrm
416  ELSE IF( anrm.GT.rmax ) THEN
417  iscale = 1
418  sigma = rmax / anrm
419  END IF
420  IF( iscale.EQ.1 ) THEN
421  IF( lower ) THEN
422  CALL zlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
423  ELSE
424  CALL zlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
425  END IF
426  IF( abstol.GT.0 )
427  $ abstll = abstol*sigma
428  IF( valeig ) THEN
429  vll = vl*sigma
430  vuu = vu*sigma
431  END IF
432  END IF
433 *
434 * Call ZHBTRD to reduce Hermitian band matrix to tridiagonal form.
435 *
436  indd = 1
437  inde = indd + n
438  indrwk = inde + n
439  indwrk = 1
440  CALL zhbtrd( jobz, uplo, n, kd, ab, ldab, rwork( indd ),
441  $ rwork( inde ), q, ldq, work( indwrk ), iinfo )
442 *
443 * If all eigenvalues are desired and ABSTOL is less than or equal
444 * to zero, then call DSTERF or ZSTEQR. If this fails for some
445 * eigenvalue, then try DSTEBZ.
446 *
447  test = .false.
448  IF (indeig) THEN
449  IF (il.EQ.1 .AND. iu.EQ.n) THEN
450  test = .true.
451  END IF
452  END IF
453  IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
454  CALL dcopy( n, rwork( indd ), 1, w, 1 )
455  indee = indrwk + 2*n
456  IF( .NOT.wantz ) THEN
457  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
458  CALL dsterf( n, w, rwork( indee ), info )
459  ELSE
460  CALL zlacpy( 'A', n, n, q, ldq, z, ldz )
461  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
462  CALL zsteqr( jobz, n, w, rwork( indee ), z, ldz,
463  $ rwork( indrwk ), info )
464  IF( info.EQ.0 ) THEN
465  DO 10 i = 1, n
466  ifail( i ) = 0
467  10 CONTINUE
468  END IF
469  END IF
470  IF( info.EQ.0 ) THEN
471  m = n
472  GO TO 30
473  END IF
474  info = 0
475  END IF
476 *
477 * Otherwise, call DSTEBZ and, if eigenvectors are desired, ZSTEIN.
478 *
479  IF( wantz ) THEN
480  order = 'B'
481  ELSE
482  order = 'E'
483  END IF
484  indibl = 1
485  indisp = indibl + n
486  indiwk = indisp + n
487  CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
488  $ rwork( indd ), rwork( inde ), m, nsplit, w,
489  $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
490  $ iwork( indiwk ), info )
491 *
492  IF( wantz ) THEN
493  CALL zstein( n, rwork( indd ), rwork( inde ), m, w,
494  $ iwork( indibl ), iwork( indisp ), z, ldz,
495  $ rwork( indrwk ), iwork( indiwk ), ifail, info )
496 *
497 * Apply unitary matrix used in reduction to tridiagonal
498 * form to eigenvectors returned by ZSTEIN.
499 *
500  DO 20 j = 1, m
501  CALL zcopy( n, z( 1, j ), 1, work( 1 ), 1 )
502  CALL zgemv( 'N', n, n, cone, q, ldq, work, 1, czero,
503  $ z( 1, j ), 1 )
504  20 CONTINUE
505  END IF
506 *
507 * If matrix was scaled, then rescale eigenvalues appropriately.
508 *
509  30 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 dscal( 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 50 j = 1, m - 1
524  i = 0
525  tmp1 = w( j )
526  DO 40 jj = j + 1, m
527  IF( w( jj ).LT.tmp1 ) THEN
528  i = jj
529  tmp1 = w( jj )
530  END IF
531  40 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 zswap( 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  50 CONTINUE
547  END IF
548 *
549  RETURN
550 *
551 * End of ZHBEVX
552 *
553  END
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:88
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:275
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zhbevx(JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, RWORK, IWORK, IFAIL, INFO)
ZHBEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matric...
Definition: zhbevx.f:269
subroutine zsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
ZSTEQR
Definition: zsteqr.f:134
subroutine zstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
ZSTEIN
Definition: zstein.f:184
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:145
subroutine zhbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
ZHBTRD
Definition: zhbtrd.f:165