LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
chbevx.f
Go to the documentation of this file.
1 *> \brief <b> CHBEVX 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 CHBEVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chbevx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chbevx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chbevx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CHBEVX( 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 * REAL ABSTOL, VL, VU
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IFAIL( * ), IWORK( * )
32 * REAL RWORK( * ), W( * )
33 * COMPLEX AB( LDAB, * ), Q( LDQ, * ), WORK( * ),
34 * $ Z( LDZ, * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> CHBEVX 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 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 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 REAL
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 REAL
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 REAL
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*SLAMCH('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*SLAMCH('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 REAL 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 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 array, dimension (N)
223 *> \endverbatim
224 *>
225 *> \param[out] RWORK
226 *> \verbatim
227 *> RWORK is REAL 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 *> \ingroup complexOTHEReigen
262 *
263 * =====================================================================
264  SUBROUTINE chbevx( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ, VL,
265  $ VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, RWORK,
266  $ IWORK, IFAIL, INFO )
267 *
268 * -- LAPACK driver routine --
269 * -- LAPACK is a software package provided by Univ. of Tennessee, --
270 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
271 *
272 * .. Scalar Arguments ..
273  CHARACTER JOBZ, RANGE, UPLO
274  INTEGER IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N
275  REAL ABSTOL, VL, VU
276 * ..
277 * .. Array Arguments ..
278  INTEGER IFAIL( * ), IWORK( * )
279  REAL RWORK( * ), W( * )
280  COMPLEX AB( LDAB, * ), Q( LDQ, * ), WORK( * ),
281  $ z( ldz, * )
282 * ..
283 *
284 * =====================================================================
285 *
286 * .. Parameters ..
287  REAL ZERO, ONE
288  PARAMETER ( ZERO = 0.0e0, one = 1.0e0 )
289  COMPLEX CZERO, CONE
290  parameter( czero = ( 0.0e0, 0.0e0 ),
291  $ cone = ( 1.0e0, 0.0e0 ) )
292 * ..
293 * .. Local Scalars ..
294  LOGICAL ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ
295  CHARACTER ORDER
296  INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
297  $ indisp, indiwk, indrwk, indwrk, iscale, itmp1,
298  $ j, jj, nsplit
299  REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
300  $ SIGMA, SMLNUM, TMP1, VLL, VUU
301  COMPLEX CTMP1
302 * ..
303 * .. External Functions ..
304  LOGICAL LSAME
305  REAL CLANHB, SLAMCH
306  EXTERNAL lsame, clanhb, slamch
307 * ..
308 * .. External Subroutines ..
309  EXTERNAL ccopy, cgemv, chbtrd, clacpy, clascl, cstein,
311  $ xerbla
312 * ..
313 * .. Intrinsic Functions ..
314  INTRINSIC max, min, real, sqrt
315 * ..
316 * .. Executable Statements ..
317 *
318 * Test the input parameters.
319 *
320  wantz = lsame( jobz, 'V' )
321  alleig = lsame( range, 'A' )
322  valeig = lsame( range, 'V' )
323  indeig = lsame( range, 'I' )
324  lower = lsame( uplo, 'L' )
325 *
326  info = 0
327  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
328  info = -1
329  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
330  info = -2
331  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
332  info = -3
333  ELSE IF( n.LT.0 ) THEN
334  info = -4
335  ELSE IF( kd.LT.0 ) THEN
336  info = -5
337  ELSE IF( ldab.LT.kd+1 ) THEN
338  info = -7
339  ELSE IF( wantz .AND. ldq.LT.max( 1, n ) ) THEN
340  info = -9
341  ELSE
342  IF( valeig ) THEN
343  IF( n.GT.0 .AND. vu.LE.vl )
344  $ info = -11
345  ELSE IF( indeig ) THEN
346  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
347  info = -12
348  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
349  info = -13
350  END IF
351  END IF
352  END IF
353  IF( info.EQ.0 ) THEN
354  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
355  $ info = -18
356  END IF
357 *
358  IF( info.NE.0 ) THEN
359  CALL xerbla( 'CHBEVX', -info )
360  RETURN
361  END IF
362 *
363 * Quick return if possible
364 *
365  m = 0
366  IF( n.EQ.0 )
367  $ RETURN
368 *
369  IF( n.EQ.1 ) THEN
370  m = 1
371  IF( lower ) THEN
372  ctmp1 = ab( 1, 1 )
373  ELSE
374  ctmp1 = ab( kd+1, 1 )
375  END IF
376  tmp1 = real( ctmp1 )
377  IF( valeig ) THEN
378  IF( .NOT.( vl.LT.tmp1 .AND. vu.GE.tmp1 ) )
379  $ m = 0
380  END IF
381  IF( m.EQ.1 ) THEN
382  w( 1 ) = real( ctmp1 )
383  IF( wantz )
384  $ z( 1, 1 ) = cone
385  END IF
386  RETURN
387  END IF
388 *
389 * Get machine constants.
390 *
391  safmin = slamch( 'Safe minimum' )
392  eps = slamch( 'Precision' )
393  smlnum = safmin / eps
394  bignum = one / smlnum
395  rmin = sqrt( smlnum )
396  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
397 *
398 * Scale matrix to allowable range, if necessary.
399 *
400  iscale = 0
401  abstll = abstol
402  IF ( valeig ) THEN
403  vll = vl
404  vuu = vu
405  ELSE
406  vll = zero
407  vuu = zero
408  ENDIF
409  anrm = clanhb( 'M', uplo, n, kd, ab, ldab, rwork )
410  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
411  iscale = 1
412  sigma = rmin / anrm
413  ELSE IF( anrm.GT.rmax ) THEN
414  iscale = 1
415  sigma = rmax / anrm
416  END IF
417  IF( iscale.EQ.1 ) THEN
418  IF( lower ) THEN
419  CALL clascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
420  ELSE
421  CALL clascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
422  END IF
423  IF( abstol.GT.0 )
424  $ abstll = abstol*sigma
425  IF( valeig ) THEN
426  vll = vl*sigma
427  vuu = vu*sigma
428  END IF
429  END IF
430 *
431 * Call CHBTRD to reduce Hermitian band matrix to tridiagonal form.
432 *
433  indd = 1
434  inde = indd + n
435  indrwk = inde + n
436  indwrk = 1
437  CALL chbtrd( jobz, uplo, n, kd, ab, ldab, rwork( indd ),
438  $ rwork( inde ), q, ldq, work( indwrk ), iinfo )
439 *
440 * If all eigenvalues are desired and ABSTOL is less than or equal
441 * to zero, then call SSTERF or CSTEQR. If this fails for some
442 * eigenvalue, then try SSTEBZ.
443 *
444  test = .false.
445  IF (indeig) THEN
446  IF (il.EQ.1 .AND. iu.EQ.n) THEN
447  test = .true.
448  END IF
449  END IF
450  IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
451  CALL scopy( n, rwork( indd ), 1, w, 1 )
452  indee = indrwk + 2*n
453  IF( .NOT.wantz ) THEN
454  CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
455  CALL ssterf( n, w, rwork( indee ), info )
456  ELSE
457  CALL clacpy( 'A', n, n, q, ldq, z, ldz )
458  CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
459  CALL csteqr( jobz, n, w, rwork( indee ), z, ldz,
460  $ rwork( indrwk ), info )
461  IF( info.EQ.0 ) THEN
462  DO 10 i = 1, n
463  ifail( i ) = 0
464  10 CONTINUE
465  END IF
466  END IF
467  IF( info.EQ.0 ) THEN
468  m = n
469  GO TO 30
470  END IF
471  info = 0
472  END IF
473 *
474 * Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN.
475 *
476  IF( wantz ) THEN
477  order = 'B'
478  ELSE
479  order = 'E'
480  END IF
481  indibl = 1
482  indisp = indibl + n
483  indiwk = indisp + n
484  CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
485  $ rwork( indd ), rwork( inde ), m, nsplit, w,
486  $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
487  $ iwork( indiwk ), info )
488 *
489  IF( wantz ) THEN
490  CALL cstein( n, rwork( indd ), rwork( inde ), m, w,
491  $ iwork( indibl ), iwork( indisp ), z, ldz,
492  $ rwork( indrwk ), iwork( indiwk ), ifail, info )
493 *
494 * Apply unitary matrix used in reduction to tridiagonal
495 * form to eigenvectors returned by CSTEIN.
496 *
497  DO 20 j = 1, m
498  CALL ccopy( n, z( 1, j ), 1, work( 1 ), 1 )
499  CALL cgemv( 'N', n, n, cone, q, ldq, work, 1, czero,
500  $ z( 1, j ), 1 )
501  20 CONTINUE
502  END IF
503 *
504 * If matrix was scaled, then rescale eigenvalues appropriately.
505 *
506  30 CONTINUE
507  IF( iscale.EQ.1 ) THEN
508  IF( info.EQ.0 ) THEN
509  imax = m
510  ELSE
511  imax = info - 1
512  END IF
513  CALL sscal( imax, one / sigma, w, 1 )
514  END IF
515 *
516 * If eigenvalues are not in order, then sort them, along with
517 * eigenvectors.
518 *
519  IF( wantz ) THEN
520  DO 50 j = 1, m - 1
521  i = 0
522  tmp1 = w( j )
523  DO 40 jj = j + 1, m
524  IF( w( jj ).LT.tmp1 ) THEN
525  i = jj
526  tmp1 = w( jj )
527  END IF
528  40 CONTINUE
529 *
530  IF( i.NE.0 ) THEN
531  itmp1 = iwork( indibl+i-1 )
532  w( i ) = w( j )
533  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
534  w( j ) = tmp1
535  iwork( indibl+j-1 ) = itmp1
536  CALL cswap( n, z( 1, i ), 1, z( 1, j ), 1 )
537  IF( info.NE.0 ) THEN
538  itmp1 = ifail( i )
539  ifail( i ) = ifail( j )
540  ifail( j ) = itmp1
541  END IF
542  END IF
543  50 CONTINUE
544  END IF
545 *
546  RETURN
547 *
548 * End of CHBEVX
549 *
550  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:86
subroutine sstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
SSTEBZ
Definition: sstebz.f:273
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
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 cstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
CSTEIN
Definition: cstein.f:182
subroutine chbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
CHBTRD
Definition: chbtrd.f:163
subroutine csteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
CSTEQR
Definition: csteqr.f:132
subroutine chbevx(JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, RWORK, IWORK, IFAIL, INFO)
CHBEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition: chbevx.f:267
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79