LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
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 *> \endverbatim
127 *>
128 *> \param[in] VU
129 *> \verbatim
130 *> VU is REAL
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 REAL
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*SLAMCH('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*SLAMCH('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 REAL 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 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 array, dimension (N)
216 *> \endverbatim
217 *>
218 *> \param[out] RWORK
219 *> \verbatim
220 *> RWORK is REAL 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 complexOTHEReigen
257 *
258 * =====================================================================
259  SUBROUTINE chbevx( 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  REAL abstol, vl, vu
272 * ..
273 * .. Array Arguments ..
274  INTEGER ifail( * ), iwork( * )
275  REAL rwork( * ), w( * )
276  COMPLEX ab( ldab, * ), q( ldq, * ), work( * ),
277  $ z( ldz, * )
278 * ..
279 *
280 * =====================================================================
281 *
282 * .. Parameters ..
283  REAL zero, one
284  parameter( zero = 0.0e0, one = 1.0e0 )
285  COMPLEX czero, cone
286  parameter( czero = ( 0.0e0, 0.0e0 ),
287  $ cone = ( 1.0e0, 0.0e0 ) )
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  REAL abstll, anrm, bignum, eps, rmax, rmin, safmin,
296  $ sigma, smlnum, tmp1, vll, vuu
297  COMPLEX ctmp1
298 * ..
299 * .. External Functions ..
300  LOGICAL lsame
301  REAL clanhb, slamch
302  EXTERNAL lsame, clanhb, slamch
303 * ..
304 * .. External Subroutines ..
305  EXTERNAL ccopy, cgemv, chbtrd, clacpy, clascl, cstein,
307  $ xerbla
308 * ..
309 * .. Intrinsic Functions ..
310  INTRINSIC max, min, REAL, 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( 'CHBEVX', -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 = REAL( 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 = slamch( 'Safe minimum' )
388  eps = slamch( '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  ENDIF
405  anrm = clanhb( '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 clascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
416  ELSE
417  CALL clascl( '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 CHBTRD 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 chbtrd( 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 SSTERF or CSTEQR. If this fails for some
438 * eigenvalue, then try SSTEBZ.
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 scopy( n, rwork( indd ), 1, w, 1 )
448  indee = indrwk + 2*n
449  IF( .NOT.wantz ) THEN
450  CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
451  CALL ssterf( n, w, rwork( indee ), info )
452  ELSE
453  CALL clacpy( 'A', n, n, q, ldq, z, ldz )
454  CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
455  CALL csteqr( 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 SSTEBZ and, if eigenvectors are desired, CSTEIN.
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 sstebz( 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 cstein( 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 CSTEIN.
492 *
493  DO 20 j = 1, m
494  CALL ccopy( n, z( 1, j ), 1, work( 1 ), 1 )
495  CALL cgemv( '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 sscal( 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 cswap( 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 CHBEVX
545 *
546  END
LOGICAL function lsame(CA, CB)
LSAME
Definition: lsame.f:54
REAL function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:61
subroutine chbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
CHBTRD
Definition: chbtrd.f:163
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:51
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:52
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:159
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
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:140
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
Definition: xerbla-fortran:9
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:51
REAL function clanhb(NORM, UPLO, N, K, AB, LDAB, WORK)
CLANHB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a Hermitian band matrix.
Definition: clanhb.f:132
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:87
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 matric...
Definition: chbevx.f:259
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:54