LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \endverbatim
127 *>
128 *> \param[in] VU
129 *> \verbatim
130 *> VU is DOUBLE PRECISION
131 *> If RANGE='V', the lower and upper bounds of the interval to
132 *> be searched for eigenvalues. VL < VU.
133 *> Not referenced if RANGE = 'A' or 'I'.
134 *> \endverbatim
135 *>
136 *> \param[in] IL
137 *> \verbatim
138 *> IL is INTEGER
139 *> \endverbatim
140 *>
141 *> \param[in] IU
142 *> \verbatim
143 *> IU is INTEGER
144 *> If RANGE='I', the indices (in ascending order) of the
145 *> smallest and largest eigenvalues to be returned.
146 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
147 *> Not referenced if RANGE = 'A' or 'V'.
148 *> \endverbatim
149 *>
150 *> \param[in] ABSTOL
151 *> \verbatim
152 *> ABSTOL is DOUBLE PRECISION
153 *> The absolute error tolerance for the eigenvalues.
154 *> An approximate eigenvalue is accepted as converged
155 *> when it is determined to lie in an interval [a,b]
156 *> of width less than or equal to
157 *>
158 *> ABSTOL + EPS * max( |a|,|b| ) ,
159 *>
160 *> where EPS is the machine precision. If ABSTOL is less than
161 *> or equal to zero, then EPS*|T| will be used in its place,
162 *> where |T| is the 1-norm of the tridiagonal matrix obtained
163 *> by reducing AB to tridiagonal form.
164 *>
165 *> Eigenvalues will be computed most accurately when ABSTOL is
166 *> set to twice the underflow threshold 2*DLAMCH('S'), not zero.
167 *> If this routine returns with INFO>0, indicating that some
168 *> eigenvectors did not converge, try setting ABSTOL to
169 *> 2*DLAMCH('S').
170 *>
171 *> See "Computing Small Singular Values of Bidiagonal Matrices
172 *> with Guaranteed High Relative Accuracy," by Demmel and
173 *> Kahan, LAPACK Working Note #3.
174 *> \endverbatim
175 *>
176 *> \param[out] M
177 *> \verbatim
178 *> M is INTEGER
179 *> The total number of eigenvalues found. 0 <= M <= N.
180 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
181 *> \endverbatim
182 *>
183 *> \param[out] W
184 *> \verbatim
185 *> W is DOUBLE PRECISION array, dimension (N)
186 *> The first M elements contain the selected eigenvalues in
187 *> ascending order.
188 *> \endverbatim
189 *>
190 *> \param[out] Z
191 *> \verbatim
192 *> Z is COMPLEX*16 array, dimension (LDZ, max(1,M))
193 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
194 *> contain the orthonormal eigenvectors of the matrix A
195 *> corresponding to the selected eigenvalues, with the i-th
196 *> column of Z holding the eigenvector associated with W(i).
197 *> If an eigenvector fails to converge, then that column of Z
198 *> contains the latest approximation to the eigenvector, and the
199 *> index of the eigenvector is returned in IFAIL.
200 *> If JOBZ = 'N', then Z is not referenced.
201 *> Note: the user must ensure that at least max(1,M) columns are
202 *> supplied in the array Z; if RANGE = 'V', the exact value of M
203 *> is not known in advance and an upper bound must be used.
204 *> \endverbatim
205 *>
206 *> \param[in] LDZ
207 *> \verbatim
208 *> LDZ is INTEGER
209 *> The leading dimension of the array Z. LDZ >= 1, and if
210 *> JOBZ = 'V', LDZ >= max(1,N).
211 *> \endverbatim
212 *>
213 *> \param[out] WORK
214 *> \verbatim
215 *> WORK is COMPLEX*16 array, dimension (N)
216 *> \endverbatim
217 *>
218 *> \param[out] RWORK
219 *> \verbatim
220 *> RWORK is DOUBLE PRECISION array, dimension (7*N)
221 *> \endverbatim
222 *>
223 *> \param[out] IWORK
224 *> \verbatim
225 *> IWORK is INTEGER array, dimension (5*N)
226 *> \endverbatim
227 *>
228 *> \param[out] IFAIL
229 *> \verbatim
230 *> IFAIL is INTEGER array, dimension (N)
231 *> If JOBZ = 'V', then if INFO = 0, the first M elements of
232 *> IFAIL are zero. If INFO > 0, then IFAIL contains the
233 *> indices of the eigenvectors that failed to converge.
234 *> If JOBZ = 'N', then IFAIL is not referenced.
235 *> \endverbatim
236 *>
237 *> \param[out] INFO
238 *> \verbatim
239 *> INFO is INTEGER
240 *> = 0: successful exit
241 *> < 0: if INFO = -i, the i-th argument had an illegal value
242 *> > 0: if INFO = i, then i eigenvectors failed to converge.
243 *> Their indices are stored in array IFAIL.
244 *> \endverbatim
245 *
246 * Authors:
247 * ========
248 *
249 *> \author Univ. of Tennessee
250 *> \author Univ. of California Berkeley
251 *> \author Univ. of Colorado Denver
252 *> \author NAG Ltd.
253 *
254 *> \date November 2011
255 *
256 *> \ingroup complex16OTHEReigen
257 *
258 * =====================================================================
259  SUBROUTINE zhbevx( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ, VL,
260  $ vu, il, iu, abstol, m, w, z, ldz, work, rwork,
261  $ iwork, ifail, info )
262 *
263 * -- LAPACK driver routine (version 3.4.0) --
264 * -- LAPACK is a software package provided by Univ. of Tennessee, --
265 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
266 * November 2011
267 *
268 * .. Scalar Arguments ..
269  CHARACTER jobz, range, uplo
270  INTEGER il, info, iu, kd, ldab, ldq, ldz, m, n
271  DOUBLE PRECISION abstol, vl, vu
272 * ..
273 * .. Array Arguments ..
274  INTEGER ifail( * ), iwork( * )
275  DOUBLE PRECISION rwork( * ), w( * )
276  COMPLEX*16 ab( ldab, * ), q( ldq, * ), work( * ),
277  $ z( ldz, * )
278 * ..
279 *
280 * =====================================================================
281 *
282 * .. Parameters ..
283  DOUBLE PRECISION zero, one
284  parameter( zero = 0.0d0, one = 1.0d0 )
285  COMPLEX*16 czero, cone
286  parameter( czero = ( 0.0d0, 0.0d0 ),
287  $ cone = ( 1.0d0, 0.0d0 ) )
288 * ..
289 * .. Local Scalars ..
290  LOGICAL alleig, indeig, lower, test, valeig, wantz
291  CHARACTER order
292  INTEGER i, iinfo, imax, indd, inde, indee, indibl,
293  $ indisp, indiwk, indrwk, indwrk, iscale, itmp1,
294  $ j, jj, nsplit
295  DOUBLE PRECISION abstll, anrm, bignum, eps, rmax, rmin, safmin,
296  $ sigma, smlnum, tmp1, vll, vuu
297  COMPLEX*16 ctmp1
298 * ..
299 * .. External Functions ..
300  LOGICAL lsame
301  DOUBLE PRECISION dlamch, zlanhb
302  EXTERNAL lsame, dlamch, zlanhb
303 * ..
304 * .. External Subroutines ..
305  EXTERNAL dcopy, dscal, dstebz, dsterf, xerbla, zcopy,
307  $ zswap
308 * ..
309 * .. Intrinsic Functions ..
310  INTRINSIC dble, max, min, sqrt
311 * ..
312 * .. Executable Statements ..
313 *
314 * Test the input parameters.
315 *
316  wantz = lsame( jobz, 'V' )
317  alleig = lsame( range, 'A' )
318  valeig = lsame( range, 'V' )
319  indeig = lsame( range, 'I' )
320  lower = lsame( uplo, 'L' )
321 *
322  info = 0
323  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
324  info = -1
325  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
326  info = -2
327  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
328  info = -3
329  ELSE IF( n.LT.0 ) THEN
330  info = -4
331  ELSE IF( kd.LT.0 ) THEN
332  info = -5
333  ELSE IF( ldab.LT.kd+1 ) THEN
334  info = -7
335  ELSE IF( wantz .AND. ldq.LT.max( 1, n ) ) THEN
336  info = -9
337  ELSE
338  IF( valeig ) THEN
339  IF( n.GT.0 .AND. vu.LE.vl )
340  $ info = -11
341  ELSE IF( indeig ) THEN
342  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
343  info = -12
344  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
345  info = -13
346  END IF
347  END IF
348  END IF
349  IF( info.EQ.0 ) THEN
350  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
351  $ info = -18
352  END IF
353 *
354  IF( info.NE.0 ) THEN
355  CALL xerbla( 'ZHBEVX', -info )
356  return
357  END IF
358 *
359 * Quick return if possible
360 *
361  m = 0
362  IF( n.EQ.0 )
363  $ return
364 *
365  IF( n.EQ.1 ) THEN
366  m = 1
367  IF( lower ) THEN
368  ctmp1 = ab( 1, 1 )
369  ELSE
370  ctmp1 = ab( kd+1, 1 )
371  END IF
372  tmp1 = dble( ctmp1 )
373  IF( valeig ) THEN
374  IF( .NOT.( vl.LT.tmp1 .AND. vu.GE.tmp1 ) )
375  $ m = 0
376  END IF
377  IF( m.EQ.1 ) THEN
378  w( 1 ) = ctmp1
379  IF( wantz )
380  $ z( 1, 1 ) = cone
381  END IF
382  return
383  END IF
384 *
385 * Get machine constants.
386 *
387  safmin = dlamch( 'Safe minimum' )
388  eps = dlamch( 'Precision' )
389  smlnum = safmin / eps
390  bignum = one / smlnum
391  rmin = sqrt( smlnum )
392  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
393 *
394 * Scale matrix to allowable range, if necessary.
395 *
396  iscale = 0
397  abstll = abstol
398  IF( valeig ) THEN
399  vll = vl
400  vuu = vu
401  ELSE
402  vll = zero
403  vuu = zero
404  END IF
405  anrm = zlanhb( 'M', uplo, n, kd, ab, ldab, rwork )
406  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
407  iscale = 1
408  sigma = rmin / anrm
409  ELSE IF( anrm.GT.rmax ) THEN
410  iscale = 1
411  sigma = rmax / anrm
412  END IF
413  IF( iscale.EQ.1 ) THEN
414  IF( lower ) THEN
415  CALL zlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
416  ELSE
417  CALL zlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
418  END IF
419  IF( abstol.GT.0 )
420  $ abstll = abstol*sigma
421  IF( valeig ) THEN
422  vll = vl*sigma
423  vuu = vu*sigma
424  END IF
425  END IF
426 *
427 * Call ZHBTRD to reduce Hermitian band matrix to tridiagonal form.
428 *
429  indd = 1
430  inde = indd + n
431  indrwk = inde + n
432  indwrk = 1
433  CALL zhbtrd( jobz, uplo, n, kd, ab, ldab, rwork( indd ),
434  $ rwork( inde ), q, ldq, work( indwrk ), iinfo )
435 *
436 * If all eigenvalues are desired and ABSTOL is less than or equal
437 * to zero, then call DSTERF or ZSTEQR. If this fails for some
438 * eigenvalue, then try DSTEBZ.
439 *
440  test = .false.
441  IF (indeig) THEN
442  IF (il.EQ.1 .AND. iu.EQ.n) THEN
443  test = .true.
444  END IF
445  END IF
446  IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
447  CALL dcopy( n, rwork( indd ), 1, w, 1 )
448  indee = indrwk + 2*n
449  IF( .NOT.wantz ) THEN
450  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
451  CALL dsterf( n, w, rwork( indee ), info )
452  ELSE
453  CALL zlacpy( 'A', n, n, q, ldq, z, ldz )
454  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
455  CALL zsteqr( jobz, n, w, rwork( indee ), z, ldz,
456  $ rwork( indrwk ), info )
457  IF( info.EQ.0 ) THEN
458  DO 10 i = 1, n
459  ifail( i ) = 0
460  10 continue
461  END IF
462  END IF
463  IF( info.EQ.0 ) THEN
464  m = n
465  go to 30
466  END IF
467  info = 0
468  END IF
469 *
470 * Otherwise, call DSTEBZ and, if eigenvectors are desired, ZSTEIN.
471 *
472  IF( wantz ) THEN
473  order = 'B'
474  ELSE
475  order = 'E'
476  END IF
477  indibl = 1
478  indisp = indibl + n
479  indiwk = indisp + n
480  CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
481  $ rwork( indd ), rwork( inde ), m, nsplit, w,
482  $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
483  $ iwork( indiwk ), info )
484 *
485  IF( wantz ) THEN
486  CALL zstein( n, rwork( indd ), rwork( inde ), m, w,
487  $ iwork( indibl ), iwork( indisp ), z, ldz,
488  $ rwork( indrwk ), iwork( indiwk ), ifail, info )
489 *
490 * Apply unitary matrix used in reduction to tridiagonal
491 * form to eigenvectors returned by ZSTEIN.
492 *
493  DO 20 j = 1, m
494  CALL zcopy( n, z( 1, j ), 1, work( 1 ), 1 )
495  CALL zgemv( 'N', n, n, cone, q, ldq, work, 1, czero,
496  $ z( 1, j ), 1 )
497  20 continue
498  END IF
499 *
500 * If matrix was scaled, then rescale eigenvalues appropriately.
501 *
502  30 continue
503  IF( iscale.EQ.1 ) THEN
504  IF( info.EQ.0 ) THEN
505  imax = m
506  ELSE
507  imax = info - 1
508  END IF
509  CALL dscal( imax, one / sigma, w, 1 )
510  END IF
511 *
512 * If eigenvalues are not in order, then sort them, along with
513 * eigenvectors.
514 *
515  IF( wantz ) THEN
516  DO 50 j = 1, m - 1
517  i = 0
518  tmp1 = w( j )
519  DO 40 jj = j + 1, m
520  IF( w( jj ).LT.tmp1 ) THEN
521  i = jj
522  tmp1 = w( jj )
523  END IF
524  40 continue
525 *
526  IF( i.NE.0 ) THEN
527  itmp1 = iwork( indibl+i-1 )
528  w( i ) = w( j )
529  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
530  w( j ) = tmp1
531  iwork( indibl+j-1 ) = itmp1
532  CALL zswap( n, z( 1, i ), 1, z( 1, j ), 1 )
533  IF( info.NE.0 ) THEN
534  itmp1 = ifail( i )
535  ifail( i ) = ifail( j )
536  ifail( j ) = itmp1
537  END IF
538  END IF
539  50 continue
540  END IF
541 *
542  return
543 *
544 * End of ZHBEVX
545 *
546  END