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

## ◆ ssyevx_2stage()

 subroutine ssyevx_2stage ( character jobz, character range, character uplo, integer n, real, dimension( lda, * ) a, integer lda, 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 )

SSYEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices

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

Purpose:
``` SSYEVX_2STAGE computes selected eigenvalues and, optionally, eigenvectors
of a real symmetric 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,out] A ``` A is REAL array, dimension (LDA, N) On entry, the symmetric matrix A. If UPLO = 'U', the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A. If UPLO = 'L', the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A. On exit, the lower triangle (if UPLO='L') or the upper triangle (if UPLO='U') of A, including the diagonal, is destroyed.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= 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 A 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) On normal exit, 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 (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal 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, 8*N, dimension) where dimension = max(stage1,stage2) + (KD+1)*N + 3*N = N*KD + N*max(KD+1,FACTOPTNB) + max(2*KD*KD, KD*NTHREADS) + (KD+1)*N + 3*N where KD is the blocking size of the reduction, FACTOPTNB is the blocking used by the QR or LQ algorithm, usually FACTOPTNB=128 is a good choice 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).
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 297 of file ssyevx_2stage.f.

300*
301 IMPLICIT NONE
302*
303* -- LAPACK driver routine --
304* -- LAPACK is a software package provided by Univ. of Tennessee, --
305* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
306*
307* .. Scalar Arguments ..
308 CHARACTER JOBZ, RANGE, UPLO
309 INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N
310 REAL ABSTOL, VL, VU
311* ..
312* .. Array Arguments ..
313 INTEGER IFAIL( * ), IWORK( * )
314 REAL A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
315* ..
316*
317* =====================================================================
318*
319* .. Parameters ..
320 REAL ZERO, ONE
321 parameter( zero = 0.0e+0, one = 1.0e+0 )
322* ..
323* .. Local Scalars ..
324 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
325 \$ WANTZ
326 CHARACTER ORDER
327 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
328 \$ INDISP, INDIWO, INDTAU, INDWKN, INDWRK, ISCALE,
329 \$ ITMP1, J, JJ, LLWORK, LLWRKN,
330 \$ NSPLIT, LWMIN, LHTRD, LWTRD, KD, IB, INDHOUS
331 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
332 \$ SIGMA, SMLNUM, TMP1, VLL, VUU
333* ..
334* .. External Functions ..
335 LOGICAL LSAME
336 INTEGER ILAENV2STAGE
337 REAL SLAMCH, SLANSY, SROUNDUP_LWORK
338 EXTERNAL lsame, slamch, slansy, ilaenv2stage,
340* ..
341* .. External Subroutines ..
342 EXTERNAL scopy, slacpy, sorgtr, sormtr, sscal, sstebz,
345* ..
346* .. Intrinsic Functions ..
347 INTRINSIC max, min, sqrt
348* ..
349* .. Executable Statements ..
350*
351* Test the input parameters.
352*
353 lower = lsame( uplo, 'L' )
354 wantz = lsame( jobz, 'V' )
355 alleig = lsame( range, 'A' )
356 valeig = lsame( range, 'V' )
357 indeig = lsame( range, 'I' )
358 lquery = ( lwork.EQ.-1 )
359*
360 info = 0
361 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
362 info = -1
363 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
364 info = -2
365 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
366 info = -3
367 ELSE IF( n.LT.0 ) THEN
368 info = -4
369 ELSE IF( lda.LT.max( 1, n ) ) THEN
370 info = -6
371 ELSE
372 IF( valeig ) THEN
373 IF( n.GT.0 .AND. vu.LE.vl )
374 \$ info = -8
375 ELSE IF( indeig ) THEN
376 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
377 info = -9
378 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
379 info = -10
380 END IF
381 END IF
382 END IF
383 IF( info.EQ.0 ) THEN
384 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
385 info = -15
386 END IF
387 END IF
388*
389 IF( info.EQ.0 ) THEN
390 IF( n.LE.1 ) THEN
391 lwmin = 1
392 work( 1 ) = sroundup_lwork(lwmin)
393 ELSE
394 kd = ilaenv2stage( 1, 'SSYTRD_2STAGE', jobz,
395 \$ n, -1, -1, -1 )
396 ib = ilaenv2stage( 2, 'SSYTRD_2STAGE', jobz,
397 \$ n, kd, -1, -1 )
398 lhtrd = ilaenv2stage( 3, 'SSYTRD_2STAGE', jobz,
399 \$ n, kd, ib, -1 )
400 lwtrd = ilaenv2stage( 4, 'SSYTRD_2STAGE', jobz,
401 \$ n, kd, ib, -1 )
402 lwmin = max( 8*n, 3*n + lhtrd + lwtrd )
403 work( 1 ) = lwmin
404 END IF
405*
406 IF( lwork.LT.lwmin .AND. .NOT.lquery )
407 \$ info = -17
408 END IF
409*
410 IF( info.NE.0 ) THEN
411 CALL xerbla( 'SSYEVX_2STAGE', -info )
412 RETURN
413 ELSE IF( lquery ) THEN
414 RETURN
415 END IF
416*
417* Quick return if possible
418*
419 m = 0
420 IF( n.EQ.0 ) THEN
421 RETURN
422 END IF
423*
424 IF( n.EQ.1 ) THEN
425 IF( alleig .OR. indeig ) THEN
426 m = 1
427 w( 1 ) = a( 1, 1 )
428 ELSE
429 IF( vl.LT.a( 1, 1 ) .AND. vu.GE.a( 1, 1 ) ) THEN
430 m = 1
431 w( 1 ) = a( 1, 1 )
432 END IF
433 END IF
434 IF( wantz )
435 \$ z( 1, 1 ) = one
436 RETURN
437 END IF
438*
439* Get machine constants.
440*
441 safmin = slamch( 'Safe minimum' )
442 eps = slamch( 'Precision' )
443 smlnum = safmin / eps
444 bignum = one / smlnum
445 rmin = sqrt( smlnum )
446 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
447*
448* Scale matrix to allowable range, if necessary.
449*
450 iscale = 0
451 abstll = abstol
452 IF( valeig ) THEN
453 vll = vl
454 vuu = vu
455 END IF
456 anrm = slansy( 'M', uplo, n, a, lda, work )
457 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
458 iscale = 1
459 sigma = rmin / anrm
460 ELSE IF( anrm.GT.rmax ) THEN
461 iscale = 1
462 sigma = rmax / anrm
463 END IF
464 IF( iscale.EQ.1 ) THEN
465 IF( lower ) THEN
466 DO 10 j = 1, n
467 CALL sscal( n-j+1, sigma, a( j, j ), 1 )
468 10 CONTINUE
469 ELSE
470 DO 20 j = 1, n
471 CALL sscal( j, sigma, a( 1, j ), 1 )
472 20 CONTINUE
473 END IF
474 IF( abstol.GT.0 )
475 \$ abstll = abstol*sigma
476 IF( valeig ) THEN
477 vll = vl*sigma
478 vuu = vu*sigma
479 END IF
480 END IF
481*
482* Call SSYTRD_2STAGE to reduce symmetric matrix to tridiagonal form.
483*
484 indtau = 1
485 inde = indtau + n
486 indd = inde + n
487 indhous = indd + n
488 indwrk = indhous + lhtrd
489 llwork = lwork - indwrk + 1
490*
491 CALL ssytrd_2stage( jobz, uplo, n, a, lda, work( indd ),
492 \$ work( inde ), work( indtau ), work( indhous ),
493 \$ lhtrd, work( indwrk ), llwork, iinfo )
494*
495* If all eigenvalues are desired and ABSTOL is less than or equal to
496* zero, then call SSTERF or SORGTR and SSTEQR. If this fails for
497* some eigenvalue, then try SSTEBZ.
498*
499 test = .false.
500 IF( indeig ) THEN
501 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
502 test = .true.
503 END IF
504 END IF
505 IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
506 CALL scopy( n, work( indd ), 1, w, 1 )
507 indee = indwrk + 2*n
508 IF( .NOT.wantz ) THEN
509 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
510 CALL ssterf( n, w, work( indee ), info )
511 ELSE
512 CALL slacpy( 'A', n, n, a, lda, z, ldz )
513 CALL sorgtr( uplo, n, z, ldz, work( indtau ),
514 \$ work( indwrk ), llwork, iinfo )
515 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
516 CALL ssteqr( jobz, n, w, work( indee ), z, ldz,
517 \$ work( indwrk ), info )
518 IF( info.EQ.0 ) THEN
519 DO 30 i = 1, n
520 ifail( i ) = 0
521 30 CONTINUE
522 END IF
523 END IF
524 IF( info.EQ.0 ) THEN
525 m = n
526 GO TO 40
527 END IF
528 info = 0
529 END IF
530*
531* Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
532*
533 IF( wantz ) THEN
534 order = 'B'
535 ELSE
536 order = 'E'
537 END IF
538 indibl = 1
539 indisp = indibl + n
540 indiwo = indisp + n
541 CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
542 \$ work( indd ), work( inde ), m, nsplit, w,
543 \$ iwork( indibl ), iwork( indisp ), work( indwrk ),
544 \$ iwork( indiwo ), info )
545*
546 IF( wantz ) THEN
547 CALL sstein( n, work( indd ), work( inde ), m, w,
548 \$ iwork( indibl ), iwork( indisp ), z, ldz,
549 \$ work( indwrk ), iwork( indiwo ), ifail, info )
550*
551* Apply orthogonal matrix used in reduction to tridiagonal
552* form to eigenvectors returned by SSTEIN.
553*
554 indwkn = inde
555 llwrkn = lwork - indwkn + 1
556 CALL sormtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
557 \$ ldz, work( indwkn ), llwrkn, iinfo )
558 END IF
559*
560* If matrix was scaled, then rescale eigenvalues appropriately.
561*
562 40 CONTINUE
563 IF( iscale.EQ.1 ) THEN
564 IF( info.EQ.0 ) THEN
565 imax = m
566 ELSE
567 imax = info - 1
568 END IF
569 CALL sscal( imax, one / sigma, w, 1 )
570 END IF
571*
572* If eigenvalues are not in order, then sort them, along with
573* eigenvectors.
574*
575 IF( wantz ) THEN
576 DO 60 j = 1, m - 1
577 i = 0
578 tmp1 = w( j )
579 DO 50 jj = j + 1, m
580 IF( w( jj ).LT.tmp1 ) THEN
581 i = jj
582 tmp1 = w( jj )
583 END IF
584 50 CONTINUE
585*
586 IF( i.NE.0 ) THEN
587 itmp1 = iwork( indibl+i-1 )
588 w( i ) = w( j )
589 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
590 w( j ) = tmp1
591 iwork( indibl+j-1 ) = itmp1
592 CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
593 IF( info.NE.0 ) THEN
594 itmp1 = ifail( i )
595 ifail( i ) = ifail( j )
596 ifail( j ) = itmp1
597 END IF
598 END IF
599 60 CONTINUE
600 END IF
601*
602* Set WORK(1) to optimal workspace size.
603*
604 work( 1 ) = sroundup_lwork(lwmin)
605*
606 RETURN
607*
608* End of SSYEVX_2STAGE
609*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine ssytrd_2stage(vect, uplo, n, a, lda, d, e, tau, hous2, lhous2, work, lwork, info)
SSYTRD_2STAGE
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 slansy(norm, uplo, n, a, lda, work)
SLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slansy.f:122
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
subroutine sorgtr(uplo, n, a, lda, tau, work, lwork, info)
SORGTR
Definition sorgtr.f:123
subroutine sormtr(side, uplo, trans, m, n, a, lda, tau, c, ldc, work, lwork, info)
SORMTR
Definition sormtr.f:172
Here is the call graph for this function:
Here is the caller graph for this function: