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

◆ ssbevx_2stage()

subroutine ssbevx_2stage ( character  jobz,
character  range,
character  uplo,
integer  n,
integer  kd,
real, dimension( ldab, * )  ab,
integer  ldab,
real, dimension( ldq, * )  q,
integer  ldq,
real  vl,
real  vu,
integer  il,
integer  iu,
real  abstol,
integer  m,
real, dimension( * )  w,
real, dimension( ldz, * )  z,
integer  ldz,
real, dimension( * )  work,
integer  lwork,
integer, dimension( * )  iwork,
integer, dimension( * )  ifail,
integer  info 
)

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

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

Purpose:
 SSBEVX_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 REAL 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 REAL 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 REAL
          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 REAL
          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 REAL
          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*SLAMCH('S'), not zero.
          If this routine returns with INFO>0, indicating that some
          eigenvectors did not converge, try setting ABSTOL to
          2*SLAMCH('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 REAL array, dimension (N)
          The first M elements contain the selected eigenvalues in
          ascending order.
[out]Z
          Z is REAL 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 REAL 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 ssbevx_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 REAL ABSTOL, VL, VU
333* ..
334* .. Array Arguments ..
335 INTEGER IFAIL( * ), IWORK( * )
336 REAL AB( LDAB, * ), Q( LDQ, * ), W( * ), WORK( * ),
337 $ Z( LDZ, * )
338* ..
339*
340* =====================================================================
341*
342* .. Parameters ..
343 REAL ZERO, ONE
344 parameter( zero = 0.0e0, one = 1.0e0 )
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 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
355 $ SIGMA, SMLNUM, TMP1, VLL, VUU
356* ..
357* .. External Functions ..
358 LOGICAL LSAME
359 INTEGER ILAENV2STAGE
360 REAL SLAMCH, SLANSB, SROUNDUP_LWORK
361 EXTERNAL lsame, slamch, slansb, ilaenv2stage,
363* ..
364* .. External Subroutines ..
365 EXTERNAL scopy, sgemv, slacpy, slascl, sscal,
368* ..
369* .. Intrinsic Functions ..
370 INTRINSIC max, min, sqrt
371* ..
372* .. Executable Statements ..
373*
374* Test the input parameters.
375*
376 wantz = lsame( jobz, 'V' )
377 alleig = lsame( range, 'A' )
378 valeig = lsame( range, 'V' )
379 indeig = lsame( range, 'I' )
380 lower = lsame( uplo, 'L' )
381 lquery = ( lwork.EQ.-1 )
382*
383 info = 0
384 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
385 info = -1
386 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
387 info = -2
388 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
389 info = -3
390 ELSE IF( n.LT.0 ) THEN
391 info = -4
392 ELSE IF( kd.LT.0 ) THEN
393 info = -5
394 ELSE IF( ldab.LT.kd+1 ) THEN
395 info = -7
396 ELSE IF( wantz .AND. ldq.LT.max( 1, n ) ) THEN
397 info = -9
398 ELSE
399 IF( valeig ) THEN
400 IF( n.GT.0 .AND. vu.LE.vl )
401 $ info = -11
402 ELSE IF( indeig ) THEN
403 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
404 info = -12
405 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
406 info = -13
407 END IF
408 END IF
409 END IF
410 IF( info.EQ.0 ) THEN
411 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
412 $ info = -18
413 END IF
414*
415 IF( info.EQ.0 ) THEN
416 IF( n.LE.1 ) THEN
417 lwmin = 1
418 work( 1 ) = sroundup_lwork(lwmin)
419 ELSE
420 ib = ilaenv2stage( 2, 'SSYTRD_SB2ST', jobz,
421 $ n, kd, -1, -1 )
422 lhtrd = ilaenv2stage( 3, 'SSYTRD_SB2ST', jobz,
423 $ n, kd, ib, -1 )
424 lwtrd = ilaenv2stage( 4, 'SSYTRD_SB2ST', jobz,
425 $ n, kd, ib, -1 )
426 lwmin = 2*n + lhtrd + lwtrd
427 work( 1 ) = sroundup_lwork(lwmin)
428 ENDIF
429*
430 IF( lwork.LT.lwmin .AND. .NOT.lquery )
431 $ info = -20
432 END IF
433*
434 IF( info.NE.0 ) THEN
435 CALL xerbla( 'SSBEVX_2STAGE ', -info )
436 RETURN
437 ELSE IF( lquery ) THEN
438 RETURN
439 END IF
440*
441* Quick return if possible
442*
443 m = 0
444 IF( n.EQ.0 )
445 $ RETURN
446*
447 IF( n.EQ.1 ) THEN
448 m = 1
449 IF( lower ) THEN
450 tmp1 = ab( 1, 1 )
451 ELSE
452 tmp1 = ab( kd+1, 1 )
453 END IF
454 IF( valeig ) THEN
455 IF( .NOT.( vl.LT.tmp1 .AND. vu.GE.tmp1 ) )
456 $ m = 0
457 END IF
458 IF( m.EQ.1 ) THEN
459 w( 1 ) = tmp1
460 IF( wantz )
461 $ z( 1, 1 ) = one
462 END IF
463 RETURN
464 END IF
465*
466* Get machine constants.
467*
468 safmin = slamch( 'Safe minimum' )
469 eps = slamch( 'Precision' )
470 smlnum = safmin / eps
471 bignum = one / smlnum
472 rmin = sqrt( smlnum )
473 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
474*
475* Scale matrix to allowable range, if necessary.
476*
477 iscale = 0
478 abstll = abstol
479 IF( valeig ) THEN
480 vll = vl
481 vuu = vu
482 ELSE
483 vll = zero
484 vuu = zero
485 END IF
486 anrm = slansb( 'M', uplo, n, kd, ab, ldab, work )
487 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
488 iscale = 1
489 sigma = rmin / anrm
490 ELSE IF( anrm.GT.rmax ) THEN
491 iscale = 1
492 sigma = rmax / anrm
493 END IF
494 IF( iscale.EQ.1 ) THEN
495 IF( lower ) THEN
496 CALL slascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
497 ELSE
498 CALL slascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
499 END IF
500 IF( abstol.GT.0 )
501 $ abstll = abstol*sigma
502 IF( valeig ) THEN
503 vll = vl*sigma
504 vuu = vu*sigma
505 END IF
506 END IF
507*
508* Call SSYTRD_SB2ST to reduce symmetric band matrix to tridiagonal form.
509*
510 indd = 1
511 inde = indd + n
512 indhous = inde + n
513 indwrk = indhous + lhtrd
514 llwork = lwork - indwrk + 1
515*
516 CALL ssytrd_sb2st( "N", jobz, uplo, n, kd, ab, ldab, work( indd ),
517 $ work( inde ), work( indhous ), lhtrd,
518 $ work( indwrk ), llwork, iinfo )
519*
520* If all eigenvalues are desired and ABSTOL is less than or equal
521* to zero, then call SSTERF or SSTEQR. If this fails for some
522* eigenvalue, then try SSTEBZ.
523*
524 test = .false.
525 IF (indeig) THEN
526 IF (il.EQ.1 .AND. iu.EQ.n) THEN
527 test = .true.
528 END IF
529 END IF
530 IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
531 CALL scopy( n, work( indd ), 1, w, 1 )
532 indee = indwrk + 2*n
533 IF( .NOT.wantz ) THEN
534 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
535 CALL ssterf( n, w, work( indee ), info )
536 ELSE
537 CALL slacpy( 'A', n, n, q, ldq, z, ldz )
538 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
539 CALL ssteqr( jobz, n, w, work( indee ), z, ldz,
540 $ work( indwrk ), info )
541 IF( info.EQ.0 ) THEN
542 DO 10 i = 1, n
543 ifail( i ) = 0
544 10 CONTINUE
545 END IF
546 END IF
547 IF( info.EQ.0 ) THEN
548 m = n
549 GO TO 30
550 END IF
551 info = 0
552 END IF
553*
554* Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
555*
556 IF( wantz ) THEN
557 order = 'B'
558 ELSE
559 order = 'E'
560 END IF
561 indibl = 1
562 indisp = indibl + n
563 indiwo = indisp + n
564 CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
565 $ work( indd ), work( inde ), m, nsplit, w,
566 $ iwork( indibl ), iwork( indisp ), work( indwrk ),
567 $ iwork( indiwo ), info )
568*
569 IF( wantz ) THEN
570 CALL sstein( n, work( indd ), work( inde ), m, w,
571 $ iwork( indibl ), iwork( indisp ), z, ldz,
572 $ work( indwrk ), iwork( indiwo ), ifail, info )
573*
574* Apply orthogonal matrix used in reduction to tridiagonal
575* form to eigenvectors returned by SSTEIN.
576*
577 DO 20 j = 1, m
578 CALL scopy( n, z( 1, j ), 1, work( 1 ), 1 )
579 CALL sgemv( 'N', n, n, one, q, ldq, work, 1, zero,
580 $ z( 1, j ), 1 )
581 20 CONTINUE
582 END IF
583*
584* If matrix was scaled, then rescale eigenvalues appropriately.
585*
586 30 CONTINUE
587 IF( iscale.EQ.1 ) THEN
588 IF( info.EQ.0 ) THEN
589 imax = m
590 ELSE
591 imax = info - 1
592 END IF
593 CALL sscal( imax, one / sigma, w, 1 )
594 END IF
595*
596* If eigenvalues are not in order, then sort them, along with
597* eigenvectors.
598*
599 IF( wantz ) THEN
600 DO 50 j = 1, m - 1
601 i = 0
602 tmp1 = w( j )
603 DO 40 jj = j + 1, m
604 IF( w( jj ).LT.tmp1 ) THEN
605 i = jj
606 tmp1 = w( jj )
607 END IF
608 40 CONTINUE
609*
610 IF( i.NE.0 ) THEN
611 itmp1 = iwork( indibl+i-1 )
612 w( i ) = w( j )
613 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
614 w( j ) = tmp1
615 iwork( indibl+j-1 ) = itmp1
616 CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
617 IF( info.NE.0 ) THEN
618 itmp1 = ifail( i )
619 ifail( i ) = ifail( j )
620 ifail( j ) = itmp1
621 END IF
622 END IF
623 50 CONTINUE
624 END IF
625*
626* Set WORK(1) to optimal workspace size.
627*
628 work( 1 ) = sroundup_lwork(lwmin)
629*
630 RETURN
631*
632* End of SSBEVX_2STAGE
633*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:158
subroutine ssytrd_sb2st(stage1, vect, uplo, n, kd, ab, ldab, d, e, hous, lhous, work, lwork, info)
SSYTRD_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 slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
real function slansb(norm, uplo, n, k, ab, ldab, work)
SLANSB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slansb.f:129
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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
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 sstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
SSTEIN
Definition sstein.f:174
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 sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: