LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ dsbevx_2stage()

 subroutine dsbevx_2stage ( character JOBZ, character RANGE, character UPLO, integer N, integer KD, double precision, dimension( ldab, * ) AB, integer LDAB, double precision, dimension( ldq, * ) Q, integer LDQ, 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 LWORK, integer, dimension( * ) IWORK, integer, dimension( * ) IFAIL, integer INFO )

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

Purpose:
``` DSBEVX_2STAGE computes selected eigenvalues and, optionally, eigenvectors
of a real symmetric band matrix A using the 2stage technique for
the reduction to tridiagonal. 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. Not available in this release.``` [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] UPLO ``` UPLO is CHARACTER*1 = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored.``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in] KD ``` KD is INTEGER The number of superdiagonals of the matrix A if UPLO = 'U', or the number of subdiagonals if UPLO = 'L'. KD >= 0.``` [in,out] AB ``` AB is DOUBLE PRECISION array, dimension (LDAB, N) On entry, the upper or lower triangle of the symmetric band matrix A, stored in the first KD+1 rows of the array. The j-th column of A is stored in the j-th column of the array AB as follows: if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). On exit, AB is overwritten by values generated during the reduction to tridiagonal form. If UPLO = 'U', the first superdiagonal and the diagonal of the tridiagonal matrix T are returned in rows KD and KD+1 of AB, and if UPLO = 'L', the diagonal and first subdiagonal of T are returned in the first two rows of AB.``` [in] LDAB ``` LDAB is INTEGER The leading dimension of the array AB. LDAB >= KD + 1.``` [out] Q ``` Q is DOUBLE PRECISION array, dimension (LDQ, N) If JOBZ = 'V', the N-by-N orthogonal matrix used in the reduction to tridiagonal form. If JOBZ = 'N', the array Q is not referenced.``` [in] LDQ ``` LDQ is INTEGER The leading dimension of the array Q. If JOBZ = 'V', then LDQ >= max(1,N).``` [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 obtained by reducing AB to tridiagonal form. 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 "Computing Small Singular Values of Bidiagonal Matrices with Guaranteed High Relative Accuracy," 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, 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 (LWORK)` [in] LWORK ``` LWORK is INTEGER The length of the array WORK. LWORK >= 1, when N <= 1; otherwise If JOBZ = 'N' and N > 1, LWORK must be queried. LWORK = MAX(1, 7*N, dimension) where dimension = (2KD+1)*N + KD*NTHREADS + 2*N where KD is the size of the band. NTHREADS is the number of threads used when openMP compilation is enabled, otherwise =1. If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA.``` [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.```
Further Details:
```  All details about the 2stage techniques are available in:

Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
Parallel reduction to condensed forms for symmetric eigenvalue problems
using aggregated fine-grained and memory-aware kernels. In Proceedings
of 2011 International Conference for High Performance Computing,
Networking, Storage and Analysis (SC '11), New York, NY, USA,
Article 8 , 11 pages.
http://doi.acm.org/10.1145/2063384.2063394

A. Haidar, J. Kurzak, P. Luszczek, 2013.
An improved parallel singular value algorithm and its implementation
for multicore hardware, In Proceedings of 2013 International Conference
for High Performance Computing, Networking, Storage and Analysis (SC '13).
Article 90, 12 pages.
http://doi.acm.org/10.1145/2503210.2503292

A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
A novel hybrid CPU-GPU generalized eigensolver for electronic structure
calculations based on fine-grained memory aware tasks.
International Journal of High Performance Computing Applications.
Volume 28 Issue 2, Pages 196-209, May 2014.
http://hpc.sagepub.com/content/28/2/196 ```

Definition at line 319 of file dsbevx_2stage.f.

322 *
323  IMPLICIT NONE
324 *
325 * -- LAPACK driver routine --
326 * -- LAPACK is a software package provided by Univ. of Tennessee, --
327 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
328 *
329 * .. Scalar Arguments ..
330  CHARACTER JOBZ, RANGE, UPLO
331  INTEGER IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N, LWORK
332  DOUBLE PRECISION ABSTOL, VL, VU
333 * ..
334 * .. Array Arguments ..
335  INTEGER IFAIL( * ), IWORK( * )
336  DOUBLE PRECISION AB( LDAB, * ), Q( LDQ, * ), W( * ), WORK( * ),
337  \$ Z( LDZ, * )
338 * ..
339 *
340 * =====================================================================
341 *
342 * .. Parameters ..
343  DOUBLE PRECISION ZERO, ONE
344  parameter( zero = 0.0d0, one = 1.0d0 )
345 * ..
346 * .. Local Scalars ..
347  LOGICAL ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ,
348  \$ LQUERY
349  CHARACTER ORDER
350  INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
351  \$ INDISP, INDIWO, INDWRK, ISCALE, ITMP1, J, JJ,
352  \$ LLWORK, LWMIN, LHTRD, LWTRD, IB, INDHOUS,
353  \$ NSPLIT
354  DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
355  \$ SIGMA, SMLNUM, TMP1, VLL, VUU
356 * ..
357 * .. External Functions ..
358  LOGICAL LSAME
359  INTEGER ILAENV2STAGE
360  DOUBLE PRECISION DLAMCH, DLANSB
361  EXTERNAL lsame, dlamch, dlansb, ilaenv2stage
362 * ..
363 * .. External Subroutines ..
364  EXTERNAL dcopy, dgemv, dlacpy, dlascl, dscal,
366  \$ dsytrd_sb2st
367 * ..
368 * .. Intrinsic Functions ..
369  INTRINSIC max, min, sqrt
370 * ..
371 * .. Executable Statements ..
372 *
373 * Test the input parameters.
374 *
375  wantz = lsame( jobz, 'V' )
376  alleig = lsame( range, 'A' )
377  valeig = lsame( range, 'V' )
378  indeig = lsame( range, 'I' )
379  lower = lsame( uplo, 'L' )
380  lquery = ( lwork.EQ.-1 )
381 *
382  info = 0
383  IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
384  info = -1
385  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
386  info = -2
387  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
388  info = -3
389  ELSE IF( n.LT.0 ) THEN
390  info = -4
391  ELSE IF( kd.LT.0 ) THEN
392  info = -5
393  ELSE IF( ldab.LT.kd+1 ) THEN
394  info = -7
395  ELSE IF( wantz .AND. ldq.LT.max( 1, n ) ) THEN
396  info = -9
397  ELSE
398  IF( valeig ) THEN
399  IF( n.GT.0 .AND. vu.LE.vl )
400  \$ info = -11
401  ELSE IF( indeig ) THEN
402  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
403  info = -12
404  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
405  info = -13
406  END IF
407  END IF
408  END IF
409  IF( info.EQ.0 ) THEN
410  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
411  \$ info = -18
412  END IF
413 *
414  IF( info.EQ.0 ) THEN
415  IF( n.LE.1 ) THEN
416  lwmin = 1
417  work( 1 ) = lwmin
418  ELSE
419  ib = ilaenv2stage( 2, 'DSYTRD_SB2ST', jobz,
420  \$ n, kd, -1, -1 )
421  lhtrd = ilaenv2stage( 3, 'DSYTRD_SB2ST', jobz,
422  \$ n, kd, ib, -1 )
423  lwtrd = ilaenv2stage( 4, 'DSYTRD_SB2ST', jobz,
424  \$ n, kd, ib, -1 )
425  lwmin = 2*n + lhtrd + lwtrd
426  work( 1 ) = lwmin
427  ENDIF
428 *
429  IF( lwork.LT.lwmin .AND. .NOT.lquery )
430  \$ info = -20
431  END IF
432 *
433  IF( info.NE.0 ) THEN
434  CALL xerbla( 'DSBEVX_2STAGE ', -info )
435  RETURN
436  ELSE IF( lquery ) THEN
437  RETURN
438  END IF
439 *
440 * Quick return if possible
441 *
442  m = 0
443  IF( n.EQ.0 )
444  \$ RETURN
445 *
446  IF( n.EQ.1 ) THEN
447  m = 1
448  IF( lower ) THEN
449  tmp1 = ab( 1, 1 )
450  ELSE
451  tmp1 = ab( kd+1, 1 )
452  END IF
453  IF( valeig ) THEN
454  IF( .NOT.( vl.LT.tmp1 .AND. vu.GE.tmp1 ) )
455  \$ m = 0
456  END IF
457  IF( m.EQ.1 ) THEN
458  w( 1 ) = tmp1
459  IF( wantz )
460  \$ z( 1, 1 ) = one
461  END IF
462  RETURN
463  END IF
464 *
465 * Get machine constants.
466 *
467  safmin = dlamch( 'Safe minimum' )
468  eps = dlamch( 'Precision' )
469  smlnum = safmin / eps
470  bignum = one / smlnum
471  rmin = sqrt( smlnum )
472  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
473 *
474 * Scale matrix to allowable range, if necessary.
475 *
476  iscale = 0
477  abstll = abstol
478  IF( valeig ) THEN
479  vll = vl
480  vuu = vu
481  ELSE
482  vll = zero
483  vuu = zero
484  END IF
485  anrm = dlansb( 'M', uplo, n, kd, ab, ldab, work )
486  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
487  iscale = 1
488  sigma = rmin / anrm
489  ELSE IF( anrm.GT.rmax ) THEN
490  iscale = 1
491  sigma = rmax / anrm
492  END IF
493  IF( iscale.EQ.1 ) THEN
494  IF( lower ) THEN
495  CALL dlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
496  ELSE
497  CALL dlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
498  END IF
499  IF( abstol.GT.0 )
500  \$ abstll = abstol*sigma
501  IF( valeig ) THEN
502  vll = vl*sigma
503  vuu = vu*sigma
504  END IF
505  END IF
506 *
507 * Call DSYTRD_SB2ST to reduce symmetric band matrix to tridiagonal form.
508 *
509  indd = 1
510  inde = indd + n
511  indhous = inde + n
512  indwrk = indhous + lhtrd
513  llwork = lwork - indwrk + 1
514 *
515  CALL dsytrd_sb2st( "N", jobz, uplo, n, kd, ab, ldab, work( indd ),
516  \$ work( inde ), work( indhous ), lhtrd,
517  \$ work( indwrk ), llwork, iinfo )
518 *
519 * If all eigenvalues are desired and ABSTOL is less than or equal
520 * to zero, then call DSTERF or SSTEQR. If this fails for some
521 * eigenvalue, then try DSTEBZ.
522 *
523  test = .false.
524  IF (indeig) THEN
525  IF (il.EQ.1 .AND. iu.EQ.n) THEN
526  test = .true.
527  END IF
528  END IF
529  IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
530  CALL dcopy( n, work( indd ), 1, w, 1 )
531  indee = indwrk + 2*n
532  IF( .NOT.wantz ) THEN
533  CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
534  CALL dsterf( n, w, work( indee ), info )
535  ELSE
536  CALL dlacpy( 'A', n, n, q, ldq, z, ldz )
537  CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
538  CALL dsteqr( jobz, n, w, work( indee ), z, ldz,
539  \$ work( indwrk ), info )
540  IF( info.EQ.0 ) THEN
541  DO 10 i = 1, n
542  ifail( i ) = 0
543  10 CONTINUE
544  END IF
545  END IF
546  IF( info.EQ.0 ) THEN
547  m = n
548  GO TO 30
549  END IF
550  info = 0
551  END IF
552 *
553 * Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
554 *
555  IF( wantz ) THEN
556  order = 'B'
557  ELSE
558  order = 'E'
559  END IF
560  indibl = 1
561  indisp = indibl + n
562  indiwo = indisp + n
563  CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
564  \$ work( indd ), work( inde ), m, nsplit, w,
565  \$ iwork( indibl ), iwork( indisp ), work( indwrk ),
566  \$ iwork( indiwo ), info )
567 *
568  IF( wantz ) THEN
569  CALL dstein( n, work( indd ), work( inde ), m, w,
570  \$ iwork( indibl ), iwork( indisp ), z, ldz,
571  \$ work( indwrk ), iwork( indiwo ), ifail, info )
572 *
573 * Apply orthogonal matrix used in reduction to tridiagonal
574 * form to eigenvectors returned by DSTEIN.
575 *
576  DO 20 j = 1, m
577  CALL dcopy( n, z( 1, j ), 1, work( 1 ), 1 )
578  CALL dgemv( 'N', n, n, one, q, ldq, work, 1, zero,
579  \$ z( 1, j ), 1 )
580  20 CONTINUE
581  END IF
582 *
583 * If matrix was scaled, then rescale eigenvalues appropriately.
584 *
585  30 CONTINUE
586  IF( iscale.EQ.1 ) THEN
587  IF( info.EQ.0 ) THEN
588  imax = m
589  ELSE
590  imax = info - 1
591  END IF
592  CALL dscal( imax, one / sigma, w, 1 )
593  END IF
594 *
595 * If eigenvalues are not in order, then sort them, along with
596 * eigenvectors.
597 *
598  IF( wantz ) THEN
599  DO 50 j = 1, m - 1
600  i = 0
601  tmp1 = w( j )
602  DO 40 jj = j + 1, m
603  IF( w( jj ).LT.tmp1 ) THEN
604  i = jj
605  tmp1 = w( jj )
606  END IF
607  40 CONTINUE
608 *
609  IF( i.NE.0 ) THEN
610  itmp1 = iwork( indibl+i-1 )
611  w( i ) = w( j )
612  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
613  w( j ) = tmp1
614  iwork( indibl+j-1 ) = itmp1
615  CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
616  IF( info.NE.0 ) THEN
617  itmp1 = ifail( i )
618  ifail( i ) = ifail( j )
619  ifail( j ) = itmp1
620  END IF
621  END IF
622  50 CONTINUE
623  END IF
624 *
625 * Set WORK(1) to optimal workspace size.
626 *
627  work( 1 ) = lwmin
628 *
629  RETURN
630 *
631 * End of DSBEVX_2STAGE
632 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dsytrd_sb2st(STAGE1, VECT, UPLO, N, KD, AB, LDAB, D, E, HOUS, LHOUS, WORK, LWORK, INFO)
DSYTRD_SB2ST reduces a real symmetric band matrix A to real symmetric tridiagonal form T
Definition: dsytrd_sb2st.F:230
integer function ilaenv2stage(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV2STAGE
Definition: ilaenv2stage.f:149
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:143
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:273
subroutine dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DSTEQR
Definition: dsteqr.f:131
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:86
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:156
double precision function dlansb(NORM, UPLO, N, K, AB, LDAB, WORK)
DLANSB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: dlansb.f:129
subroutine dstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
DSTEIN
Definition: dstein.f:174
Here is the call graph for this function:
Here is the caller graph for this function: