LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
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
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