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