LAPACK  3.8.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 *> \date June 2016
260 *
261 *> \ingroup realOTHEReigen
262 *
263 * =====================================================================
264  SUBROUTINE ssbevx( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ, VL,
265  $ VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, IWORK,
266  $ IFAIL, INFO )
267 *
268 * -- LAPACK driver routine (version 3.7.0) --
269 * -- LAPACK is a software package provided by Univ. of Tennessee, --
270 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
271 * June 2016
272 *
273 * .. Scalar Arguments ..
274  CHARACTER JOBZ, RANGE, UPLO
275  INTEGER IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N
276  REAL ABSTOL, VL, VU
277 * ..
278 * .. Array Arguments ..
279  INTEGER IFAIL( * ), IWORK( * )
280  REAL AB( ldab, * ), Q( ldq, * ), W( * ), WORK( * ),
281  $ z( ldz, * )
282 * ..
283 *
284 * =====================================================================
285 *
286 * .. Parameters ..
287  REAL ZERO, ONE
288  parameter( zero = 0.0e0, one = 1.0e0 )
289 * ..
290 * .. Local Scalars ..
291  LOGICAL ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ
292  CHARACTER ORDER
293  INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
294  $ indisp, indiwo, indwrk, iscale, itmp1, j, jj,
295  $ nsplit
296  REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
297  $ sigma, smlnum, tmp1, vll, vuu
298 * ..
299 * .. External Functions ..
300  LOGICAL LSAME
301  REAL SLAMCH, SLANSB
302  EXTERNAL lsame, slamch, slansb
303 * ..
304 * .. External Subroutines ..
305  EXTERNAL scopy, sgemv, slacpy, slascl, ssbtrd, sscal,
307 * ..
308 * .. Intrinsic Functions ..
309  INTRINSIC max, min, sqrt
310 * ..
311 * .. Executable Statements ..
312 *
313 * Test the input parameters.
314 *
315  wantz = lsame( jobz, 'V' )
316  alleig = lsame( range, 'A' )
317  valeig = lsame( range, 'V' )
318  indeig = lsame( range, 'I' )
319  lower = lsame( uplo, 'L' )
320 *
321  info = 0
322  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
323  info = -1
324  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
325  info = -2
326  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
327  info = -3
328  ELSE IF( n.LT.0 ) THEN
329  info = -4
330  ELSE IF( kd.LT.0 ) THEN
331  info = -5
332  ELSE IF( ldab.LT.kd+1 ) THEN
333  info = -7
334  ELSE IF( wantz .AND. ldq.LT.max( 1, n ) ) THEN
335  info = -9
336  ELSE
337  IF( valeig ) THEN
338  IF( n.GT.0 .AND. vu.LE.vl )
339  $ info = -11
340  ELSE IF( indeig ) THEN
341  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
342  info = -12
343  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
344  info = -13
345  END IF
346  END IF
347  END IF
348  IF( info.EQ.0 ) THEN
349  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
350  $ info = -18
351  END IF
352 *
353  IF( info.NE.0 ) THEN
354  CALL xerbla( 'SSBEVX', -info )
355  RETURN
356  END IF
357 *
358 * Quick return if possible
359 *
360  m = 0
361  IF( n.EQ.0 )
362  $ RETURN
363 *
364  IF( n.EQ.1 ) THEN
365  m = 1
366  IF( lower ) THEN
367  tmp1 = ab( 1, 1 )
368  ELSE
369  tmp1 = ab( kd+1, 1 )
370  END IF
371  IF( valeig ) THEN
372  IF( .NOT.( vl.LT.tmp1 .AND. vu.GE.tmp1 ) )
373  $ m = 0
374  END IF
375  IF( m.EQ.1 ) THEN
376  w( 1 ) = tmp1
377  IF( wantz )
378  $ z( 1, 1 ) = one
379  END IF
380  RETURN
381  END IF
382 *
383 * Get machine constants.
384 *
385  safmin = slamch( 'Safe minimum' )
386  eps = slamch( 'Precision' )
387  smlnum = safmin / eps
388  bignum = one / smlnum
389  rmin = sqrt( smlnum )
390  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
391 *
392 * Scale matrix to allowable range, if necessary.
393 *
394  iscale = 0
395  abstll = abstol
396  IF ( valeig ) THEN
397  vll = vl
398  vuu = vu
399  ELSE
400  vll = zero
401  vuu = zero
402  ENDIF
403  anrm = slansb( 'M', uplo, n, kd, ab, ldab, work )
404  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
405  iscale = 1
406  sigma = rmin / anrm
407  ELSE IF( anrm.GT.rmax ) THEN
408  iscale = 1
409  sigma = rmax / anrm
410  END IF
411  IF( iscale.EQ.1 ) THEN
412  IF( lower ) THEN
413  CALL slascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
414  ELSE
415  CALL slascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
416  END IF
417  IF( abstol.GT.0 )
418  $ abstll = abstol*sigma
419  IF( valeig ) THEN
420  vll = vl*sigma
421  vuu = vu*sigma
422  END IF
423  END IF
424 *
425 * Call SSBTRD to reduce symmetric band matrix to tridiagonal form.
426 *
427  indd = 1
428  inde = indd + n
429  indwrk = inde + n
430  CALL ssbtrd( jobz, uplo, n, kd, ab, ldab, work( indd ),
431  $ work( inde ), q, ldq, work( indwrk ), iinfo )
432 *
433 * If all eigenvalues are desired and ABSTOL is less than or equal
434 * to zero, then call SSTERF or SSTEQR. If this fails for some
435 * eigenvalue, then try SSTEBZ.
436 *
437  test = .false.
438  IF (indeig) THEN
439  IF (il.EQ.1 .AND. iu.EQ.n) THEN
440  test = .true.
441  END IF
442  END IF
443  IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
444  CALL scopy( n, work( indd ), 1, w, 1 )
445  indee = indwrk + 2*n
446  IF( .NOT.wantz ) THEN
447  CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
448  CALL ssterf( n, w, work( indee ), info )
449  ELSE
450  CALL slacpy( 'A', n, n, q, ldq, z, ldz )
451  CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
452  CALL ssteqr( jobz, n, w, work( indee ), z, ldz,
453  $ work( indwrk ), info )
454  IF( info.EQ.0 ) THEN
455  DO 10 i = 1, n
456  ifail( i ) = 0
457  10 CONTINUE
458  END IF
459  END IF
460  IF( info.EQ.0 ) THEN
461  m = n
462  GO TO 30
463  END IF
464  info = 0
465  END IF
466 *
467 * Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
468 *
469  IF( wantz ) THEN
470  order = 'B'
471  ELSE
472  order = 'E'
473  END IF
474  indibl = 1
475  indisp = indibl + n
476  indiwo = indisp + n
477  CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
478  $ work( indd ), work( inde ), m, nsplit, w,
479  $ iwork( indibl ), iwork( indisp ), work( indwrk ),
480  $ iwork( indiwo ), info )
481 *
482  IF( wantz ) THEN
483  CALL sstein( n, work( indd ), work( inde ), m, w,
484  $ iwork( indibl ), iwork( indisp ), z, ldz,
485  $ work( indwrk ), iwork( indiwo ), ifail, info )
486 *
487 * Apply orthogonal matrix used in reduction to tridiagonal
488 * form to eigenvectors returned by SSTEIN.
489 *
490  DO 20 j = 1, m
491  CALL scopy( n, z( 1, j ), 1, work( 1 ), 1 )
492  CALL sgemv( 'N', n, n, one, q, ldq, work, 1, zero,
493  $ z( 1, j ), 1 )
494  20 CONTINUE
495  END IF
496 *
497 * If matrix was scaled, then rescale eigenvalues appropriately.
498 *
499  30 CONTINUE
500  IF( iscale.EQ.1 ) THEN
501  IF( info.EQ.0 ) THEN
502  imax = m
503  ELSE
504  imax = info - 1
505  END IF
506  CALL sscal( imax, one / sigma, w, 1 )
507  END IF
508 *
509 * If eigenvalues are not in order, then sort them, along with
510 * eigenvectors.
511 *
512  IF( wantz ) THEN
513  DO 50 j = 1, m - 1
514  i = 0
515  tmp1 = w( j )
516  DO 40 jj = j + 1, m
517  IF( w( jj ).LT.tmp1 ) THEN
518  i = jj
519  tmp1 = w( jj )
520  END IF
521  40 CONTINUE
522 *
523  IF( i.NE.0 ) THEN
524  itmp1 = iwork( indibl+i-1 )
525  w( i ) = w( j )
526  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
527  w( j ) = tmp1
528  iwork( indibl+j-1 ) = itmp1
529  CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
530  IF( info.NE.0 ) THEN
531  itmp1 = ifail( i )
532  ifail( i ) = ifail( j )
533  ifail( j ) = itmp1
534  END IF
535  END IF
536  50 CONTINUE
537  END IF
538 *
539  RETURN
540 *
541 * End of SSBEVX
542 *
543  END
subroutine sstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
SSTEBZ
Definition: sstebz.f:275
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 matric...
Definition: ssbevx.f:267
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:133
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:158
subroutine ssbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
SSBTRD
Definition: ssbtrd.f:165
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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:145
subroutine sstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
SSTEIN
Definition: sstein.f:176
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:81
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:84
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:88
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:84