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

◆ 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

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

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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).
  Denver, Colorado, USA, 2013.
  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,
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*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
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
integer function ilaenv2stage(ispec, name, opts, n1, n2, n3, n4)
ILAENV2STAGE
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
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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 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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
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 dstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
DSTEIN
Definition dstein.f:174
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 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: