LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dstevx()

subroutine dstevx ( character jobz,
character range,
integer n,
double precision, dimension( * ) d,
double precision, dimension( * ) e,
double precision vl,
double precision vu,
integer il,
integer iu,
double precision abstol,
integer m,
double precision, dimension( * ) w,
double precision, dimension( ldz, * ) z,
integer ldz,
double precision, dimension( * ) work,
integer, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

DSTEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices

Download DSTEVX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DSTEVX computes selected eigenvalues and, optionally, eigenvectors
!> of a real symmetric tridiagonal matrix A.  Eigenvalues and
!> eigenvectors can be selected by specifying either a range of values
!> or a range of indices for the desired eigenvalues.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]RANGE
!>          RANGE is CHARACTER*1
!>          = 'A': all eigenvalues will be found.
!>          = 'V': all eigenvalues in the half-open interval (VL,VU]
!>                 will be found.
!>          = 'I': the IL-th through IU-th eigenvalues will be found.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix.  N >= 0.
!>
[in,out]D
!>          D is DOUBLE PRECISION array, dimension (N)
!>          On entry, the n diagonal elements of the tridiagonal matrix
!>          A.
!>          On exit, D may be multiplied by a constant factor chosen
!>          to avoid over/underflow in computing the eigenvalues.
!>
[in,out]E
!>          E is DOUBLE PRECISION array, dimension (max(1,N-1))
!>          On entry, the (n-1) subdiagonal elements of the tridiagonal
!>          matrix A in elements 1 to N-1 of E.
!>          On exit, E may be multiplied by a constant factor chosen
!>          to avoid over/underflow in computing the eigenvalues.
!>
[in]VL
!>          VL is DOUBLE PRECISION
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less
!>          than or equal to zero, then  EPS*|T|  will be used in
!>          its place, where |T| is the 1-norm of the tridiagonal
!>          matrix.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*DLAMCH('S').
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          The first M elements contain the selected eigenvalues in
!>          ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          If an eigenvector fails to converge (INFO > 0), then that
!>          column of Z contains the latest approximation to the
!>          eigenvector, and the index of the eigenvector is returned
!>          in IFAIL.  If JOBZ = 'N', then Z is not referenced.
!>          Note: the user must ensure that at least max(1,M) columns are
!>          supplied in the array Z; if RANGE = 'V', the exact value of M
!>          is not known in advance and an upper bound must be used.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (5*N)
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (5*N)
!>
[out]IFAIL
!>          IFAIL is INTEGER array, dimension (N)
!>          If JOBZ = 'V', then if INFO = 0, the first M elements of
!>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
!>          indices of the eigenvectors that failed to converge.
!>          If JOBZ = 'N', then IFAIL is not referenced.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, then i eigenvectors failed to converge.
!>                Their indices are stored in array IFAIL.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 223 of file dstevx.f.

226*
227* -- LAPACK driver routine --
228* -- LAPACK is a software package provided by Univ. of Tennessee, --
229* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
230*
231* .. Scalar Arguments ..
232 CHARACTER JOBZ, RANGE
233 INTEGER IL, INFO, IU, LDZ, M, N
234 DOUBLE PRECISION ABSTOL, VL, VU
235* ..
236* .. Array Arguments ..
237 INTEGER IFAIL( * ), IWORK( * )
238 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
239* ..
240*
241* =====================================================================
242*
243* .. Parameters ..
244 DOUBLE PRECISION ZERO, ONE
245 parameter( zero = 0.0d0, one = 1.0d0 )
246* ..
247* .. Local Scalars ..
248 LOGICAL ALLEIG, INDEIG, TEST, VALEIG, WANTZ
249 CHARACTER ORDER
250 INTEGER I, IMAX, INDISP, INDIWO, INDWRK,
251 $ ISCALE, ITMP1, J, JJ, NSPLIT
252 DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
253 $ TMP1, TNRM, VLL, VUU
254* ..
255* .. External Functions ..
256 LOGICAL LSAME
257 DOUBLE PRECISION DLAMCH, DLANST
258 EXTERNAL lsame, dlamch, dlanst
259* ..
260* .. External Subroutines ..
261 EXTERNAL dcopy, dscal, dstebz, dstein, dsteqr,
262 $ dsterf,
263 $ dswap, xerbla
264* ..
265* .. Intrinsic Functions ..
266 INTRINSIC max, min, sqrt
267* ..
268* .. Executable Statements ..
269*
270* Test the input parameters.
271*
272 wantz = lsame( jobz, 'V' )
273 alleig = lsame( range, 'A' )
274 valeig = lsame( range, 'V' )
275 indeig = lsame( range, 'I' )
276*
277 info = 0
278 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
279 info = -1
280 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
281 info = -2
282 ELSE IF( n.LT.0 ) THEN
283 info = -3
284 ELSE
285 IF( valeig ) THEN
286 IF( n.GT.0 .AND. vu.LE.vl )
287 $ info = -7
288 ELSE IF( indeig ) THEN
289 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
290 info = -8
291 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
292 info = -9
293 END IF
294 END IF
295 END IF
296 IF( info.EQ.0 ) THEN
297 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
298 $ info = -14
299 END IF
300*
301 IF( info.NE.0 ) THEN
302 CALL xerbla( 'DSTEVX', -info )
303 RETURN
304 END IF
305*
306* Quick return if possible
307*
308 m = 0
309 IF( n.EQ.0 )
310 $ RETURN
311*
312 IF( n.EQ.1 ) THEN
313 IF( alleig .OR. indeig ) THEN
314 m = 1
315 w( 1 ) = d( 1 )
316 ELSE
317 IF( vl.LT.d( 1 ) .AND. vu.GE.d( 1 ) ) THEN
318 m = 1
319 w( 1 ) = d( 1 )
320 END IF
321 END IF
322 IF( wantz )
323 $ z( 1, 1 ) = one
324 RETURN
325 END IF
326*
327* Get machine constants.
328*
329 safmin = dlamch( 'Safe minimum' )
330 eps = dlamch( 'Precision' )
331 smlnum = safmin / eps
332 bignum = one / smlnum
333 rmin = sqrt( smlnum )
334 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
335*
336* Scale matrix to allowable range, if necessary.
337*
338 iscale = 0
339 IF( valeig ) THEN
340 vll = vl
341 vuu = vu
342 ELSE
343 vll = zero
344 vuu = zero
345 END IF
346 tnrm = dlanst( 'M', n, d, e )
347 IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
348 iscale = 1
349 sigma = rmin / tnrm
350 ELSE IF( tnrm.GT.rmax ) THEN
351 iscale = 1
352 sigma = rmax / tnrm
353 END IF
354 IF( iscale.EQ.1 ) THEN
355 CALL dscal( n, sigma, d, 1 )
356 CALL dscal( n-1, sigma, e( 1 ), 1 )
357 IF( valeig ) THEN
358 vll = vl*sigma
359 vuu = vu*sigma
360 END IF
361 END IF
362*
363* If all eigenvalues are desired and ABSTOL is less than zero, then
364* call DSTERF or SSTEQR. If this fails for some eigenvalue, then
365* try DSTEBZ.
366*
367 test = .false.
368 IF( indeig ) THEN
369 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
370 test = .true.
371 END IF
372 END IF
373 IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
374 CALL dcopy( n, d, 1, w, 1 )
375 CALL dcopy( n-1, e( 1 ), 1, work( 1 ), 1 )
376 indwrk = n + 1
377 IF( .NOT.wantz ) THEN
378 CALL dsterf( n, w, work, info )
379 ELSE
380 CALL dsteqr( 'I', n, w, work, z, ldz, work( indwrk ),
381 $ info )
382 IF( info.EQ.0 ) THEN
383 DO 10 i = 1, n
384 ifail( i ) = 0
385 10 CONTINUE
386 END IF
387 END IF
388 IF( info.EQ.0 ) THEN
389 m = n
390 GO TO 20
391 END IF
392 info = 0
393 END IF
394*
395* Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
396*
397 IF( wantz ) THEN
398 order = 'B'
399 ELSE
400 order = 'E'
401 END IF
402 indwrk = 1
403 indisp = 1 + n
404 indiwo = indisp + n
405 CALL dstebz( range, order, n, vll, vuu, il, iu, abstol, d, e,
406 $ m,
407 $ nsplit, w, iwork( 1 ), iwork( indisp ),
408 $ work( indwrk ), iwork( indiwo ), info )
409*
410 IF( wantz ) THEN
411 CALL dstein( n, d, e, m, w, iwork( 1 ), iwork( indisp ),
412 $ z, ldz, work( indwrk ), iwork( indiwo ), ifail,
413 $ info )
414 END IF
415*
416* If matrix was scaled, then rescale eigenvalues appropriately.
417*
418 20 CONTINUE
419 IF( iscale.EQ.1 ) THEN
420 IF( info.EQ.0 ) THEN
421 imax = m
422 ELSE
423 imax = info - 1
424 END IF
425 CALL dscal( imax, one / sigma, w, 1 )
426 END IF
427*
428* If eigenvalues are not in order, then sort them, along with
429* eigenvectors.
430*
431 IF( wantz ) THEN
432 DO 40 j = 1, m - 1
433 i = 0
434 tmp1 = w( j )
435 DO 30 jj = j + 1, m
436 IF( w( jj ).LT.tmp1 ) THEN
437 i = jj
438 tmp1 = w( jj )
439 END IF
440 30 CONTINUE
441*
442 IF( i.NE.0 ) THEN
443 itmp1 = iwork( 1 + i-1 )
444 w( i ) = w( j )
445 iwork( 1 + i-1 ) = iwork( 1 + j-1 )
446 w( j ) = tmp1
447 iwork( 1 + j-1 ) = itmp1
448 CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
449 IF( info.NE.0 ) THEN
450 itmp1 = ifail( i )
451 ifail( i ) = ifail( j )
452 ifail( j ) = itmp1
453 END IF
454 END IF
455 40 CONTINUE
456 END IF
457*
458 RETURN
459*
460* End of DSTEVX
461*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlanst
double precision function dlanst(norm, n, d, e)
DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition dlanst.f:98
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dstebz
subroutine dstebz(range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, iwork, info)
DSTEBZ
Definition dstebz.f:272
dstein
subroutine dstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
DSTEIN
Definition dstein.f:172
dsteqr
subroutine dsteqr(compz, n, d, e, z, ldz, work, info)
DSTEQR
Definition dsteqr.f:129
dsterf
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:84
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: