OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
real
LAPACK » Eigenvalue

Functions

subroutine sbdsvdx (uplo, jobz, range, n, d, e, vl, vu, il, iu, ns, s, z, ldz, work, iwork, info)
 SBDSVDX
subroutine sggglm (n, m, p, a, lda, b, ldb, d, x, y, work, lwork, info)
 SGGGLM
subroutine ssbev (jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, info)
  SSBEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine ssbev_2stage (jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, lwork, info)
  SSBEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine ssbevd (jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, lwork, iwork, liwork, info)
  SSBEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine ssbevd_2stage (jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, lwork, iwork, liwork, info)
  SSBEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine ssbevx (jobz, range, uplo, n, kd, ab, ldab, q, ldq, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info)
  SSBEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine ssbevx_2stage (jobz, range, uplo, n, kd, ab, ldab, q, ldq, vl, vu, il, iu, abstol, m, w, z, ldz, work, lwork, iwork, ifail, info)
  SSBEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine ssbgv (jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, w, z, ldz, work, info)
 SSBGV
subroutine ssbgvd (jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, w, z, ldz, work, lwork, iwork, liwork, info)
 SSBGVD
subroutine ssbgvx (jobz, range, uplo, n, ka, kb, ab, ldab, bb, ldbb, q, ldq, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info)
 SSBGVX
subroutine sspev (jobz, uplo, n, ap, w, z, ldz, work, info)
  SSPEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine sspevd (jobz, uplo, n, ap, w, z, ldz, work, lwork, iwork, liwork, info)
  SSPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine sspevx (jobz, range, uplo, n, ap, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info)
  SSPEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine sspgv (itype, jobz, uplo, n, ap, bp, w, z, ldz, work, info)
 SSPGV
subroutine sspgvd (itype, jobz, uplo, n, ap, bp, w, z, ldz, work, lwork, iwork, liwork, info)
 SSPGVD
subroutine sspgvx (itype, jobz, range, uplo, n, ap, bp, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info)
 SSPGVX
subroutine sstev (jobz, n, d, e, z, ldz, work, info)
  SSTEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine sstevd (jobz, n, d, e, z, ldz, work, lwork, iwork, liwork, info)
  SSTEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine sstevr (jobz, range, n, d, e, vl, vu, il, iu, abstol, m, w, z, ldz, isuppz, work, lwork, iwork, liwork, info)
  SSTEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine sstevx (jobz, range, n, d, e, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info)
  SSTEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices

Detailed Description

This is the group of real Other Eigenvalue routines

Function Documentation

◆ sbdsvdx()

subroutine sbdsvdx ( character uplo,
character jobz,
character range,
integer n,
real, dimension( * ) d,
real, dimension( * ) e,
real vl,
real vu,
integer il,
integer iu,
integer ns,
real, dimension( * ) s,
real, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer, dimension( * ) iwork,
integer info )

SBDSVDX

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

Purpose:
!>
!>  SBDSVDX computes the singular value decomposition (SVD) of a real
!>  N-by-N (upper or lower) bidiagonal matrix B, B = U * S * VT,
!>  where S is a diagonal matrix with non-negative diagonal elements
!>  (the singular values of B), and U and VT are orthogonal matrices
!>  of left and right singular vectors, respectively.
!>
!>  Given an upper bidiagonal B with diagonal D = [ d_1 d_2 ... d_N ]
!>  and superdiagonal E = [ e_1 e_2 ... e_N-1 ], SBDSVDX computes the
!>  singular value decompositon of B through the eigenvalues and
!>  eigenvectors of the N*2-by-N*2 tridiagonal matrix
!>
!>        |  0  d_1                |
!>        | d_1  0  e_1            |
!>  TGK = |     e_1  0  d_2        |
!>        |         d_2  .   .     |
!>        |              .   .   . |
!>
!>  If (s,u,v) is a singular triplet of B with ||u|| = ||v|| = 1, then
!>  (+/-s,q), ||q|| = 1, are eigenpairs of TGK, with q = P * ( u' +/-v' ) /
!>  sqrt(2) = ( v_1 u_1 v_2 u_2 ... v_n u_n ) / sqrt(2), and
!>  P = [ e_{n+1} e_{1} e_{n+2} e_{2} ... ].
!>
!>  Given a TGK matrix, one can either a) compute -s,-v and change signs
!>  so that the singular values (and corresponding vectors) are already in
!>  descending order (as in SGESVD/SGESDD) or b) compute s,v and reorder
!>  the values (and corresponding vectors). SBDSVDX implements a) by
!>  calling SSTEVX (bisection plus inverse iteration, to be replaced
!>  with a version of the Multiple Relative Robust Representation
!>  algorithm. (See P. Willems and B. Lang, A framework for the MR^3
!>  algorithm: theory and implementation, SIAM J. Sci. Comput.,
!>  35:740-766, 2013.)
!>
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  B is upper bidiagonal;
!>          = 'L':  B is lower bidiagonal.
!>
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute singular values only;
!>          = 'V':  Compute singular values and singular vectors.
!>
[in]RANGE
!>          RANGE is CHARACTER*1
!>          = 'A': all singular values will be found.
!>          = 'V': all singular values in the half-open interval [VL,VU)
!>                 will be found.
!>          = 'I': the IL-th through IU-th singular values will be found.
!>
[in]N
!>          N is INTEGER
!>          The order of the bidiagonal matrix.  N >= 0.
!>
[in]D
!>          D is REAL array, dimension (N)
!>          The n diagonal elements of the bidiagonal matrix B.
!>
[in]E
!>          E is REAL array, dimension (max(1,N-1))
!>          The (n-1) superdiagonal elements of the bidiagonal matrix
!>          B in elements 1 to N-1.
!>
[in]VL
!>         VL is REAL
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for singular values. VU > VL.
!>          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 singular values. VU > VL.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest singular value to be returned.
!>          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest singular value to be returned.
!>          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[out]NS
!>          NS is INTEGER
!>          The total number of singular values found.  0 <= NS <= N.
!>          If RANGE = 'A', NS = N, and if RANGE = 'I', NS = IU-IL+1.
!>
[out]S
!>          S is REAL array, dimension (N)
!>          The first NS elements contain the selected singular values in
!>          ascending order.
!>
[out]Z
!>          Z is REAL array, dimension (2*N,K)
!>          If JOBZ = 'V', then if INFO = 0 the first NS columns of Z
!>          contain the singular vectors of the matrix B corresponding to
!>          the selected singular values, with U in rows 1 to N and V
!>          in rows N+1 to N*2, i.e.
!>          Z = [ U ]
!>              [ V ]
!>          If JOBZ = 'N', then Z is not referenced.
!>          Note: The user must ensure that at least K = NS+1 columns are
!>          supplied in the array Z; if RANGE = 'V', the exact value of
!>          NS 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(2,N*2).
!>
[out]WORK
!>          WORK is REAL array, dimension (14*N)
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (12*N)
!>          If JOBZ = 'V', then if INFO = 0, the first NS elements of
!>          IWORK are zero. If INFO > 0, then IWORK contains the indices
!>          of the eigenvectors that failed to converge in DSTEVX.
!>
[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
!>                   in SSTEVX. The indices of the eigenvectors
!>                   (as returned by SSTEVX) are stored in the
!>                   array IWORK.
!>                if INFO = N*2 + 1, an internal error occurred.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 224 of file sbdsvdx.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, UPLO
233 INTEGER IL, INFO, IU, LDZ, N, NS
234 REAL VL, VU
235* ..
236* .. Array Arguments ..
237 INTEGER IWORK( * )
238 REAL D( * ), E( * ), S( * ), WORK( * ),
239 $ Z( LDZ, * )
240* ..
241*
242* =====================================================================
243*
244* .. Parameters ..
245 REAL ZERO, ONE, TEN, HNDRD, MEIGTH
246 parameter( zero = 0.0e0, one = 1.0e0, ten = 10.0e0,
247 $ hndrd = 100.0e0, meigth = -0.1250e0 )
248 REAL FUDGE
249 parameter( fudge = 2.0e0 )
250* ..
251* .. Local Scalars ..
252 CHARACTER RNGVX
253 LOGICAL ALLSV, INDSV, LOWER, SPLIT, SVEQ0, VALSV, WANTZ
254 INTEGER I, ICOLZ, IDBEG, IDEND, IDTGK, IDPTR, IEPTR,
255 $ IETGK, IIFAIL, IIWORK, ILTGK, IROWU, IROWV,
256 $ IROWZ, ISBEG, ISPLT, ITEMP, IUTGK, J, K,
257 $ NTGK, NRU, NRV, NSL
258 REAL ABSTOL, EPS, EMIN, MU, NRMU, NRMV, ORTOL, SMAX,
259 $ SMIN, SQRT2, THRESH, TOL, ULP,
260 $ VLTGK, VUTGK, ZJTJI
261* ..
262* .. External Functions ..
263 LOGICAL LSAME
264 INTEGER ISAMAX
265 REAL SDOT, SLAMCH, SNRM2
266 EXTERNAL isamax, lsame, saxpy, sdot, slamch, snrm2
267* ..
268* .. External Subroutines ..
269 EXTERNAL scopy, slaset, sscal, sswap, sstevx, xerbla
270* ..
271* .. Intrinsic Functions ..
272 INTRINSIC abs, real, sign, sqrt
273* ..
274* .. Executable Statements ..
275*
276* Test the input parameters.
277*
278 allsv = lsame( range, 'A' )
279 valsv = lsame( range, 'V' )
280 indsv = lsame( range, 'I' )
281 wantz = lsame( jobz, 'V' )
282 lower = lsame( uplo, 'L' )
283*
284 info = 0
285 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lower ) THEN
286 info = -1
287 ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
288 info = -2
289 ELSE IF( .NOT.( allsv .OR. valsv .OR. indsv ) ) THEN
290 info = -3
291 ELSE IF( n.LT.0 ) THEN
292 info = -4
293 ELSE IF( n.GT.0 ) THEN
294 IF( valsv ) THEN
295 IF( vl.LT.zero ) THEN
296 info = -7
297 ELSE IF( vu.LE.vl ) THEN
298 info = -8
299 END IF
300 ELSE IF( indsv ) THEN
301 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
302 info = -9
303 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
304 info = -10
305 END IF
306 END IF
307 END IF
308 IF( info.EQ.0 ) THEN
309 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n*2 ) ) info = -14
310 END IF
311*
312 IF( info.NE.0 ) THEN
313 CALL xerbla( 'SBDSVDX', -info )
314 RETURN
315 END IF
316*
317* Quick return if possible (N.LE.1)
318*
319 ns = 0
320 IF( n.EQ.0 ) RETURN
321*
322 IF( n.EQ.1 ) THEN
323 IF( allsv .OR. indsv ) THEN
324 ns = 1
325 s( 1 ) = abs( d( 1 ) )
326 ELSE
327 IF( vl.LT.abs( d( 1 ) ) .AND. vu.GE.abs( d( 1 ) ) ) THEN
328 ns = 1
329 s( 1 ) = abs( d( 1 ) )
330 END IF
331 END IF
332 IF( wantz ) THEN
333 z( 1, 1 ) = sign( one, d( 1 ) )
334 z( 2, 1 ) = one
335 END IF
336 RETURN
337 END IF
338*
339 abstol = 2*slamch( 'Safe Minimum' )
340 ulp = slamch( 'Precision' )
341 eps = slamch( 'Epsilon' )
342 sqrt2 = sqrt( 2.0e0 )
343 ortol = sqrt( ulp )
344*
345* Criterion for splitting is taken from SBDSQR when singular
346* values are computed to relative accuracy TOL. (See J. Demmel and
347* W. Kahan, Accurate singular values of bidiagonal matrices, SIAM
348* J. Sci. and Stat. Comput., 11:873–912, 1990.)
349*
350 tol = max( ten, min( hndrd, eps**meigth ) )*eps
351*
352* Compute approximate maximum, minimum singular values.
353*
354 i = isamax( n, d, 1 )
355 smax = abs( d( i ) )
356 i = isamax( n-1, e, 1 )
357 smax = max( smax, abs( e( i ) ) )
358*
359* Compute threshold for neglecting D's and E's.
360*
361 smin = abs( d( 1 ) )
362 IF( smin.NE.zero ) THEN
363 mu = smin
364 DO i = 2, n
365 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
366 smin = min( smin, mu )
367 IF( smin.EQ.zero ) EXIT
368 END DO
369 END IF
370 smin = smin / sqrt( real( n ) )
371 thresh = tol*smin
372*
373* Check for zeros in D and E (splits), i.e. submatrices.
374*
375 DO i = 1, n-1
376 IF( abs( d( i ) ).LE.thresh ) d( i ) = zero
377 IF( abs( e( i ) ).LE.thresh ) e( i ) = zero
378 END DO
379 IF( abs( d( n ) ).LE.thresh ) d( n ) = zero
380*
381* Pointers for arrays used by SSTEVX.
382*
383 idtgk = 1
384 ietgk = idtgk + n*2
385 itemp = ietgk + n*2
386 iifail = 1
387 iiwork = iifail + n*2
388*
389* Set RNGVX, which corresponds to RANGE for SSTEVX in TGK mode.
390* VL,VU or IL,IU are redefined to conform to implementation a)
391* described in the leading comments.
392*
393 iltgk = 0
394 iutgk = 0
395 vltgk = zero
396 vutgk = zero
397*
398 IF( allsv ) THEN
399*
400* All singular values will be found. We aim at -s (see
401* leading comments) with RNGVX = 'I'. IL and IU are set
402* later (as ILTGK and IUTGK) according to the dimension
403* of the active submatrix.
404*
405 rngvx = 'I'
406 IF( wantz ) CALL slaset( 'F', n*2, n+1, zero, zero, z, ldz )
407 ELSE IF( valsv ) THEN
408*
409* Find singular values in a half-open interval. We aim
410* at -s (see leading comments) and we swap VL and VU
411* (as VUTGK and VLTGK), changing their signs.
412*
413 rngvx = 'V'
414 vltgk = -vu
415 vutgk = -vl
416 work( idtgk:idtgk+2*n-1 ) = zero
417 CALL scopy( n, d, 1, work( ietgk ), 2 )
418 CALL scopy( n-1, e, 1, work( ietgk+1 ), 2 )
419 CALL sstevx( 'N', 'V', n*2, work( idtgk ), work( ietgk ),
420 $ vltgk, vutgk, iltgk, iltgk, abstol, ns, s,
421 $ z, ldz, work( itemp ), iwork( iiwork ),
422 $ iwork( iifail ), info )
423 IF( ns.EQ.0 ) THEN
424 RETURN
425 ELSE
426 IF( wantz ) CALL slaset( 'F', n*2, ns, zero, zero, z, ldz )
427 END IF
428 ELSE IF( indsv ) THEN
429*
430* Find the IL-th through the IU-th singular values. We aim
431* at -s (see leading comments) and indices are mapped into
432* values, therefore mimicking SSTEBZ, where
433*
434* GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN
435* GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN
436*
437 iltgk = il
438 iutgk = iu
439 rngvx = 'V'
440 work( idtgk:idtgk+2*n-1 ) = zero
441 CALL scopy( n, d, 1, work( ietgk ), 2 )
442 CALL scopy( n-1, e, 1, work( ietgk+1 ), 2 )
443 CALL sstevx( 'N', 'I', n*2, work( idtgk ), work( ietgk ),
444 $ vltgk, vltgk, iltgk, iltgk, abstol, ns, s,
445 $ z, ldz, work( itemp ), iwork( iiwork ),
446 $ iwork( iifail ), info )
447 vltgk = s( 1 ) - fudge*smax*ulp*n
448 work( idtgk:idtgk+2*n-1 ) = zero
449 CALL scopy( n, d, 1, work( ietgk ), 2 )
450 CALL scopy( n-1, e, 1, work( ietgk+1 ), 2 )
451 CALL sstevx( 'N', 'I', n*2, work( idtgk ), work( ietgk ),
452 $ vutgk, vutgk, iutgk, iutgk, abstol, ns, s,
453 $ z, ldz, work( itemp ), iwork( iiwork ),
454 $ iwork( iifail ), info )
455 vutgk = s( 1 ) + fudge*smax*ulp*n
456 vutgk = min( vutgk, zero )
457*
458* If VLTGK=VUTGK, SSTEVX returns an error message,
459* so if needed we change VUTGK slightly.
460*
461 IF( vltgk.EQ.vutgk ) vltgk = vltgk - tol
462*
463 IF( wantz ) CALL slaset( 'F', n*2, iu-il+1, zero, zero, z, ldz)
464 END IF
465*
466* Initialize variables and pointers for S, Z, and WORK.
467*
468* NRU, NRV: number of rows in U and V for the active submatrix
469* IDBEG, ISBEG: offsets for the entries of D and S
470* IROWZ, ICOLZ: offsets for the rows and columns of Z
471* IROWU, IROWV: offsets for the rows of U and V
472*
473 ns = 0
474 nru = 0
475 nrv = 0
476 idbeg = 1
477 isbeg = 1
478 irowz = 1
479 icolz = 1
480 irowu = 2
481 irowv = 1
482 split = .false.
483 sveq0 = .false.
484*
485* Form the tridiagonal TGK matrix.
486*
487 s( 1:n ) = zero
488 work( ietgk+2*n-1 ) = zero
489 work( idtgk:idtgk+2*n-1 ) = zero
490 CALL scopy( n, d, 1, work( ietgk ), 2 )
491 CALL scopy( n-1, e, 1, work( ietgk+1 ), 2 )
492*
493*
494* Check for splits in two levels, outer level
495* in E and inner level in D.
496*
497 DO ieptr = 2, n*2, 2
498 IF( work( ietgk+ieptr-1 ).EQ.zero ) THEN
499*
500* Split in E (this piece of B is square) or bottom
501* of the (input bidiagonal) matrix.
502*
503 isplt = idbeg
504 idend = ieptr - 1
505 DO idptr = idbeg, idend, 2
506 IF( work( ietgk+idptr-1 ).EQ.zero ) THEN
507*
508* Split in D (rectangular submatrix). Set the number
509* of rows in U and V (NRU and NRV) accordingly.
510*
511 IF( idptr.EQ.idbeg ) THEN
512*
513* D=0 at the top.
514*
515 sveq0 = .true.
516 IF( idbeg.EQ.idend) THEN
517 nru = 1
518 nrv = 1
519 END IF
520 ELSE IF( idptr.EQ.idend ) THEN
521*
522* D=0 at the bottom.
523*
524 sveq0 = .true.
525 nru = (idend-isplt)/2 + 1
526 nrv = nru
527 IF( isplt.NE.idbeg ) THEN
528 nru = nru + 1
529 END IF
530 ELSE
531 IF( isplt.EQ.idbeg ) THEN
532*
533* Split: top rectangular submatrix.
534*
535 nru = (idptr-idbeg)/2
536 nrv = nru + 1
537 ELSE
538*
539* Split: middle square submatrix.
540*
541 nru = (idptr-isplt)/2 + 1
542 nrv = nru
543 END IF
544 END IF
545 ELSE IF( idptr.EQ.idend ) THEN
546*
547* Last entry of D in the active submatrix.
548*
549 IF( isplt.EQ.idbeg ) THEN
550*
551* No split (trivial case).
552*
553 nru = (idend-idbeg)/2 + 1
554 nrv = nru
555 ELSE
556*
557* Split: bottom rectangular submatrix.
558*
559 nrv = (idend-isplt)/2 + 1
560 nru = nrv + 1
561 END IF
562 END IF
563*
564 ntgk = nru + nrv
565*
566 IF( ntgk.GT.0 ) THEN
567*
568* Compute eigenvalues/vectors of the active
569* submatrix according to RANGE:
570* if RANGE='A' (ALLSV) then RNGVX = 'I'
571* if RANGE='V' (VALSV) then RNGVX = 'V'
572* if RANGE='I' (INDSV) then RNGVX = 'V'
573*
574 iltgk = 1
575 iutgk = ntgk / 2
576 IF( allsv .OR. vutgk.EQ.zero ) THEN
577 IF( sveq0 .OR.
578 $ smin.LT.eps .OR.
579 $ mod(ntgk,2).GT.0 ) THEN
580* Special case: eigenvalue equal to zero or very
581* small, additional eigenvector is needed.
582 iutgk = iutgk + 1
583 END IF
584 END IF
585*
586* Workspace needed by SSTEVX:
587* WORK( ITEMP: ): 2*5*NTGK
588* IWORK( 1: ): 2*6*NTGK
589*
590 CALL sstevx( jobz, rngvx, ntgk, work( idtgk+isplt-1 ),
591 $ work( ietgk+isplt-1 ), vltgk, vutgk,
592 $ iltgk, iutgk, abstol, nsl, s( isbeg ),
593 $ z( irowz,icolz ), ldz, work( itemp ),
594 $ iwork( iiwork ), iwork( iifail ),
595 $ info )
596 IF( info.NE.0 ) THEN
597* Exit with the error code from SSTEVX.
598 RETURN
599 END IF
600 emin = abs( maxval( s( isbeg:isbeg+nsl-1 ) ) )
601*
602 IF( nsl.GT.0 .AND. wantz ) THEN
603*
604* Normalize u=Z([2,4,...],:) and v=Z([1,3,...],:),
605* changing the sign of v as discussed in the leading
606* comments. The norms of u and v may be (slightly)
607* different from 1/sqrt(2) if the corresponding
608* eigenvalues are very small or too close. We check
609* those norms and, if needed, reorthogonalize the
610* vectors.
611*
612 IF( nsl.GT.1 .AND.
613 $ vutgk.EQ.zero .AND.
614 $ mod(ntgk,2).EQ.0 .AND.
615 $ emin.EQ.0 .AND. .NOT.split ) THEN
616*
617* D=0 at the top or bottom of the active submatrix:
618* one eigenvalue is equal to zero; concatenate the
619* eigenvectors corresponding to the two smallest
620* eigenvalues.
621*
622 z( irowz:irowz+ntgk-1,icolz+nsl-2 ) =
623 $ z( irowz:irowz+ntgk-1,icolz+nsl-2 ) +
624 $ z( irowz:irowz+ntgk-1,icolz+nsl-1 )
625 z( irowz:irowz+ntgk-1,icolz+nsl-1 ) =
626 $ zero
627* IF( IUTGK*2.GT.NTGK ) THEN
628* Eigenvalue equal to zero or very small.
629* NSL = NSL - 1
630* END IF
631 END IF
632*
633 DO i = 0, min( nsl-1, nru-1 )
634 nrmu = snrm2( nru, z( irowu, icolz+i ), 2 )
635 IF( nrmu.EQ.zero ) THEN
636 info = n*2 + 1
637 RETURN
638 END IF
639 CALL sscal( nru, one/nrmu,
640 $ z( irowu,icolz+i ), 2 )
641 IF( nrmu.NE.one .AND.
642 $ abs( nrmu-ortol )*sqrt2.GT.one )
643 $ THEN
644 DO j = 0, i-1
645 zjtji = -sdot( nru, z( irowu, icolz+j ),
646 $ 2, z( irowu, icolz+i ), 2 )
647 CALL saxpy( nru, zjtji,
648 $ z( irowu, icolz+j ), 2,
649 $ z( irowu, icolz+i ), 2 )
650 END DO
651 nrmu = snrm2( nru, z( irowu, icolz+i ), 2 )
652 CALL sscal( nru, one/nrmu,
653 $ z( irowu,icolz+i ), 2 )
654 END IF
655 END DO
656 DO i = 0, min( nsl-1, nrv-1 )
657 nrmv = snrm2( nrv, z( irowv, icolz+i ), 2 )
658 IF( nrmv.EQ.zero ) THEN
659 info = n*2 + 1
660 RETURN
661 END IF
662 CALL sscal( nrv, -one/nrmv,
663 $ z( irowv,icolz+i ), 2 )
664 IF( nrmv.NE.one .AND.
665 $ abs( nrmv-ortol )*sqrt2.GT.one )
666 $ THEN
667 DO j = 0, i-1
668 zjtji = -sdot( nrv, z( irowv, icolz+j ),
669 $ 2, z( irowv, icolz+i ), 2 )
670 CALL saxpy( nru, zjtji,
671 $ z( irowv, icolz+j ), 2,
672 $ z( irowv, icolz+i ), 2 )
673 END DO
674 nrmv = snrm2( nrv, z( irowv, icolz+i ), 2 )
675 CALL sscal( nrv, one/nrmv,
676 $ z( irowv,icolz+i ), 2 )
677 END IF
678 END DO
679 IF( vutgk.EQ.zero .AND.
680 $ idptr.LT.idend .AND.
681 $ mod(ntgk,2).GT.0 ) THEN
682*
683* D=0 in the middle of the active submatrix (one
684* eigenvalue is equal to zero): save the corresponding
685* eigenvector for later use (when bottom of the
686* active submatrix is reached).
687*
688 split = .true.
689 z( irowz:irowz+ntgk-1,n+1 ) =
690 $ z( irowz:irowz+ntgk-1,ns+nsl )
691 z( irowz:irowz+ntgk-1,ns+nsl ) =
692 $ zero
693 END IF
694 END IF !** WANTZ **!
695*
696 nsl = min( nsl, nru )
697 sveq0 = .false.
698*
699* Absolute values of the eigenvalues of TGK.
700*
701 DO i = 0, nsl-1
702 s( isbeg+i ) = abs( s( isbeg+i ) )
703 END DO
704*
705* Update pointers for TGK, S and Z.
706*
707 isbeg = isbeg + nsl
708 irowz = irowz + ntgk
709 icolz = icolz + nsl
710 irowu = irowz
711 irowv = irowz + 1
712 isplt = idptr + 1
713 ns = ns + nsl
714 nru = 0
715 nrv = 0
716 END IF !** NTGK.GT.0 **!
717 IF( irowz.LT.n*2 .AND. wantz ) THEN
718 z( 1:irowz-1, icolz ) = zero
719 END IF
720 END DO !** IDPTR loop **!
721 IF( split .AND. wantz ) THEN
722*
723* Bring back eigenvector corresponding
724* to eigenvalue equal to zero.
725*
726 z( idbeg:idend-ntgk+1,isbeg-1 ) =
727 $ z( idbeg:idend-ntgk+1,isbeg-1 ) +
728 $ z( idbeg:idend-ntgk+1,n+1 )
729 z( idbeg:idend-ntgk+1,n+1 ) = 0
730 END IF
731 irowv = irowv - 1
732 irowu = irowu + 1
733 idbeg = ieptr + 1
734 sveq0 = .false.
735 split = .false.
736 END IF !** Check for split in E **!
737 END DO !** IEPTR loop **!
738*
739* Sort the singular values into decreasing order (insertion sort on
740* singular values, but only one transposition per singular vector)
741*
742 DO i = 1, ns-1
743 k = 1
744 smin = s( 1 )
745 DO j = 2, ns + 1 - i
746 IF( s( j ).LE.smin ) THEN
747 k = j
748 smin = s( j )
749 END IF
750 END DO
751 IF( k.NE.ns+1-i ) THEN
752 s( k ) = s( ns+1-i )
753 s( ns+1-i ) = smin
754 IF( wantz ) CALL sswap( n*2, z( 1,k ), 1, z( 1,ns+1-i ), 1 )
755 END IF
756 END DO
757*
758* If RANGE=I, check for singular values/vectors to be discarded.
759*
760 IF( indsv ) THEN
761 k = iu - il + 1
762 IF( k.LT.ns ) THEN
763 s( k+1:ns ) = zero
764 IF( wantz ) z( 1:n*2,k+1:ns ) = zero
765 ns = k
766 END IF
767 END IF
768*
769* Reorder Z: U = Z( 1:N,1:NS ), V = Z( N+1:N*2,1:NS ).
770* If B is a lower diagonal, swap U and V.
771*
772 IF( wantz ) THEN
773 DO i = 1, ns
774 CALL scopy( n*2, z( 1,i ), 1, work, 1 )
775 IF( lower ) THEN
776 CALL scopy( n, work( 2 ), 2, z( n+1,i ), 1 )
777 CALL scopy( n, work( 1 ), 2, z( 1 ,i ), 1 )
778 ELSE
779 CALL scopy( n, work( 2 ), 2, z( 1 ,i ), 1 )
780 CALL scopy( n, work( 1 ), 2, z( n+1,i ), 1 )
781 END IF
782 END DO
783 END IF
784*
785 RETURN
786*
787* End of SBDSVDX
788*
slaset
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:110
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
sstevx
subroutine sstevx(jobz, range, n, d, e, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info)
SSTEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition sstevx.f:227
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
sdot
real function sdot(n, sx, incx, sy, incy)
SDOT
Definition sdot.f:82
saxpy
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
snrm2
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
split
void split(mapping_t *, PORD_INT, PORD_INT, PORD_INT, PORD_INT *, PORD_INT *, FLOAT *, PORD_INT)

◆ sggglm()

subroutine sggglm ( integer n,
integer m,
integer p,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) d,
real, dimension( * ) x,
real, dimension( * ) y,
real, dimension( * ) work,
integer lwork,
integer info )

SGGGLM

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

Purpose:
!>
!> SGGGLM solves a general Gauss-Markov linear model (GLM) problem:
!>
!>         minimize || y ||_2   subject to   d = A*x + B*y
!>             x
!>
!> where A is an N-by-M matrix, B is an N-by-P matrix, and d is a
!> given N-vector. It is assumed that M <= N <= M+P, and
!>
!>            rank(A) = M    and    rank( A B ) = N.
!>
!> Under these assumptions, the constrained equation is always
!> consistent, and there is a unique solution x and a minimal 2-norm
!> solution y, which is obtained using a generalized QR factorization
!> of the matrices (A, B) given by
!>
!>    A = Q*(R),   B = Q*T*Z.
!>          (0)
!>
!> In particular, if matrix B is square nonsingular, then the problem
!> GLM is equivalent to the following weighted linear least squares
!> problem
!>
!>              minimize || inv(B)*(d-A*x) ||_2
!>                  x
!>
!> where inv(B) denotes the inverse of B.
!>
Parameters
[in]N
!>          N is INTEGER
!>          The number of rows of the matrices A and B.  N >= 0.
!>
[in]M
!>          M is INTEGER
!>          The number of columns of the matrix A.  0 <= M <= N.
!>
[in]P
!>          P is INTEGER
!>          The number of columns of the matrix B.  P >= N-M.
!>
[in,out]A
!>          A is REAL array, dimension (LDA,M)
!>          On entry, the N-by-M matrix A.
!>          On exit, the upper triangular part of the array A contains
!>          the M-by-M upper triangular matrix R.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A. LDA >= max(1,N).
!>
[in,out]B
!>          B is REAL array, dimension (LDB,P)
!>          On entry, the N-by-P matrix B.
!>          On exit, if N <= P, the upper triangle of the subarray
!>          B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T;
!>          if N > P, the elements on and above the (N-P)th subdiagonal
!>          contain the N-by-P upper trapezoidal matrix T.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B. LDB >= max(1,N).
!>
[in,out]D
!>          D is REAL array, dimension (N)
!>          On entry, D is the left hand side of the GLM equation.
!>          On exit, D is destroyed.
!>
[out]X
!>          X is REAL array, dimension (M)
!>
[out]Y
!>          Y is REAL array, dimension (P)
!>
!>          On exit, X and Y are the solutions of the GLM problem.
!>
[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 dimension of the array WORK. LWORK >= max(1,N+M+P).
!>          For optimum performance, LWORK >= M+min(N,P)+max(N,P)*NB,
!>          where NB is an upper bound for the optimal blocksizes for
!>          SGEQRF, SGERQF, SORMQR and SORMRQ.
!>
!>          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]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          = 1:  the upper triangular factor R associated with A in the
!>                generalized QR factorization of the pair (A, B) is
!>                singular, so that rank(A) < M; the least squares
!>                solution could not be computed.
!>          = 2:  the bottom (N-M) by (N-M) part of the upper trapezoidal
!>                factor T associated with B in the generalized QR
!>                factorization of the pair (A, B) is singular, so that
!>                rank( A B ) < N; the least squares solution could not
!>                be computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 183 of file sggglm.f.

185*
186* -- LAPACK driver routine --
187* -- LAPACK is a software package provided by Univ. of Tennessee, --
188* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
189*
190* .. Scalar Arguments ..
191 INTEGER INFO, LDA, LDB, LWORK, M, N, P
192* ..
193* .. Array Arguments ..
194 REAL A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
195 $ X( * ), Y( * )
196* ..
197*
198* ===================================================================
199*
200* .. Parameters ..
201 REAL ZERO, ONE
202 parameter( zero = 0.0e+0, one = 1.0e+0 )
203* ..
204* .. Local Scalars ..
205 LOGICAL LQUERY
206 INTEGER I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3,
207 $ NB4, NP
208* ..
209* .. External Subroutines ..
210 EXTERNAL scopy, sgemv, sggqrf, sormqr, sormrq, strtrs,
211 $ xerbla
212* ..
213* .. External Functions ..
214 INTEGER ILAENV
215 EXTERNAL ilaenv
216* ..
217* .. Intrinsic Functions ..
218 INTRINSIC int, max, min
219* ..
220* .. Executable Statements ..
221*
222* Test the input parameters
223*
224 info = 0
225 np = min( n, p )
226 lquery = ( lwork.EQ.-1 )
227 IF( n.LT.0 ) THEN
228 info = -1
229 ELSE IF( m.LT.0 .OR. m.GT.n ) THEN
230 info = -2
231 ELSE IF( p.LT.0 .OR. p.LT.n-m ) THEN
232 info = -3
233 ELSE IF( lda.LT.max( 1, n ) ) THEN
234 info = -5
235 ELSE IF( ldb.LT.max( 1, n ) ) THEN
236 info = -7
237 END IF
238*
239* Calculate workspace
240*
241 IF( info.EQ.0) THEN
242 IF( n.EQ.0 ) THEN
243 lwkmin = 1
244 lwkopt = 1
245 ELSE
246 nb1 = ilaenv( 1, 'SGEQRF', ' ', n, m, -1, -1 )
247 nb2 = ilaenv( 1, 'SGERQF', ' ', n, m, -1, -1 )
248 nb3 = ilaenv( 1, 'SORMQR', ' ', n, m, p, -1 )
249 nb4 = ilaenv( 1, 'SORMRQ', ' ', n, m, p, -1 )
250 nb = max( nb1, nb2, nb3, nb4 )
251 lwkmin = m + n + p
252 lwkopt = m + np + max( n, p )*nb
253 END IF
254 work( 1 ) = lwkopt
255*
256 IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
257 info = -12
258 END IF
259 END IF
260*
261 IF( info.NE.0 ) THEN
262 CALL xerbla( 'SGGGLM', -info )
263 RETURN
264 ELSE IF( lquery ) THEN
265 RETURN
266 END IF
267*
268* Quick return if possible
269*
270 IF( n.EQ.0 ) THEN
271 DO i = 1, m
272 x(i) = zero
273 END DO
274 DO i = 1, p
275 y(i) = zero
276 END DO
277 RETURN
278 END IF
279*
280* Compute the GQR factorization of matrices A and B:
281*
282* Q**T*A = ( R11 ) M, Q**T*B*Z**T = ( T11 T12 ) M
283* ( 0 ) N-M ( 0 T22 ) N-M
284* M M+P-N N-M
285*
286* where R11 and T22 are upper triangular, and Q and Z are
287* orthogonal.
288*
289 CALL sggqrf( n, m, p, a, lda, work, b, ldb, work( m+1 ),
290 $ work( m+np+1 ), lwork-m-np, info )
291 lopt = work( m+np+1 )
292*
293* Update left-hand-side vector d = Q**T*d = ( d1 ) M
294* ( d2 ) N-M
295*
296 CALL sormqr( 'Left', 'Transpose', n, 1, m, a, lda, work, d,
297 $ max( 1, n ), work( m+np+1 ), lwork-m-np, info )
298 lopt = max( lopt, int( work( m+np+1 ) ) )
299*
300* Solve T22*y2 = d2 for y2
301*
302 IF( n.GT.m ) THEN
303 CALL strtrs( 'Upper', 'No transpose', 'Non unit', n-m, 1,
304 $ b( m+1, m+p-n+1 ), ldb, d( m+1 ), n-m, info )
305*
306 IF( info.GT.0 ) THEN
307 info = 1
308 RETURN
309 END IF
310*
311 CALL scopy( n-m, d( m+1 ), 1, y( m+p-n+1 ), 1 )
312 END IF
313*
314* Set y1 = 0
315*
316 DO 10 i = 1, m + p - n
317 y( i ) = zero
318 10 CONTINUE
319*
320* Update d1 = d1 - T12*y2
321*
322 CALL sgemv( 'No transpose', m, n-m, -one, b( 1, m+p-n+1 ), ldb,
323 $ y( m+p-n+1 ), 1, one, d, 1 )
324*
325* Solve triangular system: R11*x = d1
326*
327 IF( m.GT.0 ) THEN
328 CALL strtrs( 'Upper', 'No Transpose', 'Non unit', m, 1, a, lda,
329 $ d, m, info )
330*
331 IF( info.GT.0 ) THEN
332 info = 2
333 RETURN
334 END IF
335*
336* Copy D to X
337*
338 CALL scopy( m, d, 1, x, 1 )
339 END IF
340*
341* Backward transformation y = Z**T *y
342*
343 CALL sormrq( 'Left', 'Transpose', p, 1, np,
344 $ b( max( 1, n-p+1 ), 1 ), ldb, work( m+1 ), y,
345 $ max( 1, p ), work( m+np+1 ), lwork-m-np, info )
346 work( 1 ) = m + np + max( lopt, int( work( m+np+1 ) ) )
347*
348 RETURN
349*
350* End of SGGGLM
351*
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
strtrs
subroutine strtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
STRTRS
Definition strtrs.f:140
sggqrf
subroutine sggqrf(n, m, p, a, lda, taua, b, ldb, taub, work, lwork, info)
SGGQRF
Definition sggqrf.f:215
sormrq
subroutine sormrq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMRQ
Definition sormrq.f:168
sormqr
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR
Definition sormqr.f:168
sgemv
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:156

◆ ssbev()

subroutine ssbev ( character jobz,
character uplo,
integer n,
integer kd,
real, dimension( ldab, * ) ab,
integer ldab,
real, dimension( * ) w,
real, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer info )

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

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

Purpose:
!>
!> SSBEV computes all the eigenvalues and, optionally, eigenvectors of
!> a real symmetric band matrix A.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[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]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
!>          eigenvectors of the matrix A, with the i-th column of Z
!>          holding the eigenvector associated with W(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>
[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,3*N-2))
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, the algorithm failed to converge; i
!>                off-diagonal elements of an intermediate tridiagonal
!>                form did not converge to zero.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 144 of file ssbev.f.

146*
147* -- LAPACK driver routine --
148* -- LAPACK is a software package provided by Univ. of Tennessee, --
149* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
150*
151* .. Scalar Arguments ..
152 CHARACTER JOBZ, UPLO
153 INTEGER INFO, KD, LDAB, LDZ, N
154* ..
155* .. Array Arguments ..
156 REAL AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
157* ..
158*
159* =====================================================================
160*
161* .. Parameters ..
162 REAL ZERO, ONE
163 parameter( zero = 0.0e0, one = 1.0e0 )
164* ..
165* .. Local Scalars ..
166 LOGICAL LOWER, WANTZ
167 INTEGER IINFO, IMAX, INDE, INDWRK, ISCALE
168 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
169 $ SMLNUM
170* ..
171* .. External Functions ..
172 LOGICAL LSAME
173 REAL SLAMCH, SLANSB
174 EXTERNAL lsame, slamch, slansb
175* ..
176* .. External Subroutines ..
177 EXTERNAL slascl, ssbtrd, sscal, ssteqr, ssterf, xerbla
178* ..
179* .. Intrinsic Functions ..
180 INTRINSIC sqrt
181* ..
182* .. Executable Statements ..
183*
184* Test the input parameters.
185*
186 wantz = lsame( jobz, 'V' )
187 lower = lsame( uplo, 'L' )
188*
189 info = 0
190 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
191 info = -1
192 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
193 info = -2
194 ELSE IF( n.LT.0 ) THEN
195 info = -3
196 ELSE IF( kd.LT.0 ) THEN
197 info = -4
198 ELSE IF( ldab.LT.kd+1 ) THEN
199 info = -6
200 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
201 info = -9
202 END IF
203*
204 IF( info.NE.0 ) THEN
205 CALL xerbla( 'SSBEV ', -info )
206 RETURN
207 END IF
208*
209* Quick return if possible
210*
211 IF( n.EQ.0 )
212 $ RETURN
213*
214 IF( n.EQ.1 ) THEN
215 IF( lower ) THEN
216 w( 1 ) = ab( 1, 1 )
217 ELSE
218 w( 1 ) = ab( kd+1, 1 )
219 END IF
220 IF( wantz )
221 $ z( 1, 1 ) = one
222 RETURN
223 END IF
224*
225* Get machine constants.
226*
227 safmin = slamch( 'Safe minimum' )
228 eps = slamch( 'Precision' )
229 smlnum = safmin / eps
230 bignum = one / smlnum
231 rmin = sqrt( smlnum )
232 rmax = sqrt( bignum )
233*
234* Scale matrix to allowable range, if necessary.
235*
236 anrm = slansb( 'M', uplo, n, kd, ab, ldab, work )
237 iscale = 0
238 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
239 iscale = 1
240 sigma = rmin / anrm
241 ELSE IF( anrm.GT.rmax ) THEN
242 iscale = 1
243 sigma = rmax / anrm
244 END IF
245 IF( iscale.EQ.1 ) THEN
246 IF( lower ) THEN
247 CALL slascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
248 ELSE
249 CALL slascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
250 END IF
251 END IF
252*
253* Call SSBTRD to reduce symmetric band matrix to tridiagonal form.
254*
255 inde = 1
256 indwrk = inde + n
257 CALL ssbtrd( jobz, uplo, n, kd, ab, ldab, w, work( inde ), z, ldz,
258 $ work( indwrk ), iinfo )
259*
260* For eigenvalues only, call SSTERF. For eigenvectors, call SSTEQR.
261*
262 IF( .NOT.wantz ) THEN
263 CALL ssterf( n, w, work( inde ), info )
264 ELSE
265 CALL ssteqr( jobz, n, w, work( inde ), z, ldz, work( indwrk ),
266 $ info )
267 END IF
268*
269* If matrix was scaled, then rescale eigenvalues appropriately.
270*
271 IF( iscale.EQ.1 ) THEN
272 IF( info.EQ.0 ) THEN
273 imax = n
274 ELSE
275 imax = info - 1
276 END IF
277 CALL sscal( imax, one / sigma, w, 1 )
278 END IF
279*
280 RETURN
281*
282* End of SSBEV
283*
slascl
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
ssteqr
subroutine ssteqr(compz, n, d, e, z, ldz, work, info)
SSTEQR
Definition ssteqr.f:131
ssterf
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:86
slansb
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
ssbtrd
subroutine ssbtrd(vect, uplo, n, kd, ab, ldab, d, e, q, ldq, work, info)
SSBTRD
Definition ssbtrd.f:163

◆ ssbev_2stage()

subroutine ssbev_2stage ( character jobz,
character uplo,
integer n,
integer kd,
real, dimension( ldab, * ) ab,
integer ldab,
real, dimension( * ) w,
real, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer lwork,
integer info )

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

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

Purpose:
!>
!> SSBEV_2STAGE computes all the eigenvalues and, optionally, eigenvectors of
!> a real symmetric band matrix A using the 2stage technique for
!> the reduction to tridiagonal.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>                  Not available in this release.
!>
[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]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
!>          eigenvectors of the matrix A, with the i-th column of Z
!>          holding the eigenvector associated with W(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>
[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
!>          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, dimension) where
!>                                   dimension = (2KD+1)*N + KD*NTHREADS + 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]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, the algorithm failed to converge; i
!>                off-diagonal elements of an intermediate tridiagonal
!>                form did not converge to zero.
!>
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 202 of file ssbev_2stage.f.

204*
205 IMPLICIT NONE
206*
207* -- LAPACK driver routine --
208* -- LAPACK is a software package provided by Univ. of Tennessee, --
209* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
210*
211* .. Scalar Arguments ..
212 CHARACTER JOBZ, UPLO
213 INTEGER INFO, KD, LDAB, LDZ, N, LWORK
214* ..
215* .. Array Arguments ..
216 REAL AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
217* ..
218*
219* =====================================================================
220*
221* .. Parameters ..
222 REAL ZERO, ONE
223 parameter( zero = 0.0e0, one = 1.0e0 )
224* ..
225* .. Local Scalars ..
226 LOGICAL LOWER, WANTZ, LQUERY
227 INTEGER IINFO, IMAX, INDE, INDWRK, ISCALE,
228 $ LLWORK, LWMIN, LHTRD, LWTRD, IB, INDHOUS
229 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
230 $ SMLNUM
231* ..
232* .. External Functions ..
233 LOGICAL LSAME
234 INTEGER ILAENV2STAGE
235 REAL SLAMCH, SLANSB
236 EXTERNAL lsame, slamch, slansb, ilaenv2stage
237* ..
238* .. External Subroutines ..
239 EXTERNAL slascl, sscal, ssteqr, ssterf, xerbla,
240 $ ssytrd_sb2st
241* ..
242* .. Intrinsic Functions ..
243 INTRINSIC sqrt
244* ..
245* .. Executable Statements ..
246*
247* Test the input parameters.
248*
249 wantz = lsame( jobz, 'V' )
250 lower = lsame( uplo, 'L' )
251 lquery = ( lwork.EQ.-1 )
252*
253 info = 0
254 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
255 info = -1
256 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
257 info = -2
258 ELSE IF( n.LT.0 ) THEN
259 info = -3
260 ELSE IF( kd.LT.0 ) THEN
261 info = -4
262 ELSE IF( ldab.LT.kd+1 ) THEN
263 info = -6
264 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
265 info = -9
266 END IF
267*
268 IF( info.EQ.0 ) THEN
269 IF( n.LE.1 ) THEN
270 lwmin = 1
271 work( 1 ) = lwmin
272 ELSE
273 ib = ilaenv2stage( 2, 'SSYTRD_SB2ST', jobz,
274 $ n, kd, -1, -1 )
275 lhtrd = ilaenv2stage( 3, 'SSYTRD_SB2ST', jobz,
276 $ n, kd, ib, -1 )
277 lwtrd = ilaenv2stage( 4, 'SSYTRD_SB2ST', jobz,
278 $ n, kd, ib, -1 )
279 lwmin = n + lhtrd + lwtrd
280 work( 1 ) = lwmin
281 ENDIF
282*
283 IF( lwork.LT.lwmin .AND. .NOT.lquery )
284 $ info = -11
285 END IF
286*
287 IF( info.NE.0 ) THEN
288 CALL xerbla( 'SSBEV_2STAGE ', -info )
289 RETURN
290 ELSE IF( lquery ) THEN
291 RETURN
292 END IF
293*
294* Quick return if possible
295*
296 IF( n.EQ.0 )
297 $ RETURN
298*
299 IF( n.EQ.1 ) THEN
300 IF( lower ) THEN
301 w( 1 ) = ab( 1, 1 )
302 ELSE
303 w( 1 ) = ab( kd+1, 1 )
304 END IF
305 IF( wantz )
306 $ z( 1, 1 ) = one
307 RETURN
308 END IF
309*
310* Get machine constants.
311*
312 safmin = slamch( 'Safe minimum' )
313 eps = slamch( 'Precision' )
314 smlnum = safmin / eps
315 bignum = one / smlnum
316 rmin = sqrt( smlnum )
317 rmax = sqrt( bignum )
318*
319* Scale matrix to allowable range, if necessary.
320*
321 anrm = slansb( 'M', uplo, n, kd, ab, ldab, work )
322 iscale = 0
323 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
324 iscale = 1
325 sigma = rmin / anrm
326 ELSE IF( anrm.GT.rmax ) THEN
327 iscale = 1
328 sigma = rmax / anrm
329 END IF
330 IF( iscale.EQ.1 ) THEN
331 IF( lower ) THEN
332 CALL slascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
333 ELSE
334 CALL slascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
335 END IF
336 END IF
337*
338* Call SSYTRD_SB2ST to reduce symmetric band matrix to tridiagonal form.
339*
340 inde = 1
341 indhous = inde + n
342 indwrk = indhous + lhtrd
343 llwork = lwork - indwrk + 1
344*
345 CALL ssytrd_sb2st( "N", jobz, uplo, n, kd, ab, ldab, w,
346 $ work( inde ), work( indhous ), lhtrd,
347 $ work( indwrk ), llwork, iinfo )
348*
349* For eigenvalues only, call SSTERF. For eigenvectors, call SSTEQR.
350*
351 IF( .NOT.wantz ) THEN
352 CALL ssterf( n, w, work( inde ), info )
353 ELSE
354 CALL ssteqr( jobz, n, w, work( inde ), z, ldz, work( indwrk ),
355 $ info )
356 END IF
357*
358* If matrix was scaled, then rescale eigenvalues appropriately.
359*
360 IF( iscale.EQ.1 ) THEN
361 IF( info.EQ.0 ) THEN
362 imax = n
363 ELSE
364 imax = info - 1
365 END IF
366 CALL sscal( imax, one / sigma, w, 1 )
367 END IF
368*
369* Set WORK(1) to optimal workspace size.
370*
371 work( 1 ) = lwmin
372*
373 RETURN
374*
375* End of SSBEV_2STAGE
376*
ilaenv2stage
integer function ilaenv2stage(ispec, name, opts, n1, n2, n3, n4)
ILAENV2STAGE
Definition ilaenv2stage.f:149
ssytrd_sb2st
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
Definition ssytrd_sb2st.F:230

◆ ssbevd()

subroutine ssbevd ( character jobz,
character uplo,
integer n,
integer kd,
real, dimension( ldab, * ) ab,
integer ldab,
real, dimension( * ) w,
real, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

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

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

Purpose:
!>
!> SSBEVD computes all the eigenvalues and, optionally, eigenvectors of
!> a real symmetric band matrix A. If eigenvectors are desired, it uses
!> a divide and conquer algorithm.
!>
!> The divide and conquer algorithm makes very mild assumptions about
!> floating point arithmetic. It will work on machines with a guard
!> digit in add/subtract, or on those binary machines without guard
!> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
!> Cray-2. It could conceivably fail on hexadecimal or decimal machines
!> without guard digits, but we know of none.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[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]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
!>          eigenvectors of the matrix A, with the i-th column of Z
!>          holding the eigenvector associated with W(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>
[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)
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          IF N <= 1,                LWORK must be at least 1.
!>          If JOBZ  = 'N' and N > 2, LWORK must be at least 2*N.
!>          If JOBZ  = 'V' and N > 2, LWORK must be at least
!>                         ( 1 + 5*N + 2*N**2 ).
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal sizes of the WORK and IWORK
!>          arrays, returns these values as the first entries of the WORK
!>          and IWORK arrays, and no error message related to LWORK or
!>          LIWORK is issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.
!>          If JOBZ  = 'N' or N <= 1, LIWORK must be at least 1.
!>          If JOBZ  = 'V' and N > 2, LIWORK must be at least 3 + 5*N.
!>
!>          If LIWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal sizes of the WORK and
!>          IWORK arrays, returns these values as the first entries of
!>          the WORK and IWORK arrays, and no error message related to
!>          LWORK or LIWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, the algorithm failed to converge; i
!>                off-diagonal elements of an intermediate tridiagonal
!>                form did not converge to zero.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 191 of file ssbevd.f.

193*
194* -- LAPACK driver routine --
195* -- LAPACK is a software package provided by Univ. of Tennessee, --
196* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
197*
198* .. Scalar Arguments ..
199 CHARACTER JOBZ, UPLO
200 INTEGER INFO, KD, LDAB, LDZ, LIWORK, LWORK, N
201* ..
202* .. Array Arguments ..
203 INTEGER IWORK( * )
204 REAL AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
205* ..
206*
207* =====================================================================
208*
209* .. Parameters ..
210 REAL ZERO, ONE
211 parameter( zero = 0.0e+0, one = 1.0e+0 )
212* ..
213* .. Local Scalars ..
214 LOGICAL LOWER, LQUERY, WANTZ
215 INTEGER IINFO, INDE, INDWK2, INDWRK, ISCALE, LIWMIN,
216 $ LLWRK2, LWMIN
217 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
218 $ SMLNUM
219* ..
220* .. External Functions ..
221 LOGICAL LSAME
222 REAL SLAMCH, SLANSB
223 EXTERNAL lsame, slamch, slansb
224* ..
225* .. External Subroutines ..
226 EXTERNAL sgemm, slacpy, slascl, ssbtrd, sscal, sstedc,
227 $ ssterf, xerbla
228* ..
229* .. Intrinsic Functions ..
230 INTRINSIC sqrt
231* ..
232* .. Executable Statements ..
233*
234* Test the input parameters.
235*
236 wantz = lsame( jobz, 'V' )
237 lower = lsame( uplo, 'L' )
238 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
239*
240 info = 0
241 IF( n.LE.1 ) THEN
242 liwmin = 1
243 lwmin = 1
244 ELSE
245 IF( wantz ) THEN
246 liwmin = 3 + 5*n
247 lwmin = 1 + 5*n + 2*n**2
248 ELSE
249 liwmin = 1
250 lwmin = 2*n
251 END IF
252 END IF
253 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
254 info = -1
255 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
256 info = -2
257 ELSE IF( n.LT.0 ) THEN
258 info = -3
259 ELSE IF( kd.LT.0 ) THEN
260 info = -4
261 ELSE IF( ldab.LT.kd+1 ) THEN
262 info = -6
263 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
264 info = -9
265 END IF
266*
267 IF( info.EQ.0 ) THEN
268 work( 1 ) = lwmin
269 iwork( 1 ) = liwmin
270*
271 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
272 info = -11
273 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
274 info = -13
275 END IF
276 END IF
277*
278 IF( info.NE.0 ) THEN
279 CALL xerbla( 'SSBEVD', -info )
280 RETURN
281 ELSE IF( lquery ) THEN
282 RETURN
283 END IF
284*
285* Quick return if possible
286*
287 IF( n.EQ.0 )
288 $ RETURN
289*
290 IF( n.EQ.1 ) THEN
291 w( 1 ) = ab( 1, 1 )
292 IF( wantz )
293 $ z( 1, 1 ) = one
294 RETURN
295 END IF
296*
297* Get machine constants.
298*
299 safmin = slamch( 'Safe minimum' )
300 eps = slamch( 'Precision' )
301 smlnum = safmin / eps
302 bignum = one / smlnum
303 rmin = sqrt( smlnum )
304 rmax = sqrt( bignum )
305*
306* Scale matrix to allowable range, if necessary.
307*
308 anrm = slansb( 'M', uplo, n, kd, ab, ldab, work )
309 iscale = 0
310 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
311 iscale = 1
312 sigma = rmin / anrm
313 ELSE IF( anrm.GT.rmax ) THEN
314 iscale = 1
315 sigma = rmax / anrm
316 END IF
317 IF( iscale.EQ.1 ) THEN
318 IF( lower ) THEN
319 CALL slascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
320 ELSE
321 CALL slascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
322 END IF
323 END IF
324*
325* Call SSBTRD to reduce symmetric band matrix to tridiagonal form.
326*
327 inde = 1
328 indwrk = inde + n
329 indwk2 = indwrk + n*n
330 llwrk2 = lwork - indwk2 + 1
331 CALL ssbtrd( jobz, uplo, n, kd, ab, ldab, w, work( inde ), z, ldz,
332 $ work( indwrk ), iinfo )
333*
334* For eigenvalues only, call SSTERF. For eigenvectors, call SSTEDC.
335*
336 IF( .NOT.wantz ) THEN
337 CALL ssterf( n, w, work( inde ), info )
338 ELSE
339 CALL sstedc( 'I', n, w, work( inde ), work( indwrk ), n,
340 $ work( indwk2 ), llwrk2, iwork, liwork, info )
341 CALL sgemm( 'N', 'N', n, n, n, one, z, ldz, work( indwrk ), n,
342 $ zero, work( indwk2 ), n )
343 CALL slacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
344 END IF
345*
346* If matrix was scaled, then rescale eigenvalues appropriately.
347*
348 IF( iscale.EQ.1 )
349 $ CALL sscal( n, one / sigma, w, 1 )
350*
351 work( 1 ) = lwmin
352 iwork( 1 ) = liwmin
353 RETURN
354*
355* End of SSBEVD
356*
slacpy
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
sstedc
subroutine sstedc(compz, n, d, e, z, ldz, work, lwork, iwork, liwork, info)
SSTEDC
Definition sstedc.f:188
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:187

◆ ssbevd_2stage()

subroutine ssbevd_2stage ( character jobz,
character uplo,
integer n,
integer kd,
real, dimension( ldab, * ) ab,
integer ldab,
real, dimension( * ) w,
real, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

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

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

Purpose:
!>
!> SSBEVD_2STAGE computes all the eigenvalues and, optionally, eigenvectors of
!> a real symmetric band matrix A using the 2stage technique for
!> the reduction to tridiagonal. If eigenvectors are desired, it uses
!> a divide and conquer algorithm.
!>
!> The divide and conquer algorithm makes very mild assumptions about
!> floating point arithmetic. It will work on machines with a guard
!> digit in add/subtract, or on those binary machines without guard
!> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
!> Cray-2. It could conceivably fail on hexadecimal or decimal machines
!> without guard digits, but we know of none.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>                  Not available in this release.
!>
[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]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
!>          eigenvectors of the matrix A, with the i-th column of Z
!>          holding the eigenvector associated with W(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>
[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
!>          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, dimension) where
!>                                   dimension = (2KD+1)*N + KD*NTHREADS + 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 sizes of the WORK and IWORK
!>          arrays, returns these values as the first entries of the WORK
!>          and IWORK arrays, and no error message related to LWORK or
!>          LIWORK is issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.
!>          If JOBZ  = 'N' or N <= 1, LIWORK must be at least 1.
!>          If JOBZ  = 'V' and N > 2, LIWORK must be at least 3 + 5*N.
!>
!>          If LIWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal sizes of the WORK and
!>          IWORK arrays, returns these values as the first entries of
!>          the WORK and IWORK arrays, and no error message related to
!>          LWORK or LIWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, the algorithm failed to converge; i
!>                off-diagonal elements of an intermediate tridiagonal
!>                form did not converge to zero.
!>
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 232 of file ssbevd_2stage.f.

234*
235 IMPLICIT NONE
236*
237* -- LAPACK driver routine --
238* -- LAPACK is a software package provided by Univ. of Tennessee, --
239* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
240*
241* .. Scalar Arguments ..
242 CHARACTER JOBZ, UPLO
243 INTEGER INFO, KD, LDAB, LDZ, LIWORK, LWORK, N
244* ..
245* .. Array Arguments ..
246 INTEGER IWORK( * )
247 REAL AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
248* ..
249*
250* =====================================================================
251*
252* .. Parameters ..
253 REAL ZERO, ONE
254 parameter( zero = 0.0e+0, one = 1.0e+0 )
255* ..
256* .. Local Scalars ..
257 LOGICAL LOWER, LQUERY, WANTZ
258 INTEGER IINFO, INDE, INDWK2, INDWRK, ISCALE, LIWMIN,
259 $ LLWORK, LWMIN, LHTRD, LWTRD, IB, INDHOUS,
260 $ LLWRK2
261 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
262 $ SMLNUM
263* ..
264* .. External Functions ..
265 LOGICAL LSAME
266 INTEGER ILAENV2STAGE
267 REAL SLAMCH, SLANSB
268 EXTERNAL lsame, slamch, slansb, ilaenv2stage
269* ..
270* .. External Subroutines ..
271 EXTERNAL sgemm, slacpy, slascl, sscal, sstedc,
272 $ ssterf, xerbla, ssytrd_sb2st
273* ..
274* .. Intrinsic Functions ..
275 INTRINSIC sqrt
276* ..
277* .. Executable Statements ..
278*
279* Test the input parameters.
280*
281 wantz = lsame( jobz, 'V' )
282 lower = lsame( uplo, 'L' )
283 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
284*
285 info = 0
286 IF( n.LE.1 ) THEN
287 liwmin = 1
288 lwmin = 1
289 ELSE
290 ib = ilaenv2stage( 2, 'SSYTRD_SB2ST', jobz, n, kd, -1, -1 )
291 lhtrd = ilaenv2stage( 3, 'SSYTRD_SB2ST', jobz, n, kd, ib, -1 )
292 lwtrd = ilaenv2stage( 4, 'SSYTRD_SB2ST', jobz, n, kd, ib, -1 )
293 IF( wantz ) THEN
294 liwmin = 3 + 5*n
295 lwmin = 1 + 5*n + 2*n**2
296 ELSE
297 liwmin = 1
298 lwmin = max( 2*n, n+lhtrd+lwtrd )
299 END IF
300 END IF
301 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
302 info = -1
303 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
304 info = -2
305 ELSE IF( n.LT.0 ) THEN
306 info = -3
307 ELSE IF( kd.LT.0 ) THEN
308 info = -4
309 ELSE IF( ldab.LT.kd+1 ) THEN
310 info = -6
311 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
312 info = -9
313 END IF
314*
315 IF( info.EQ.0 ) THEN
316 work( 1 ) = lwmin
317 iwork( 1 ) = liwmin
318*
319 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
320 info = -11
321 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
322 info = -13
323 END IF
324 END IF
325*
326 IF( info.NE.0 ) THEN
327 CALL xerbla( 'SSBEVD_2STAGE', -info )
328 RETURN
329 ELSE IF( lquery ) THEN
330 RETURN
331 END IF
332*
333* Quick return if possible
334*
335 IF( n.EQ.0 )
336 $ RETURN
337*
338 IF( n.EQ.1 ) THEN
339 w( 1 ) = ab( 1, 1 )
340 IF( wantz )
341 $ z( 1, 1 ) = one
342 RETURN
343 END IF
344*
345* Get machine constants.
346*
347 safmin = slamch( 'Safe minimum' )
348 eps = slamch( 'Precision' )
349 smlnum = safmin / eps
350 bignum = one / smlnum
351 rmin = sqrt( smlnum )
352 rmax = sqrt( bignum )
353*
354* Scale matrix to allowable range, if necessary.
355*
356 anrm = slansb( 'M', uplo, n, kd, ab, ldab, work )
357 iscale = 0
358 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
359 iscale = 1
360 sigma = rmin / anrm
361 ELSE IF( anrm.GT.rmax ) THEN
362 iscale = 1
363 sigma = rmax / anrm
364 END IF
365 IF( iscale.EQ.1 ) THEN
366 IF( lower ) THEN
367 CALL slascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
368 ELSE
369 CALL slascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
370 END IF
371 END IF
372*
373* Call SSYTRD_SB2ST to reduce band symmetric matrix to tridiagonal form.
374*
375 inde = 1
376 indhous = inde + n
377 indwrk = indhous + lhtrd
378 llwork = lwork - indwrk + 1
379 indwk2 = indwrk + n*n
380 llwrk2 = lwork - indwk2 + 1
381*
382 CALL ssytrd_sb2st( "N", jobz, uplo, n, kd, ab, ldab, w,
383 $ work( inde ), work( indhous ), lhtrd,
384 $ work( indwrk ), llwork, iinfo )
385*
386* For eigenvalues only, call SSTERF. For eigenvectors, call SSTEDC.
387*
388 IF( .NOT.wantz ) THEN
389 CALL ssterf( n, w, work( inde ), info )
390 ELSE
391 CALL sstedc( 'I', n, w, work( inde ), work( indwrk ), n,
392 $ work( indwk2 ), llwrk2, iwork, liwork, info )
393 CALL sgemm( 'N', 'N', n, n, n, one, z, ldz, work( indwrk ), n,
394 $ zero, work( indwk2 ), n )
395 CALL slacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
396 END IF
397*
398* If matrix was scaled, then rescale eigenvalues appropriately.
399*
400 IF( iscale.EQ.1 )
401 $ CALL sscal( n, one / sigma, w, 1 )
402*
403 work( 1 ) = lwmin
404 iwork( 1 ) = liwmin
405 RETURN
406*
407* End of SSBEVD_2STAGE
408*

◆ ssbevx()

subroutine ssbevx ( 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, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

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

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

Purpose:
!>
!> SSBEVX computes selected eigenvalues and, optionally, eigenvectors
!> of a real symmetric band 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]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  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 (7*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 262 of file ssbevx.f.

265*
266* -- LAPACK driver routine --
267* -- LAPACK is a software package provided by Univ. of Tennessee, --
268* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
269*
270* .. Scalar Arguments ..
271 CHARACTER JOBZ, RANGE, UPLO
272 INTEGER IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N
273 REAL ABSTOL, VL, VU
274* ..
275* .. Array Arguments ..
276 INTEGER IFAIL( * ), IWORK( * )
277 REAL AB( LDAB, * ), Q( LDQ, * ), W( * ), WORK( * ),
278 $ Z( LDZ, * )
279* ..
280*
281* =====================================================================
282*
283* .. Parameters ..
284 REAL ZERO, ONE
285 parameter( zero = 0.0e0, one = 1.0e0 )
286* ..
287* .. Local Scalars ..
288 LOGICAL ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ
289 CHARACTER ORDER
290 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
291 $ INDISP, INDIWO, INDWRK, ISCALE, ITMP1, J, JJ,
292 $ NSPLIT
293 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
294 $ SIGMA, SMLNUM, TMP1, VLL, VUU
295* ..
296* .. External Functions ..
297 LOGICAL LSAME
298 REAL SLAMCH, SLANSB
299 EXTERNAL lsame, slamch, slansb
300* ..
301* .. External Subroutines ..
302 EXTERNAL scopy, sgemv, slacpy, slascl, ssbtrd, sscal,
303 $ sstebz, sstein, ssteqr, ssterf, sswap, xerbla
304* ..
305* .. Intrinsic Functions ..
306 INTRINSIC max, min, sqrt
307* ..
308* .. Executable Statements ..
309*
310* Test the input parameters.
311*
312 wantz = lsame( jobz, 'V' )
313 alleig = lsame( range, 'A' )
314 valeig = lsame( range, 'V' )
315 indeig = lsame( range, 'I' )
316 lower = lsame( uplo, 'L' )
317*
318 info = 0
319 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
320 info = -1
321 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
322 info = -2
323 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
324 info = -3
325 ELSE IF( n.LT.0 ) THEN
326 info = -4
327 ELSE IF( kd.LT.0 ) THEN
328 info = -5
329 ELSE IF( ldab.LT.kd+1 ) THEN
330 info = -7
331 ELSE IF( wantz .AND. ldq.LT.max( 1, n ) ) THEN
332 info = -9
333 ELSE
334 IF( valeig ) THEN
335 IF( n.GT.0 .AND. vu.LE.vl )
336 $ info = -11
337 ELSE IF( indeig ) THEN
338 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
339 info = -12
340 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
341 info = -13
342 END IF
343 END IF
344 END IF
345 IF( info.EQ.0 ) THEN
346 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
347 $ info = -18
348 END IF
349*
350 IF( info.NE.0 ) THEN
351 CALL xerbla( 'SSBEVX', -info )
352 RETURN
353 END IF
354*
355* Quick return if possible
356*
357 m = 0
358 IF( n.EQ.0 )
359 $ RETURN
360*
361 IF( n.EQ.1 ) THEN
362 m = 1
363 IF( lower ) THEN
364 tmp1 = ab( 1, 1 )
365 ELSE
366 tmp1 = ab( kd+1, 1 )
367 END IF
368 IF( valeig ) THEN
369 IF( .NOT.( vl.LT.tmp1 .AND. vu.GE.tmp1 ) )
370 $ m = 0
371 END IF
372 IF( m.EQ.1 ) THEN
373 w( 1 ) = tmp1
374 IF( wantz )
375 $ z( 1, 1 ) = one
376 END IF
377 RETURN
378 END IF
379*
380* Get machine constants.
381*
382 safmin = slamch( 'Safe minimum' )
383 eps = slamch( 'Precision' )
384 smlnum = safmin / eps
385 bignum = one / smlnum
386 rmin = sqrt( smlnum )
387 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
388*
389* Scale matrix to allowable range, if necessary.
390*
391 iscale = 0
392 abstll = abstol
393 IF ( valeig ) THEN
394 vll = vl
395 vuu = vu
396 ELSE
397 vll = zero
398 vuu = zero
399 ENDIF
400 anrm = slansb( 'M', uplo, n, kd, ab, ldab, work )
401 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
402 iscale = 1
403 sigma = rmin / anrm
404 ELSE IF( anrm.GT.rmax ) THEN
405 iscale = 1
406 sigma = rmax / anrm
407 END IF
408 IF( iscale.EQ.1 ) THEN
409 IF( lower ) THEN
410 CALL slascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
411 ELSE
412 CALL slascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
413 END IF
414 IF( abstol.GT.0 )
415 $ abstll = abstol*sigma
416 IF( valeig ) THEN
417 vll = vl*sigma
418 vuu = vu*sigma
419 END IF
420 END IF
421*
422* Call SSBTRD to reduce symmetric band matrix to tridiagonal form.
423*
424 indd = 1
425 inde = indd + n
426 indwrk = inde + n
427 CALL ssbtrd( jobz, uplo, n, kd, ab, ldab, work( indd ),
428 $ work( inde ), q, ldq, work( indwrk ), iinfo )
429*
430* If all eigenvalues are desired and ABSTOL is less than or equal
431* to zero, then call SSTERF or SSTEQR. If this fails for some
432* eigenvalue, then try SSTEBZ.
433*
434 test = .false.
435 IF (indeig) THEN
436 IF (il.EQ.1 .AND. iu.EQ.n) THEN
437 test = .true.
438 END IF
439 END IF
440 IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
441 CALL scopy( n, work( indd ), 1, w, 1 )
442 indee = indwrk + 2*n
443 IF( .NOT.wantz ) THEN
444 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
445 CALL ssterf( n, w, work( indee ), info )
446 ELSE
447 CALL slacpy( 'A', n, n, q, ldq, z, ldz )
448 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
449 CALL ssteqr( jobz, n, w, work( indee ), z, ldz,
450 $ work( indwrk ), info )
451 IF( info.EQ.0 ) THEN
452 DO 10 i = 1, n
453 ifail( i ) = 0
454 10 CONTINUE
455 END IF
456 END IF
457 IF( info.EQ.0 ) THEN
458 m = n
459 GO TO 30
460 END IF
461 info = 0
462 END IF
463*
464* Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
465*
466 IF( wantz ) THEN
467 order = 'B'
468 ELSE
469 order = 'E'
470 END IF
471 indibl = 1
472 indisp = indibl + n
473 indiwo = indisp + n
474 CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
475 $ work( indd ), work( inde ), m, nsplit, w,
476 $ iwork( indibl ), iwork( indisp ), work( indwrk ),
477 $ iwork( indiwo ), info )
478*
479 IF( wantz ) THEN
480 CALL sstein( n, work( indd ), work( inde ), m, w,
481 $ iwork( indibl ), iwork( indisp ), z, ldz,
482 $ work( indwrk ), iwork( indiwo ), ifail, info )
483*
484* Apply orthogonal matrix used in reduction to tridiagonal
485* form to eigenvectors returned by SSTEIN.
486*
487 DO 20 j = 1, m
488 CALL scopy( n, z( 1, j ), 1, work( 1 ), 1 )
489 CALL sgemv( 'N', n, n, one, q, ldq, work, 1, zero,
490 $ z( 1, j ), 1 )
491 20 CONTINUE
492 END IF
493*
494* If matrix was scaled, then rescale eigenvalues appropriately.
495*
496 30 CONTINUE
497 IF( iscale.EQ.1 ) THEN
498 IF( info.EQ.0 ) THEN
499 imax = m
500 ELSE
501 imax = info - 1
502 END IF
503 CALL sscal( imax, one / sigma, w, 1 )
504 END IF
505*
506* If eigenvalues are not in order, then sort them, along with
507* eigenvectors.
508*
509 IF( wantz ) THEN
510 DO 50 j = 1, m - 1
511 i = 0
512 tmp1 = w( j )
513 DO 40 jj = j + 1, m
514 IF( w( jj ).LT.tmp1 ) THEN
515 i = jj
516 tmp1 = w( jj )
517 END IF
518 40 CONTINUE
519*
520 IF( i.NE.0 ) THEN
521 itmp1 = iwork( indibl+i-1 )
522 w( i ) = w( j )
523 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
524 w( j ) = tmp1
525 iwork( indibl+j-1 ) = itmp1
526 CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
527 IF( info.NE.0 ) THEN
528 itmp1 = ifail( i )
529 ifail( i ) = ifail( j )
530 ifail( j ) = itmp1
531 END IF
532 END IF
533 50 CONTINUE
534 END IF
535*
536 RETURN
537*
538* End of SSBEVX
539*
sstebz
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
sstein
subroutine sstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
SSTEIN
Definition sstein.f:174

◆ 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  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
361 EXTERNAL lsame, slamch, slansb, ilaenv2stage
362* ..
363* .. External Subroutines ..
364 EXTERNAL scopy, sgemv, slacpy, slascl, sscal,
365 $ sstebz, sstein, ssteqr, ssterf, sswap, xerbla,
366 $ ssytrd_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, 'SSYTRD_SB2ST', jobz,
420 $ n, kd, -1, -1 )
421 lhtrd = ilaenv2stage( 3, 'SSYTRD_SB2ST', jobz,
422 $ n, kd, ib, -1 )
423 lwtrd = ilaenv2stage( 4, 'SSYTRD_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( 'SSBEVX_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 = slamch( 'Safe minimum' )
468 eps = slamch( '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 = slansb( '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 slascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
496 ELSE
497 CALL slascl( '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 SSYTRD_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 ssytrd_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 SSTERF or SSTEQR. If this fails for some
521* eigenvalue, then try SSTEBZ.
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 scopy( n, work( indd ), 1, w, 1 )
531 indee = indwrk + 2*n
532 IF( .NOT.wantz ) THEN
533 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
534 CALL ssterf( n, w, work( indee ), info )
535 ELSE
536 CALL slacpy( 'A', n, n, q, ldq, z, ldz )
537 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
538 CALL ssteqr( 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 SSTEBZ 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 sstebz( 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 sstein( 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 SSTEIN.
575*
576 DO 20 j = 1, m
577 CALL scopy( n, z( 1, j ), 1, work( 1 ), 1 )
578 CALL sgemv( '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 sscal( 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 sswap( 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 SSBEVX_2STAGE
632*

◆ ssbgv()

subroutine ssbgv ( character jobz,
character uplo,
integer n,
integer ka,
integer kb,
real, dimension( ldab, * ) ab,
integer ldab,
real, dimension( ldbb, * ) bb,
integer ldbb,
real, dimension( * ) w,
real, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer info )

SSBGV

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

Purpose:
!>
!> SSBGV computes all the eigenvalues, and optionally, the eigenvectors
!> of a real generalized symmetric-definite banded eigenproblem, of
!> the form A*x=(lambda)*B*x. Here A and B are assumed to be symmetric
!> and banded, and B is also positive definite.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangles of A and B are stored;
!>          = 'L':  Lower triangles of A and B are stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A and B.  N >= 0.
!>
[in]KA
!>          KA is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'. KA >= 0.
!>
[in]KB
!>          KB is INTEGER
!>          The number of superdiagonals of the matrix B if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'. KB >= 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 ka+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(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).
!>
!>          On exit, the contents of AB are destroyed.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KA+1.
!>
[in,out]BB
!>          BB is REAL array, dimension (LDBB, N)
!>          On entry, the upper or lower triangle of the symmetric band
!>          matrix B, stored in the first kb+1 rows of the array.  The
!>          j-th column of B is stored in the j-th column of the array BB
!>          as follows:
!>          if UPLO = 'U', BB(kb+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j;
!>          if UPLO = 'L', BB(1+i-j,j)    = B(i,j) for j<=i<=min(n,j+kb).
!>
!>          On exit, the factor S from the split Cholesky factorization
!>          B = S**T*S, as returned by SPBSTF.
!>
[in]LDBB
!>          LDBB is INTEGER
!>          The leading dimension of the array BB.  LDBB >= KB+1.
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
!>          eigenvectors, with the i-th column of Z holding the
!>          eigenvector associated with W(i). The eigenvectors are
!>          normalized so that Z**T*B*Z = I.
!>          If JOBZ = 'N', then Z is not referenced.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= N.
!>
[out]WORK
!>          WORK is REAL array, dimension (3*N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, and i is:
!>             <= N:  the algorithm failed to converge:
!>                    i off-diagonal elements of an intermediate
!>                    tridiagonal form did not converge to zero;
!>             > N:   if INFO = N + i, for 1 <= i <= N, then SPBSTF
!>                    returned INFO = i: B is not positive definite.
!>                    The factorization of B could not be completed and
!>                    no eigenvalues or eigenvectors were computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 175 of file ssbgv.f.

177*
178* -- LAPACK driver routine --
179* -- LAPACK is a software package provided by Univ. of Tennessee, --
180* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
181*
182* .. Scalar Arguments ..
183 CHARACTER JOBZ, UPLO
184 INTEGER INFO, KA, KB, LDAB, LDBB, LDZ, N
185* ..
186* .. Array Arguments ..
187 REAL AB( LDAB, * ), BB( LDBB, * ), W( * ),
188 $ WORK( * ), Z( LDZ, * )
189* ..
190*
191* =====================================================================
192*
193* .. Local Scalars ..
194 LOGICAL UPPER, WANTZ
195 CHARACTER VECT
196 INTEGER IINFO, INDE, INDWRK
197* ..
198* .. External Functions ..
199 LOGICAL LSAME
200 EXTERNAL lsame
201* ..
202* .. External Subroutines ..
203 EXTERNAL spbstf, ssbgst, ssbtrd, ssteqr, ssterf, xerbla
204* ..
205* .. Executable Statements ..
206*
207* Test the input parameters.
208*
209 wantz = lsame( jobz, 'V' )
210 upper = lsame( uplo, 'U' )
211*
212 info = 0
213 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
214 info = -1
215 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
216 info = -2
217 ELSE IF( n.LT.0 ) THEN
218 info = -3
219 ELSE IF( ka.LT.0 ) THEN
220 info = -4
221 ELSE IF( kb.LT.0 .OR. kb.GT.ka ) THEN
222 info = -5
223 ELSE IF( ldab.LT.ka+1 ) THEN
224 info = -7
225 ELSE IF( ldbb.LT.kb+1 ) THEN
226 info = -9
227 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
228 info = -12
229 END IF
230 IF( info.NE.0 ) THEN
231 CALL xerbla( 'SSBGV ', -info )
232 RETURN
233 END IF
234*
235* Quick return if possible
236*
237 IF( n.EQ.0 )
238 $ RETURN
239*
240* Form a split Cholesky factorization of B.
241*
242 CALL spbstf( uplo, n, kb, bb, ldbb, info )
243 IF( info.NE.0 ) THEN
244 info = n + info
245 RETURN
246 END IF
247*
248* Transform problem to standard eigenvalue problem.
249*
250 inde = 1
251 indwrk = inde + n
252 CALL ssbgst( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, z, ldz,
253 $ work( indwrk ), iinfo )
254*
255* Reduce to tridiagonal form.
256*
257 IF( wantz ) THEN
258 vect = 'U'
259 ELSE
260 vect = 'N'
261 END IF
262 CALL ssbtrd( vect, uplo, n, ka, ab, ldab, w, work( inde ), z, ldz,
263 $ work( indwrk ), iinfo )
264*
265* For eigenvalues only, call SSTERF. For eigenvectors, call SSTEQR.
266*
267 IF( .NOT.wantz ) THEN
268 CALL ssterf( n, w, work( inde ), info )
269 ELSE
270 CALL ssteqr( jobz, n, w, work( inde ), z, ldz, work( indwrk ),
271 $ info )
272 END IF
273 RETURN
274*
275* End of SSBGV
276*
spbstf
subroutine spbstf(uplo, n, kd, ab, ldab, info)
SPBSTF
Definition spbstf.f:152
ssbgst
subroutine ssbgst(vect, uplo, n, ka, kb, ab, ldab, bb, ldbb, x, ldx, work, info)
SSBGST
Definition ssbgst.f:159

◆ ssbgvd()

subroutine ssbgvd ( character jobz,
character uplo,
integer n,
integer ka,
integer kb,
real, dimension( ldab, * ) ab,
integer ldab,
real, dimension( ldbb, * ) bb,
integer ldbb,
real, dimension( * ) w,
real, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

SSBGVD

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

Purpose:
!>
!> SSBGVD computes all the eigenvalues, and optionally, the eigenvectors
!> of a real generalized symmetric-definite banded eigenproblem, of the
!> form A*x=(lambda)*B*x.  Here A and B are assumed to be symmetric and
!> banded, and B is also positive definite.  If eigenvectors are
!> desired, it uses a divide and conquer algorithm.
!>
!> The divide and conquer algorithm makes very mild assumptions about
!> floating point arithmetic. It will work on machines with a guard
!> digit in add/subtract, or on those binary machines without guard
!> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
!> Cray-2. It could conceivably fail on hexadecimal or decimal machines
!> without guard digits, but we know of none.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangles of A and B are stored;
!>          = 'L':  Lower triangles of A and B are stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A and B.  N >= 0.
!>
[in]KA
!>          KA is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'.  KA >= 0.
!>
[in]KB
!>          KB is INTEGER
!>          The number of superdiagonals of the matrix B if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'.  KB >= 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 ka+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(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).
!>
!>          On exit, the contents of AB are destroyed.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KA+1.
!>
[in,out]BB
!>          BB is REAL array, dimension (LDBB, N)
!>          On entry, the upper or lower triangle of the symmetric band
!>          matrix B, stored in the first kb+1 rows of the array.  The
!>          j-th column of B is stored in the j-th column of the array BB
!>          as follows:
!>          if UPLO = 'U', BB(ka+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j;
!>          if UPLO = 'L', BB(1+i-j,j)    = B(i,j) for j<=i<=min(n,j+kb).
!>
!>          On exit, the factor S from the split Cholesky factorization
!>          B = S**T*S, as returned by SPBSTF.
!>
[in]LDBB
!>          LDBB is INTEGER
!>          The leading dimension of the array BB.  LDBB >= KB+1.
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
!>          eigenvectors, with the i-th column of Z holding the
!>          eigenvector associated with W(i).  The eigenvectors are
!>          normalized so Z**T*B*Z = I.
!>          If JOBZ = 'N', then Z is not referenced.
!>
[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 dimension of the array WORK.
!>          If N <= 1,               LWORK >= 1.
!>          If JOBZ = 'N' and N > 1, LWORK >= 3*N.
!>          If JOBZ = 'V' and N > 1, LWORK >= 1 + 5*N + 2*N**2.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal sizes of the WORK and IWORK
!>          arrays, returns these values as the first entries of the WORK
!>          and IWORK arrays, and no error message related to LWORK or
!>          LIWORK is issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          On exit, if LIWORK > 0, IWORK(1) returns the optimal LIWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.
!>          If JOBZ  = 'N' or N <= 1, LIWORK >= 1.
!>          If JOBZ  = 'V' and N > 1, LIWORK >= 3 + 5*N.
!>
!>          If LIWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal sizes of the WORK and
!>          IWORK arrays, returns these values as the first entries of
!>          the WORK and IWORK arrays, and no error message related to
!>          LWORK or LIWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, and i is:
!>             <= N:  the algorithm failed to converge:
!>                    i off-diagonal elements of an intermediate
!>                    tridiagonal form did not converge to zero;
!>             > N:   if INFO = N + i, for 1 <= i <= N, then SPBSTF
!>                    returned INFO = i: B is not positive definite.
!>                    The factorization of B could not be completed and
!>                    no eigenvalues or eigenvectors were computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA

Definition at line 225 of file ssbgvd.f.

227*
228* -- LAPACK driver routine --
229* -- LAPACK is a software package provided by Univ. of Tennessee, --
230* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
231*
232* .. Scalar Arguments ..
233 CHARACTER JOBZ, UPLO
234 INTEGER INFO, KA, KB, LDAB, LDBB, LDZ, LIWORK, LWORK, N
235* ..
236* .. Array Arguments ..
237 INTEGER IWORK( * )
238 REAL AB( LDAB, * ), BB( LDBB, * ), W( * ),
239 $ WORK( * ), Z( LDZ, * )
240* ..
241*
242* =====================================================================
243*
244* .. Parameters ..
245 REAL ONE, ZERO
246 parameter( one = 1.0e+0, zero = 0.0e+0 )
247* ..
248* .. Local Scalars ..
249 LOGICAL LQUERY, UPPER, WANTZ
250 CHARACTER VECT
251 INTEGER IINFO, INDE, INDWK2, INDWRK, LIWMIN, LLWRK2,
252 $ LWMIN
253* ..
254* .. External Functions ..
255 LOGICAL LSAME
256 EXTERNAL lsame
257* ..
258* .. External Subroutines ..
259 EXTERNAL sgemm, slacpy, spbstf, ssbgst, ssbtrd, sstedc,
260 $ ssterf, xerbla
261* ..
262* .. Executable Statements ..
263*
264* Test the input parameters.
265*
266 wantz = lsame( jobz, 'V' )
267 upper = lsame( uplo, 'U' )
268 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
269*
270 info = 0
271 IF( n.LE.1 ) THEN
272 liwmin = 1
273 lwmin = 1
274 ELSE IF( wantz ) THEN
275 liwmin = 3 + 5*n
276 lwmin = 1 + 5*n + 2*n**2
277 ELSE
278 liwmin = 1
279 lwmin = 2*n
280 END IF
281*
282 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
283 info = -1
284 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
285 info = -2
286 ELSE IF( n.LT.0 ) THEN
287 info = -3
288 ELSE IF( ka.LT.0 ) THEN
289 info = -4
290 ELSE IF( kb.LT.0 .OR. kb.GT.ka ) THEN
291 info = -5
292 ELSE IF( ldab.LT.ka+1 ) THEN
293 info = -7
294 ELSE IF( ldbb.LT.kb+1 ) THEN
295 info = -9
296 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
297 info = -12
298 END IF
299*
300 IF( info.EQ.0 ) THEN
301 work( 1 ) = lwmin
302 iwork( 1 ) = liwmin
303*
304 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
305 info = -14
306 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
307 info = -16
308 END IF
309 END IF
310*
311 IF( info.NE.0 ) THEN
312 CALL xerbla( 'SSBGVD', -info )
313 RETURN
314 ELSE IF( lquery ) THEN
315 RETURN
316 END IF
317*
318* Quick return if possible
319*
320 IF( n.EQ.0 )
321 $ RETURN
322*
323* Form a split Cholesky factorization of B.
324*
325 CALL spbstf( uplo, n, kb, bb, ldbb, info )
326 IF( info.NE.0 ) THEN
327 info = n + info
328 RETURN
329 END IF
330*
331* Transform problem to standard eigenvalue problem.
332*
333 inde = 1
334 indwrk = inde + n
335 indwk2 = indwrk + n*n
336 llwrk2 = lwork - indwk2 + 1
337 CALL ssbgst( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, z, ldz,
338 $ work, iinfo )
339*
340* Reduce to tridiagonal form.
341*
342 IF( wantz ) THEN
343 vect = 'U'
344 ELSE
345 vect = 'N'
346 END IF
347 CALL ssbtrd( vect, uplo, n, ka, ab, ldab, w, work( inde ), z, ldz,
348 $ work( indwrk ), iinfo )
349*
350* For eigenvalues only, call SSTERF. For eigenvectors, call SSTEDC.
351*
352 IF( .NOT.wantz ) THEN
353 CALL ssterf( n, w, work( inde ), info )
354 ELSE
355 CALL sstedc( 'I', n, w, work( inde ), work( indwrk ), n,
356 $ work( indwk2 ), llwrk2, iwork, liwork, info )
357 CALL sgemm( 'N', 'N', n, n, n, one, z, ldz, work( indwrk ), n,
358 $ zero, work( indwk2 ), n )
359 CALL slacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
360 END IF
361*
362 work( 1 ) = lwmin
363 iwork( 1 ) = liwmin
364*
365 RETURN
366*
367* End of SSBGVD
368*

◆ ssbgvx()

subroutine ssbgvx ( character jobz,
character range,
character uplo,
integer n,
integer ka,
integer kb,
real, dimension( ldab, * ) ab,
integer ldab,
real, dimension( ldbb, * ) bb,
integer ldbb,
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, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

SSBGVX

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

Purpose:
!>
!> SSBGVX computes selected eigenvalues, and optionally, eigenvectors
!> of a real generalized symmetric-definite banded eigenproblem, of
!> the form A*x=(lambda)*B*x.  Here A and B are assumed to be symmetric
!> and banded, and B is also positive definite.  Eigenvalues and
!> eigenvectors can be selected by specifying either all eigenvalues,
!> 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]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangles of A and B are stored;
!>          = 'L':  Lower triangles of A and B are stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A and B.  N >= 0.
!>
[in]KA
!>          KA is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'.  KA >= 0.
!>
[in]KB
!>          KB is INTEGER
!>          The number of superdiagonals of the matrix B if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'.  KB >= 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 ka+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(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).
!>
!>          On exit, the contents of AB are destroyed.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KA+1.
!>
[in,out]BB
!>          BB is REAL array, dimension (LDBB, N)
!>          On entry, the upper or lower triangle of the symmetric band
!>          matrix B, stored in the first kb+1 rows of the array.  The
!>          j-th column of B is stored in the j-th column of the array BB
!>          as follows:
!>          if UPLO = 'U', BB(ka+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j;
!>          if UPLO = 'L', BB(1+i-j,j)    = B(i,j) for j<=i<=min(n,j+kb).
!>
!>          On exit, the factor S from the split Cholesky factorization
!>          B = S**T*S, as returned by SPBSTF.
!>
[in]LDBB
!>          LDBB is INTEGER
!>          The leading dimension of the array BB.  LDBB >= KB+1.
!>
[out]Q
!>          Q is REAL array, dimension (LDQ, N)
!>          If JOBZ = 'V', the n-by-n matrix used in the reduction of
!>          A*x = (lambda)*B*x to standard form, i.e. C*x = (lambda)*x,
!>          and consequently C 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 = 'N',
!>          LDQ >= 1. If JOBZ = 'V', 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 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').
!>
[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)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
!>          eigenvectors, with the i-th column of Z holding the
!>          eigenvector associated with W(i).  The eigenvectors are
!>          normalized so Z**T*B*Z = I.
!>          If JOBZ = 'N', then Z is not referenced.
!>
[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 (7*N)
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (5*N)
!>
[out]IFAIL
!>          IFAIL is INTEGER array, dimension (M)
!>          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 eigenvalues 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
!>          <= N: if INFO = i, then i eigenvectors failed to converge.
!>                  Their indices are stored in IFAIL.
!>          > N:  SPBSTF returned an error code; i.e.,
!>                if INFO = N + i, for 1 <= i <= N, then the leading
!>                minor of order i of B is not positive definite.
!>                The factorization of B could not be completed and
!>                no eigenvalues or eigenvectors were computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA

Definition at line 291 of file ssbgvx.f.

294*
295* -- LAPACK driver routine --
296* -- LAPACK is a software package provided by Univ. of Tennessee, --
297* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
298*
299* .. Scalar Arguments ..
300 CHARACTER JOBZ, RANGE, UPLO
301 INTEGER IL, INFO, IU, KA, KB, LDAB, LDBB, LDQ, LDZ, M,
302 $ N
303 REAL ABSTOL, VL, VU
304* ..
305* .. Array Arguments ..
306 INTEGER IFAIL( * ), IWORK( * )
307 REAL AB( LDAB, * ), BB( LDBB, * ), Q( LDQ, * ),
308 $ W( * ), WORK( * ), Z( LDZ, * )
309* ..
310*
311* =====================================================================
312*
313* .. Parameters ..
314 REAL ZERO, ONE
315 parameter( zero = 0.0e+0, one = 1.0e+0 )
316* ..
317* .. Local Scalars ..
318 LOGICAL ALLEIG, INDEIG, TEST, UPPER, VALEIG, WANTZ
319 CHARACTER ORDER, VECT
320 INTEGER I, IINFO, INDD, INDE, INDEE, INDIBL, INDISP,
321 $ INDIWO, INDWRK, ITMP1, J, JJ, NSPLIT
322 REAL TMP1
323* ..
324* .. External Functions ..
325 LOGICAL LSAME
326 EXTERNAL lsame
327* ..
328* .. External Subroutines ..
329 EXTERNAL scopy, sgemv, slacpy, spbstf, ssbgst, ssbtrd,
330 $ sstebz, sstein, ssteqr, ssterf, sswap, xerbla
331* ..
332* .. Intrinsic Functions ..
333 INTRINSIC min
334* ..
335* .. Executable Statements ..
336*
337* Test the input parameters.
338*
339 wantz = lsame( jobz, 'V' )
340 upper = lsame( uplo, 'U' )
341 alleig = lsame( range, 'A' )
342 valeig = lsame( range, 'V' )
343 indeig = lsame( range, 'I' )
344*
345 info = 0
346 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
347 info = -1
348 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
349 info = -2
350 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
351 info = -3
352 ELSE IF( n.LT.0 ) THEN
353 info = -4
354 ELSE IF( ka.LT.0 ) THEN
355 info = -5
356 ELSE IF( kb.LT.0 .OR. kb.GT.ka ) THEN
357 info = -6
358 ELSE IF( ldab.LT.ka+1 ) THEN
359 info = -8
360 ELSE IF( ldbb.LT.kb+1 ) THEN
361 info = -10
362 ELSE IF( ldq.LT.1 .OR. ( wantz .AND. ldq.LT.n ) ) THEN
363 info = -12
364 ELSE
365 IF( valeig ) THEN
366 IF( n.GT.0 .AND. vu.LE.vl )
367 $ info = -14
368 ELSE IF( indeig ) THEN
369 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
370 info = -15
371 ELSE IF ( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
372 info = -16
373 END IF
374 END IF
375 END IF
376 IF( info.EQ.0) THEN
377 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
378 info = -21
379 END IF
380 END IF
381*
382 IF( info.NE.0 ) THEN
383 CALL xerbla( 'SSBGVX', -info )
384 RETURN
385 END IF
386*
387* Quick return if possible
388*
389 m = 0
390 IF( n.EQ.0 )
391 $ RETURN
392*
393* Form a split Cholesky factorization of B.
394*
395 CALL spbstf( uplo, n, kb, bb, ldbb, info )
396 IF( info.NE.0 ) THEN
397 info = n + info
398 RETURN
399 END IF
400*
401* Transform problem to standard eigenvalue problem.
402*
403 CALL ssbgst( jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, q, ldq,
404 $ work, iinfo )
405*
406* Reduce symmetric band matrix to tridiagonal form.
407*
408 indd = 1
409 inde = indd + n
410 indwrk = inde + n
411 IF( wantz ) THEN
412 vect = 'U'
413 ELSE
414 vect = 'N'
415 END IF
416 CALL ssbtrd( vect, uplo, n, ka, ab, ldab, work( indd ),
417 $ work( inde ), q, ldq, work( indwrk ), iinfo )
418*
419* If all eigenvalues are desired and ABSTOL is less than or equal
420* to zero, then call SSTERF or SSTEQR. If this fails for some
421* eigenvalue, then try SSTEBZ.
422*
423 test = .false.
424 IF( indeig ) THEN
425 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
426 test = .true.
427 END IF
428 END IF
429 IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
430 CALL scopy( n, work( indd ), 1, w, 1 )
431 indee = indwrk + 2*n
432 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
433 IF( .NOT.wantz ) THEN
434 CALL ssterf( n, w, work( indee ), info )
435 ELSE
436 CALL slacpy( 'A', n, n, q, ldq, z, ldz )
437 CALL ssteqr( jobz, n, w, work( indee ), z, ldz,
438 $ work( indwrk ), info )
439 IF( info.EQ.0 ) THEN
440 DO 10 i = 1, n
441 ifail( i ) = 0
442 10 CONTINUE
443 END IF
444 END IF
445 IF( info.EQ.0 ) THEN
446 m = n
447 GO TO 30
448 END IF
449 info = 0
450 END IF
451*
452* Otherwise, call SSTEBZ and, if eigenvectors are desired,
453* call SSTEIN.
454*
455 IF( wantz ) THEN
456 order = 'B'
457 ELSE
458 order = 'E'
459 END IF
460 indibl = 1
461 indisp = indibl + n
462 indiwo = indisp + n
463 CALL sstebz( range, order, n, vl, vu, il, iu, abstol,
464 $ work( indd ), work( inde ), m, nsplit, w,
465 $ iwork( indibl ), iwork( indisp ), work( indwrk ),
466 $ iwork( indiwo ), info )
467*
468 IF( wantz ) THEN
469 CALL sstein( n, work( indd ), work( inde ), m, w,
470 $ iwork( indibl ), iwork( indisp ), z, ldz,
471 $ work( indwrk ), iwork( indiwo ), ifail, info )
472*
473* Apply transformation matrix used in reduction to tridiagonal
474* form to eigenvectors returned by SSTEIN.
475*
476 DO 20 j = 1, m
477 CALL scopy( n, z( 1, j ), 1, work( 1 ), 1 )
478 CALL sgemv( 'N', n, n, one, q, ldq, work, 1, zero,
479 $ z( 1, j ), 1 )
480 20 CONTINUE
481 END IF
482*
483 30 CONTINUE
484*
485* If eigenvalues are not in order, then sort them, along with
486* eigenvectors.
487*
488 IF( wantz ) THEN
489 DO 50 j = 1, m - 1
490 i = 0
491 tmp1 = w( j )
492 DO 40 jj = j + 1, m
493 IF( w( jj ).LT.tmp1 ) THEN
494 i = jj
495 tmp1 = w( jj )
496 END IF
497 40 CONTINUE
498*
499 IF( i.NE.0 ) THEN
500 itmp1 = iwork( indibl+i-1 )
501 w( i ) = w( j )
502 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
503 w( j ) = tmp1
504 iwork( indibl+j-1 ) = itmp1
505 CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
506 IF( info.NE.0 ) THEN
507 itmp1 = ifail( i )
508 ifail( i ) = ifail( j )
509 ifail( j ) = itmp1
510 END IF
511 END IF
512 50 CONTINUE
513 END IF
514*
515 RETURN
516*
517* End of SSBGVX
518*

◆ sspev()

subroutine sspev ( character jobz,
character uplo,
integer n,
real, dimension( * ) ap,
real, dimension( * ) w,
real, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer info )

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

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

Purpose:
!>
!> SSPEV computes all the eigenvalues and, optionally, eigenvectors of a
!> real symmetric matrix A in packed storage.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[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]AP
!>          AP is REAL array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the symmetric matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>
!>          On exit, AP is overwritten by values generated during the
!>          reduction to tridiagonal form.  If UPLO = 'U', the diagonal
!>          and first superdiagonal of the tridiagonal matrix T overwrite
!>          the corresponding elements of A, and if UPLO = 'L', the
!>          diagonal and first subdiagonal of T overwrite the
!>          corresponding elements of A.
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
!>          eigenvectors of the matrix A, with the i-th column of Z
!>          holding the eigenvector associated with W(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>
[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 (3*N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  if INFO = i, the algorithm failed to converge; i
!>                off-diagonal elements of an intermediate tridiagonal
!>                form did not converge to zero.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 129 of file sspev.f.

130*
131* -- LAPACK driver routine --
132* -- LAPACK is a software package provided by Univ. of Tennessee, --
133* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
134*
135* .. Scalar Arguments ..
136 CHARACTER JOBZ, UPLO
137 INTEGER INFO, LDZ, N
138* ..
139* .. Array Arguments ..
140 REAL AP( * ), W( * ), WORK( * ), Z( LDZ, * )
141* ..
142*
143* =====================================================================
144*
145* .. Parameters ..
146 REAL ZERO, ONE
147 parameter( zero = 0.0e0, one = 1.0e0 )
148* ..
149* .. Local Scalars ..
150 LOGICAL WANTZ
151 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE
152 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
153 $ SMLNUM
154* ..
155* .. External Functions ..
156 LOGICAL LSAME
157 REAL SLAMCH, SLANSP
158 EXTERNAL lsame, slamch, slansp
159* ..
160* .. External Subroutines ..
161 EXTERNAL sopgtr, sscal, ssptrd, ssteqr, ssterf, xerbla
162* ..
163* .. Intrinsic Functions ..
164 INTRINSIC sqrt
165* ..
166* .. Executable Statements ..
167*
168* Test the input parameters.
169*
170 wantz = lsame( jobz, 'V' )
171*
172 info = 0
173 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
174 info = -1
175 ELSE IF( .NOT.( lsame( uplo, 'U' ) .OR. lsame( uplo, 'L' ) ) )
176 $ THEN
177 info = -2
178 ELSE IF( n.LT.0 ) THEN
179 info = -3
180 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
181 info = -7
182 END IF
183*
184 IF( info.NE.0 ) THEN
185 CALL xerbla( 'SSPEV ', -info )
186 RETURN
187 END IF
188*
189* Quick return if possible
190*
191 IF( n.EQ.0 )
192 $ RETURN
193*
194 IF( n.EQ.1 ) THEN
195 w( 1 ) = ap( 1 )
196 IF( wantz )
197 $ z( 1, 1 ) = one
198 RETURN
199 END IF
200*
201* Get machine constants.
202*
203 safmin = slamch( 'Safe minimum' )
204 eps = slamch( 'Precision' )
205 smlnum = safmin / eps
206 bignum = one / smlnum
207 rmin = sqrt( smlnum )
208 rmax = sqrt( bignum )
209*
210* Scale matrix to allowable range, if necessary.
211*
212 anrm = slansp( 'M', uplo, n, ap, work )
213 iscale = 0
214 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
215 iscale = 1
216 sigma = rmin / anrm
217 ELSE IF( anrm.GT.rmax ) THEN
218 iscale = 1
219 sigma = rmax / anrm
220 END IF
221 IF( iscale.EQ.1 ) THEN
222 CALL sscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
223 END IF
224*
225* Call SSPTRD to reduce symmetric packed matrix to tridiagonal form.
226*
227 inde = 1
228 indtau = inde + n
229 CALL ssptrd( uplo, n, ap, w, work( inde ), work( indtau ), iinfo )
230*
231* For eigenvalues only, call SSTERF. For eigenvectors, first call
232* SOPGTR to generate the orthogonal matrix, then call SSTEQR.
233*
234 IF( .NOT.wantz ) THEN
235 CALL ssterf( n, w, work( inde ), info )
236 ELSE
237 indwrk = indtau + n
238 CALL sopgtr( uplo, n, ap, work( indtau ), z, ldz,
239 $ work( indwrk ), iinfo )
240 CALL ssteqr( jobz, n, w, work( inde ), z, ldz, work( indtau ),
241 $ info )
242 END IF
243*
244* If matrix was scaled, then rescale eigenvalues appropriately.
245*
246 IF( iscale.EQ.1 ) THEN
247 IF( info.EQ.0 ) THEN
248 imax = n
249 ELSE
250 imax = info - 1
251 END IF
252 CALL sscal( imax, one / sigma, w, 1 )
253 END IF
254*
255 RETURN
256*
257* End of SSPEV
258*
slansp
real function slansp(norm, uplo, n, ap, work)
SLANSP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slansp.f:114
sopgtr
subroutine sopgtr(uplo, n, ap, tau, q, ldq, work, info)
SOPGTR
Definition sopgtr.f:114
ssptrd
subroutine ssptrd(uplo, n, ap, d, e, tau, info)
SSPTRD
Definition ssptrd.f:150

◆ sspevd()

subroutine sspevd ( character jobz,
character uplo,
integer n,
real, dimension( * ) ap,
real, dimension( * ) w,
real, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

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

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

Purpose:
!>
!> SSPEVD computes all the eigenvalues and, optionally, eigenvectors
!> of a real symmetric matrix A in packed storage. If eigenvectors are
!> desired, it uses a divide and conquer algorithm.
!>
!> The divide and conquer algorithm makes very mild assumptions about
!> floating point arithmetic. It will work on machines with a guard
!> digit in add/subtract, or on those binary machines without guard
!> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
!> Cray-2. It could conceivably fail on hexadecimal or decimal machines
!> without guard digits, but we know of none.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[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]AP
!>          AP is REAL array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the symmetric matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>
!>          On exit, AP is overwritten by values generated during the
!>          reduction to tridiagonal form.  If UPLO = 'U', the diagonal
!>          and first superdiagonal of the tridiagonal matrix T overwrite
!>          the corresponding elements of A, and if UPLO = 'L', the
!>          diagonal and first subdiagonal of T overwrite the
!>          corresponding elements of A.
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
!>          eigenvectors of the matrix A, with the i-th column of Z
!>          holding the eigenvector associated with W(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>
[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 required LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          If N <= 1,               LWORK must be at least 1.
!>          If JOBZ = 'N' and N > 1, LWORK must be at least 2*N.
!>          If JOBZ = 'V' and N > 1, LWORK must be at least
!>                                                 1 + 6*N + N**2.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the required sizes of the WORK and IWORK
!>          arrays, returns these values as the first entries of the WORK
!>          and IWORK arrays, and no error message related to LWORK or
!>          LIWORK is issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          On exit, if INFO = 0, IWORK(1) returns the required LIWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.
!>          If JOBZ  = 'N' or N <= 1, LIWORK must be at least 1.
!>          If JOBZ  = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
!>
!>          If LIWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the required sizes of the WORK and
!>          IWORK arrays, returns these values as the first entries of
!>          the WORK and IWORK arrays, and no error message related to
!>          LWORK or LIWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  if INFO = i, the algorithm failed to converge; i
!>                off-diagonal elements of an intermediate tridiagonal
!>                form did not converge to zero.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 176 of file sspevd.f.

178*
179* -- LAPACK driver routine --
180* -- LAPACK is a software package provided by Univ. of Tennessee, --
181* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
182*
183* .. Scalar Arguments ..
184 CHARACTER JOBZ, UPLO
185 INTEGER INFO, LDZ, LIWORK, LWORK, N
186* ..
187* .. Array Arguments ..
188 INTEGER IWORK( * )
189 REAL AP( * ), W( * ), WORK( * ), Z( LDZ, * )
190* ..
191*
192* =====================================================================
193*
194* .. Parameters ..
195 REAL ZERO, ONE
196 parameter( zero = 0.0e+0, one = 1.0e+0 )
197* ..
198* .. Local Scalars ..
199 LOGICAL LQUERY, WANTZ
200 INTEGER IINFO, INDE, INDTAU, INDWRK, ISCALE, LIWMIN,
201 $ LLWORK, LWMIN
202 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
203 $ SMLNUM
204* ..
205* .. External Functions ..
206 LOGICAL LSAME
207 REAL SLAMCH, SLANSP
208 EXTERNAL lsame, slamch, slansp
209* ..
210* .. External Subroutines ..
211 EXTERNAL sopmtr, sscal, ssptrd, sstedc, ssterf, xerbla
212* ..
213* .. Intrinsic Functions ..
214 INTRINSIC sqrt
215* ..
216* .. Executable Statements ..
217*
218* Test the input parameters.
219*
220 wantz = lsame( jobz, 'V' )
221 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
222*
223 info = 0
224 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
225 info = -1
226 ELSE IF( .NOT.( lsame( uplo, 'U' ) .OR. lsame( uplo, 'L' ) ) )
227 $ THEN
228 info = -2
229 ELSE IF( n.LT.0 ) THEN
230 info = -3
231 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
232 info = -7
233 END IF
234*
235 IF( info.EQ.0 ) THEN
236 IF( n.LE.1 ) THEN
237 liwmin = 1
238 lwmin = 1
239 ELSE
240 IF( wantz ) THEN
241 liwmin = 3 + 5*n
242 lwmin = 1 + 6*n + n**2
243 ELSE
244 liwmin = 1
245 lwmin = 2*n
246 END IF
247 END IF
248 iwork( 1 ) = liwmin
249 work( 1 ) = lwmin
250*
251 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
252 info = -9
253 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
254 info = -11
255 END IF
256 END IF
257*
258 IF( info.NE.0 ) THEN
259 CALL xerbla( 'SSPEVD', -info )
260 RETURN
261 ELSE IF( lquery ) THEN
262 RETURN
263 END IF
264*
265* Quick return if possible
266*
267 IF( n.EQ.0 )
268 $ RETURN
269*
270 IF( n.EQ.1 ) THEN
271 w( 1 ) = ap( 1 )
272 IF( wantz )
273 $ z( 1, 1 ) = one
274 RETURN
275 END IF
276*
277* Get machine constants.
278*
279 safmin = slamch( 'Safe minimum' )
280 eps = slamch( 'Precision' )
281 smlnum = safmin / eps
282 bignum = one / smlnum
283 rmin = sqrt( smlnum )
284 rmax = sqrt( bignum )
285*
286* Scale matrix to allowable range, if necessary.
287*
288 anrm = slansp( 'M', uplo, n, ap, work )
289 iscale = 0
290 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
291 iscale = 1
292 sigma = rmin / anrm
293 ELSE IF( anrm.GT.rmax ) THEN
294 iscale = 1
295 sigma = rmax / anrm
296 END IF
297 IF( iscale.EQ.1 ) THEN
298 CALL sscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
299 END IF
300*
301* Call SSPTRD to reduce symmetric packed matrix to tridiagonal form.
302*
303 inde = 1
304 indtau = inde + n
305 CALL ssptrd( uplo, n, ap, w, work( inde ), work( indtau ), iinfo )
306*
307* For eigenvalues only, call SSTERF. For eigenvectors, first call
308* SSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
309* tridiagonal matrix, then call SOPMTR to multiply it by the
310* Householder transformations represented in AP.
311*
312 IF( .NOT.wantz ) THEN
313 CALL ssterf( n, w, work( inde ), info )
314 ELSE
315 indwrk = indtau + n
316 llwork = lwork - indwrk + 1
317 CALL sstedc( 'I', n, w, work( inde ), z, ldz, work( indwrk ),
318 $ llwork, iwork, liwork, info )
319 CALL sopmtr( 'L', uplo, 'N', n, n, ap, work( indtau ), z, ldz,
320 $ work( indwrk ), iinfo )
321 END IF
322*
323* If matrix was scaled, then rescale eigenvalues appropriately.
324*
325 IF( iscale.EQ.1 )
326 $ CALL sscal( n, one / sigma, w, 1 )
327*
328 work( 1 ) = lwmin
329 iwork( 1 ) = liwmin
330 RETURN
331*
332* End of SSPEVD
333*
sopmtr
subroutine sopmtr(side, uplo, trans, m, n, ap, tau, c, ldc, work, info)
SOPMTR
Definition sopmtr.f:150

◆ sspevx()

subroutine sspevx ( character jobz,
character range,
character uplo,
integer n,
real, dimension( * ) ap,
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, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

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

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

Purpose:
!>
!> SSPEVX computes selected eigenvalues and, optionally, eigenvectors
!> of a real symmetric matrix A in packed storage.  Eigenvalues/vectors
!> 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]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]AP
!>          AP is REAL array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the symmetric matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>
!>          On exit, AP is overwritten by values generated during the
!>          reduction to tridiagonal form.  If UPLO = 'U', the diagonal
!>          and first superdiagonal of the tridiagonal matrix T overwrite
!>          the corresponding elements of A, and if UPLO = 'L', the
!>          diagonal and first subdiagonal of T overwrite the
!>          corresponding elements of A.
!>
[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 AP 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  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)
!>          If INFO = 0, 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 (8*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 231 of file sspevx.f.

234*
235* -- LAPACK driver routine --
236* -- LAPACK is a software package provided by Univ. of Tennessee, --
237* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
238*
239* .. Scalar Arguments ..
240 CHARACTER JOBZ, RANGE, UPLO
241 INTEGER IL, INFO, IU, LDZ, M, N
242 REAL ABSTOL, VL, VU
243* ..
244* .. Array Arguments ..
245 INTEGER IFAIL( * ), IWORK( * )
246 REAL AP( * ), W( * ), WORK( * ), Z( LDZ, * )
247* ..
248*
249* =====================================================================
250*
251* .. Parameters ..
252 REAL ZERO, ONE
253 parameter( zero = 0.0e0, one = 1.0e0 )
254* ..
255* .. Local Scalars ..
256 LOGICAL ALLEIG, INDEIG, TEST, VALEIG, WANTZ
257 CHARACTER ORDER
258 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
259 $ INDISP, INDIWO, INDTAU, INDWRK, ISCALE, ITMP1,
260 $ J, JJ, NSPLIT
261 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
262 $ SIGMA, SMLNUM, TMP1, VLL, VUU
263* ..
264* .. External Functions ..
265 LOGICAL LSAME
266 REAL SLAMCH, SLANSP
267 EXTERNAL lsame, slamch, slansp
268* ..
269* .. External Subroutines ..
270 EXTERNAL scopy, sopgtr, sopmtr, sscal, ssptrd, sstebz,
271 $ sstein, ssteqr, ssterf, sswap, xerbla
272* ..
273* .. Intrinsic Functions ..
274 INTRINSIC max, min, sqrt
275* ..
276* .. Executable Statements ..
277*
278* Test the input parameters.
279*
280 wantz = lsame( jobz, 'V' )
281 alleig = lsame( range, 'A' )
282 valeig = lsame( range, 'V' )
283 indeig = lsame( range, 'I' )
284*
285 info = 0
286 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
287 info = -1
288 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
289 info = -2
290 ELSE IF( .NOT.( lsame( uplo, 'L' ) .OR. lsame( uplo, 'U' ) ) )
291 $ THEN
292 info = -3
293 ELSE IF( n.LT.0 ) THEN
294 info = -4
295 ELSE
296 IF( valeig ) THEN
297 IF( n.GT.0 .AND. vu.LE.vl )
298 $ info = -7
299 ELSE IF( indeig ) THEN
300 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
301 info = -8
302 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
303 info = -9
304 END IF
305 END IF
306 END IF
307 IF( info.EQ.0 ) THEN
308 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
309 $ info = -14
310 END IF
311*
312 IF( info.NE.0 ) THEN
313 CALL xerbla( 'SSPEVX', -info )
314 RETURN
315 END IF
316*
317* Quick return if possible
318*
319 m = 0
320 IF( n.EQ.0 )
321 $ RETURN
322*
323 IF( n.EQ.1 ) THEN
324 IF( alleig .OR. indeig ) THEN
325 m = 1
326 w( 1 ) = ap( 1 )
327 ELSE
328 IF( vl.LT.ap( 1 ) .AND. vu.GE.ap( 1 ) ) THEN
329 m = 1
330 w( 1 ) = ap( 1 )
331 END IF
332 END IF
333 IF( wantz )
334 $ z( 1, 1 ) = one
335 RETURN
336 END IF
337*
338* Get machine constants.
339*
340 safmin = slamch( 'Safe minimum' )
341 eps = slamch( 'Precision' )
342 smlnum = safmin / eps
343 bignum = one / smlnum
344 rmin = sqrt( smlnum )
345 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
346*
347* Scale matrix to allowable range, if necessary.
348*
349 iscale = 0
350 abstll = abstol
351 IF ( valeig ) THEN
352 vll = vl
353 vuu = vu
354 ELSE
355 vll = zero
356 vuu = zero
357 ENDIF
358 anrm = slansp( 'M', uplo, n, ap, work )
359 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
360 iscale = 1
361 sigma = rmin / anrm
362 ELSE IF( anrm.GT.rmax ) THEN
363 iscale = 1
364 sigma = rmax / anrm
365 END IF
366 IF( iscale.EQ.1 ) THEN
367 CALL sscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
368 IF( abstol.GT.0 )
369 $ abstll = abstol*sigma
370 IF( valeig ) THEN
371 vll = vl*sigma
372 vuu = vu*sigma
373 END IF
374 END IF
375*
376* Call SSPTRD to reduce symmetric packed matrix to tridiagonal form.
377*
378 indtau = 1
379 inde = indtau + n
380 indd = inde + n
381 indwrk = indd + n
382 CALL ssptrd( uplo, n, ap, work( indd ), work( inde ),
383 $ work( indtau ), iinfo )
384*
385* If all eigenvalues are desired and ABSTOL is less than or equal
386* to zero, then call SSTERF or SOPGTR and SSTEQR. If this fails
387* for some eigenvalue, then try SSTEBZ.
388*
389 test = .false.
390 IF (indeig) THEN
391 IF (il.EQ.1 .AND. iu.EQ.n) THEN
392 test = .true.
393 END IF
394 END IF
395 IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
396 CALL scopy( n, work( indd ), 1, w, 1 )
397 indee = indwrk + 2*n
398 IF( .NOT.wantz ) THEN
399 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
400 CALL ssterf( n, w, work( indee ), info )
401 ELSE
402 CALL sopgtr( uplo, n, ap, work( indtau ), z, ldz,
403 $ work( indwrk ), iinfo )
404 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
405 CALL ssteqr( jobz, n, w, work( indee ), z, ldz,
406 $ work( indwrk ), info )
407 IF( info.EQ.0 ) THEN
408 DO 10 i = 1, n
409 ifail( i ) = 0
410 10 CONTINUE
411 END IF
412 END IF
413 IF( info.EQ.0 ) THEN
414 m = n
415 GO TO 20
416 END IF
417 info = 0
418 END IF
419*
420* Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
421*
422 IF( wantz ) THEN
423 order = 'B'
424 ELSE
425 order = 'E'
426 END IF
427 indibl = 1
428 indisp = indibl + n
429 indiwo = indisp + n
430 CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
431 $ work( indd ), work( inde ), m, nsplit, w,
432 $ iwork( indibl ), iwork( indisp ), work( indwrk ),
433 $ iwork( indiwo ), info )
434*
435 IF( wantz ) THEN
436 CALL sstein( n, work( indd ), work( inde ), m, w,
437 $ iwork( indibl ), iwork( indisp ), z, ldz,
438 $ work( indwrk ), iwork( indiwo ), ifail, info )
439*
440* Apply orthogonal matrix used in reduction to tridiagonal
441* form to eigenvectors returned by SSTEIN.
442*
443 CALL sopmtr( 'L', uplo, 'N', n, m, ap, work( indtau ), z, ldz,
444 $ work( indwrk ), iinfo )
445 END IF
446*
447* If matrix was scaled, then rescale eigenvalues appropriately.
448*
449 20 CONTINUE
450 IF( iscale.EQ.1 ) THEN
451 IF( info.EQ.0 ) THEN
452 imax = m
453 ELSE
454 imax = info - 1
455 END IF
456 CALL sscal( imax, one / sigma, w, 1 )
457 END IF
458*
459* If eigenvalues are not in order, then sort them, along with
460* eigenvectors.
461*
462 IF( wantz ) THEN
463 DO 40 j = 1, m - 1
464 i = 0
465 tmp1 = w( j )
466 DO 30 jj = j + 1, m
467 IF( w( jj ).LT.tmp1 ) THEN
468 i = jj
469 tmp1 = w( jj )
470 END IF
471 30 CONTINUE
472*
473 IF( i.NE.0 ) THEN
474 itmp1 = iwork( indibl+i-1 )
475 w( i ) = w( j )
476 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
477 w( j ) = tmp1
478 iwork( indibl+j-1 ) = itmp1
479 CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
480 IF( info.NE.0 ) THEN
481 itmp1 = ifail( i )
482 ifail( i ) = ifail( j )
483 ifail( j ) = itmp1
484 END IF
485 END IF
486 40 CONTINUE
487 END IF
488*
489 RETURN
490*
491* End of SSPEVX
492*

◆ sspgv()

subroutine sspgv ( integer itype,
character jobz,
character uplo,
integer n,
real, dimension( * ) ap,
real, dimension( * ) bp,
real, dimension( * ) w,
real, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer info )

SSPGV

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

Purpose:
!>
!> SSPGV computes all the eigenvalues and, optionally, the eigenvectors
!> of a real generalized symmetric-definite eigenproblem, of the form
!> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
!> Here A and B are assumed to be symmetric, stored in packed format,
!> and B is also positive definite.
!>
Parameters
[in]ITYPE
!>          ITYPE is INTEGER
!>          Specifies the problem type to be solved:
!>          = 1:  A*x = (lambda)*B*x
!>          = 2:  A*B*x = (lambda)*x
!>          = 3:  B*A*x = (lambda)*x
!>
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangles of A and B are stored;
!>          = 'L':  Lower triangles of A and B are stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A and B.  N >= 0.
!>
[in,out]AP
!>          AP is REAL array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the symmetric matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>
!>          On exit, the contents of AP are destroyed.
!>
[in,out]BP
!>          BP is REAL array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the symmetric matrix
!>          B, packed columnwise in a linear array.  The j-th column of B
!>          is stored in the array BP as follows:
!>          if UPLO = 'U', BP(i + (j-1)*j/2) = B(i,j) for 1<=i<=j;
!>          if UPLO = 'L', BP(i + (j-1)*(2*n-j)/2) = B(i,j) for j<=i<=n.
!>
!>          On exit, the triangular factor U or L from the Cholesky
!>          factorization B = U**T*U or B = L*L**T, in the same storage
!>          format as B.
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
!>          eigenvectors.  The eigenvectors are normalized as follows:
!>          if ITYPE = 1 or 2, Z**T*B*Z = I;
!>          if ITYPE = 3, Z**T*inv(B)*Z = I.
!>          If JOBZ = 'N', then Z is not referenced.
!>
[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 (3*N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  SPPTRF or SSPEV returned an error code:
!>             <= N:  if INFO = i, SSPEV failed to converge;
!>                    i off-diagonal elements of an intermediate
!>                    tridiagonal form did not converge to zero.
!>             > N:   if INFO = n + i, for 1 <= i <= n, then the leading
!>                    minor of order i of B is not positive definite.
!>                    The factorization of B could not be completed and
!>                    no eigenvalues or eigenvectors were computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 158 of file sspgv.f.

160*
161* -- LAPACK driver routine --
162* -- LAPACK is a software package provided by Univ. of Tennessee, --
163* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
164*
165* .. Scalar Arguments ..
166 CHARACTER JOBZ, UPLO
167 INTEGER INFO, ITYPE, LDZ, N
168* ..
169* .. Array Arguments ..
170 REAL AP( * ), BP( * ), W( * ), WORK( * ),
171 $ Z( LDZ, * )
172* ..
173*
174* =====================================================================
175*
176* .. Local Scalars ..
177 LOGICAL UPPER, WANTZ
178 CHARACTER TRANS
179 INTEGER J, NEIG
180* ..
181* .. External Functions ..
182 LOGICAL LSAME
183 EXTERNAL lsame
184* ..
185* .. External Subroutines ..
186 EXTERNAL spptrf, sspev, sspgst, stpmv, stpsv, xerbla
187* ..
188* .. Executable Statements ..
189*
190* Test the input parameters.
191*
192 wantz = lsame( jobz, 'V' )
193 upper = lsame( uplo, 'U' )
194*
195 info = 0
196 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
197 info = -1
198 ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
199 info = -2
200 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
201 info = -3
202 ELSE IF( n.LT.0 ) THEN
203 info = -4
204 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
205 info = -9
206 END IF
207 IF( info.NE.0 ) THEN
208 CALL xerbla( 'SSPGV ', -info )
209 RETURN
210 END IF
211*
212* Quick return if possible
213*
214 IF( n.EQ.0 )
215 $ RETURN
216*
217* Form a Cholesky factorization of B.
218*
219 CALL spptrf( uplo, n, bp, info )
220 IF( info.NE.0 ) THEN
221 info = n + info
222 RETURN
223 END IF
224*
225* Transform problem to standard eigenvalue problem and solve.
226*
227 CALL sspgst( itype, uplo, n, ap, bp, info )
228 CALL sspev( jobz, uplo, n, ap, w, z, ldz, work, info )
229*
230 IF( wantz ) THEN
231*
232* Backtransform eigenvectors to the original problem.
233*
234 neig = n
235 IF( info.GT.0 )
236 $ neig = info - 1
237 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
238*
239* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
240* backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
241*
242 IF( upper ) THEN
243 trans = 'N'
244 ELSE
245 trans = 'T'
246 END IF
247*
248 DO 10 j = 1, neig
249 CALL stpsv( uplo, trans, 'Non-unit', n, bp, z( 1, j ),
250 $ 1 )
251 10 CONTINUE
252*
253 ELSE IF( itype.EQ.3 ) THEN
254*
255* For B*A*x=(lambda)*x;
256* backtransform eigenvectors: x = L*y or U**T*y
257*
258 IF( upper ) THEN
259 trans = 'T'
260 ELSE
261 trans = 'N'
262 END IF
263*
264 DO 20 j = 1, neig
265 CALL stpmv( uplo, trans, 'Non-unit', n, bp, z( 1, j ),
266 $ 1 )
267 20 CONTINUE
268 END IF
269 END IF
270 RETURN
271*
272* End of SSPGV
273*
spptrf
subroutine spptrf(uplo, n, ap, info)
SPPTRF
Definition spptrf.f:119
sspgst
subroutine sspgst(itype, uplo, n, ap, bp, info)
SSPGST
Definition sspgst.f:113
sspev
subroutine sspev(jobz, uplo, n, ap, w, z, ldz, work, info)
SSPEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
Definition sspev.f:130
stpmv
subroutine stpmv(uplo, trans, diag, n, ap, x, incx)
STPMV
Definition stpmv.f:142
stpsv
subroutine stpsv(uplo, trans, diag, n, ap, x, incx)
STPSV
Definition stpsv.f:144

◆ sspgvd()

subroutine sspgvd ( integer itype,
character jobz,
character uplo,
integer n,
real, dimension( * ) ap,
real, dimension( * ) bp,
real, dimension( * ) w,
real, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

SSPGVD

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

Purpose:
!>
!> SSPGVD computes all the eigenvalues, and optionally, the eigenvectors
!> of a real generalized symmetric-definite eigenproblem, of the form
!> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
!> B are assumed to be symmetric, stored in packed format, and B is also
!> positive definite.
!> If eigenvectors are desired, it uses a divide and conquer algorithm.
!>
!> The divide and conquer algorithm makes very mild assumptions about
!> floating point arithmetic. It will work on machines with a guard
!> digit in add/subtract, or on those binary machines without guard
!> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
!> Cray-2. It could conceivably fail on hexadecimal or decimal machines
!> without guard digits, but we know of none.
!>
Parameters
[in]ITYPE
!>          ITYPE is INTEGER
!>          Specifies the problem type to be solved:
!>          = 1:  A*x = (lambda)*B*x
!>          = 2:  A*B*x = (lambda)*x
!>          = 3:  B*A*x = (lambda)*x
!>
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangles of A and B are stored;
!>          = 'L':  Lower triangles of A and B are stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A and B.  N >= 0.
!>
[in,out]AP
!>          AP is REAL array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the symmetric matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>
!>          On exit, the contents of AP are destroyed.
!>
[in,out]BP
!>          BP is REAL array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the symmetric matrix
!>          B, packed columnwise in a linear array.  The j-th column of B
!>          is stored in the array BP as follows:
!>          if UPLO = 'U', BP(i + (j-1)*j/2) = B(i,j) for 1<=i<=j;
!>          if UPLO = 'L', BP(i + (j-1)*(2*n-j)/2) = B(i,j) for j<=i<=n.
!>
!>          On exit, the triangular factor U or L from the Cholesky
!>          factorization B = U**T*U or B = L*L**T, in the same storage
!>          format as B.
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
!>          eigenvectors.  The eigenvectors are normalized as follows:
!>          if ITYPE = 1 or 2, Z**T*B*Z = I;
!>          if ITYPE = 3, Z**T*inv(B)*Z = I.
!>          If JOBZ = 'N', then Z is not referenced.
!>
[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 required LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          If N <= 1,               LWORK >= 1.
!>          If JOBZ = 'N' and N > 1, LWORK >= 2*N.
!>          If JOBZ = 'V' and N > 1, LWORK >= 1 + 6*N + 2*N**2.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the required sizes of the WORK and IWORK
!>          arrays, returns these values as the first entries of the WORK
!>          and IWORK arrays, and no error message related to LWORK or
!>          LIWORK is issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          On exit, if INFO = 0, IWORK(1) returns the required LIWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.
!>          If JOBZ  = 'N' or N <= 1, LIWORK >= 1.
!>          If JOBZ  = 'V' and N > 1, LIWORK >= 3 + 5*N.
!>
!>          If LIWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the required sizes of the WORK and
!>          IWORK arrays, returns these values as the first entries of
!>          the WORK and IWORK arrays, and no error message related to
!>          LWORK or LIWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  SPPTRF or SSPEVD returned an error code:
!>             <= N:  if INFO = i, SSPEVD failed to converge;
!>                    i off-diagonal elements of an intermediate
!>                    tridiagonal form did not converge to zero;
!>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
!>                    minor of order i of B is not positive definite.
!>                    The factorization of B could not be completed and
!>                    no eigenvalues or eigenvectors were computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA

Definition at line 208 of file sspgvd.f.

210*
211* -- LAPACK driver routine --
212* -- LAPACK is a software package provided by Univ. of Tennessee, --
213* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
214*
215* .. Scalar Arguments ..
216 CHARACTER JOBZ, UPLO
217 INTEGER INFO, ITYPE, LDZ, LIWORK, LWORK, N
218* ..
219* .. Array Arguments ..
220 INTEGER IWORK( * )
221 REAL AP( * ), BP( * ), W( * ), WORK( * ),
222 $ Z( LDZ, * )
223* ..
224*
225* =====================================================================
226*
227* .. Local Scalars ..
228 LOGICAL LQUERY, UPPER, WANTZ
229 CHARACTER TRANS
230 INTEGER J, LIWMIN, LWMIN, NEIG
231* ..
232* .. External Functions ..
233 LOGICAL LSAME
234 EXTERNAL lsame
235* ..
236* .. External Subroutines ..
237 EXTERNAL spptrf, sspevd, sspgst, stpmv, stpsv, xerbla
238* ..
239* .. Intrinsic Functions ..
240 INTRINSIC max, real
241* ..
242* .. Executable Statements ..
243*
244* Test the input parameters.
245*
246 wantz = lsame( jobz, 'V' )
247 upper = lsame( uplo, 'U' )
248 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
249*
250 info = 0
251 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
252 info = -1
253 ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
254 info = -2
255 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
256 info = -3
257 ELSE IF( n.LT.0 ) THEN
258 info = -4
259 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
260 info = -9
261 END IF
262*
263 IF( info.EQ.0 ) THEN
264 IF( n.LE.1 ) THEN
265 liwmin = 1
266 lwmin = 1
267 ELSE
268 IF( wantz ) THEN
269 liwmin = 3 + 5*n
270 lwmin = 1 + 6*n + 2*n**2
271 ELSE
272 liwmin = 1
273 lwmin = 2*n
274 END IF
275 END IF
276 work( 1 ) = lwmin
277 iwork( 1 ) = liwmin
278 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
279 info = -11
280 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
281 info = -13
282 END IF
283 END IF
284*
285 IF( info.NE.0 ) THEN
286 CALL xerbla( 'SSPGVD', -info )
287 RETURN
288 ELSE IF( lquery ) THEN
289 RETURN
290 END IF
291*
292* Quick return if possible
293*
294 IF( n.EQ.0 )
295 $ RETURN
296*
297* Form a Cholesky factorization of BP.
298*
299 CALL spptrf( uplo, n, bp, info )
300 IF( info.NE.0 ) THEN
301 info = n + info
302 RETURN
303 END IF
304*
305* Transform problem to standard eigenvalue problem and solve.
306*
307 CALL sspgst( itype, uplo, n, ap, bp, info )
308 CALL sspevd( jobz, uplo, n, ap, w, z, ldz, work, lwork, iwork,
309 $ liwork, info )
310 lwmin = max( real( lwmin ), real( work( 1 ) ) )
311 liwmin = max( real( liwmin ), real( iwork( 1 ) ) )
312*
313 IF( wantz ) THEN
314*
315* Backtransform eigenvectors to the original problem.
316*
317 neig = n
318 IF( info.GT.0 )
319 $ neig = info - 1
320 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
321*
322* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
323* backtransform eigenvectors: x = inv(L)**T *y or inv(U)*y
324*
325 IF( upper ) THEN
326 trans = 'N'
327 ELSE
328 trans = 'T'
329 END IF
330*
331 DO 10 j = 1, neig
332 CALL stpsv( uplo, trans, 'Non-unit', n, bp, z( 1, j ),
333 $ 1 )
334 10 CONTINUE
335*
336 ELSE IF( itype.EQ.3 ) THEN
337*
338* For B*A*x=(lambda)*x;
339* backtransform eigenvectors: x = L*y or U**T *y
340*
341 IF( upper ) THEN
342 trans = 'T'
343 ELSE
344 trans = 'N'
345 END IF
346*
347 DO 20 j = 1, neig
348 CALL stpmv( uplo, trans, 'Non-unit', n, bp, z( 1, j ),
349 $ 1 )
350 20 CONTINUE
351 END IF
352 END IF
353*
354 work( 1 ) = lwmin
355 iwork( 1 ) = liwmin
356*
357 RETURN
358*
359* End of SSPGVD
360*
sspevd
subroutine sspevd(jobz, uplo, n, ap, w, z, ldz, work, lwork, iwork, liwork, info)
SSPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition sspevd.f:178

◆ sspgvx()

subroutine sspgvx ( integer itype,
character jobz,
character range,
character uplo,
integer n,
real, dimension( * ) ap,
real, dimension( * ) bp,
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, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

SSPGVX

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

Purpose:
!>
!> SSPGVX computes selected eigenvalues, and optionally, eigenvectors
!> of a real generalized symmetric-definite eigenproblem, of the form
!> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A
!> and B are assumed to be symmetric, stored in packed storage, and B
!> is also positive definite.  Eigenvalues and eigenvectors can be
!> selected by specifying either a range of values or a range of indices
!> for the desired eigenvalues.
!>
Parameters
[in]ITYPE
!>          ITYPE is INTEGER
!>          Specifies the problem type to be solved:
!>          = 1:  A*x = (lambda)*B*x
!>          = 2:  A*B*x = (lambda)*x
!>          = 3:  B*A*x = (lambda)*x
!>
[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]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A and B are stored;
!>          = 'L':  Lower triangle of A and B are stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix pencil (A,B).  N >= 0.
!>
[in,out]AP
!>          AP is REAL array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the symmetric matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>
!>          On exit, the contents of AP are destroyed.
!>
[in,out]BP
!>          BP is REAL array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the symmetric matrix
!>          B, packed columnwise in a linear array.  The j-th column of B
!>          is stored in the array BP as follows:
!>          if UPLO = 'U', BP(i + (j-1)*j/2) = B(i,j) for 1<=i<=j;
!>          if UPLO = 'L', BP(i + (j-1)*(2*n-j)/2) = B(i,j) for j<=i<=n.
!>
!>          On exit, the triangular factor U or L from the Cholesky
!>          factorization B = U**T*U or B = L*L**T, in the same storage
!>          format as B.
!>
[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').
!>
[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 = 'N', then Z is not referenced.
!>          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).
!>          The eigenvectors are normalized as follows:
!>          if ITYPE = 1 or 2, Z**T*B*Z = I;
!>          if ITYPE = 3, Z**T*inv(B)*Z = 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.
!>          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 (8*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:  SPPTRF or SSPEVX returned an error code:
!>             <= N:  if INFO = i, SSPEVX failed to converge;
!>                    i eigenvectors failed to converge.  Their indices
!>                    are stored in array IFAIL.
!>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
!>                    minor of order i of B is not positive definite.
!>                    The factorization of B could not be completed and
!>                    no eigenvalues or eigenvectors were computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA

Definition at line 269 of file sspgvx.f.

272*
273* -- LAPACK driver routine --
274* -- LAPACK is a software package provided by Univ. of Tennessee, --
275* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
276*
277* .. Scalar Arguments ..
278 CHARACTER JOBZ, RANGE, UPLO
279 INTEGER IL, INFO, ITYPE, IU, LDZ, M, N
280 REAL ABSTOL, VL, VU
281* ..
282* .. Array Arguments ..
283 INTEGER IFAIL( * ), IWORK( * )
284 REAL AP( * ), BP( * ), W( * ), WORK( * ),
285 $ Z( LDZ, * )
286* ..
287*
288* =====================================================================
289*
290* .. Local Scalars ..
291 LOGICAL ALLEIG, INDEIG, UPPER, VALEIG, WANTZ
292 CHARACTER TRANS
293 INTEGER J
294* ..
295* .. External Functions ..
296 LOGICAL LSAME
297 EXTERNAL lsame
298* ..
299* .. External Subroutines ..
300 EXTERNAL spptrf, sspevx, sspgst, stpmv, stpsv, xerbla
301* ..
302* .. Intrinsic Functions ..
303 INTRINSIC min
304* ..
305* .. Executable Statements ..
306*
307* Test the input parameters.
308*
309 upper = lsame( uplo, 'U' )
310 wantz = lsame( jobz, 'V' )
311 alleig = lsame( range, 'A' )
312 valeig = lsame( range, 'V' )
313 indeig = lsame( range, 'I' )
314*
315 info = 0
316 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
317 info = -1
318 ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
319 info = -2
320 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
321 info = -3
322 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
323 info = -4
324 ELSE IF( n.LT.0 ) THEN
325 info = -5
326 ELSE
327 IF( valeig ) THEN
328 IF( n.GT.0 .AND. vu.LE.vl ) THEN
329 info = -9
330 END IF
331 ELSE IF( indeig ) THEN
332 IF( il.LT.1 ) THEN
333 info = -10
334 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
335 info = -11
336 END IF
337 END IF
338 END IF
339 IF( info.EQ.0 ) THEN
340 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
341 info = -16
342 END IF
343 END IF
344*
345 IF( info.NE.0 ) THEN
346 CALL xerbla( 'SSPGVX', -info )
347 RETURN
348 END IF
349*
350* Quick return if possible
351*
352 m = 0
353 IF( n.EQ.0 )
354 $ RETURN
355*
356* Form a Cholesky factorization of B.
357*
358 CALL spptrf( uplo, n, bp, info )
359 IF( info.NE.0 ) THEN
360 info = n + info
361 RETURN
362 END IF
363*
364* Transform problem to standard eigenvalue problem and solve.
365*
366 CALL sspgst( itype, uplo, n, ap, bp, info )
367 CALL sspevx( jobz, range, uplo, n, ap, vl, vu, il, iu, abstol, m,
368 $ w, z, ldz, work, iwork, ifail, info )
369*
370 IF( wantz ) THEN
371*
372* Backtransform eigenvectors to the original problem.
373*
374 IF( info.GT.0 )
375 $ m = info - 1
376 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
377*
378* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
379* backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
380*
381 IF( upper ) THEN
382 trans = 'N'
383 ELSE
384 trans = 'T'
385 END IF
386*
387 DO 10 j = 1, m
388 CALL stpsv( uplo, trans, 'Non-unit', n, bp, z( 1, j ),
389 $ 1 )
390 10 CONTINUE
391*
392 ELSE IF( itype.EQ.3 ) THEN
393*
394* For B*A*x=(lambda)*x;
395* backtransform eigenvectors: x = L*y or U**T*y
396*
397 IF( upper ) THEN
398 trans = 'T'
399 ELSE
400 trans = 'N'
401 END IF
402*
403 DO 20 j = 1, m
404 CALL stpmv( uplo, trans, 'Non-unit', n, bp, z( 1, j ),
405 $ 1 )
406 20 CONTINUE
407 END IF
408 END IF
409*
410 RETURN
411*
412* End of SSPGVX
413*
sspevx
subroutine sspevx(jobz, range, uplo, n, ap, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info)
SSPEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition sspevx.f:234

◆ sstev()

subroutine sstev ( character jobz,
integer n,
real, dimension( * ) d,
real, dimension( * ) e,
real, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer info )

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

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

Purpose:
!>
!> SSTEV computes all eigenvalues and, optionally, eigenvectors of a
!> real symmetric tridiagonal matrix A.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix.  N >= 0.
!>
[in,out]D
!>          D is REAL array, dimension (N)
!>          On entry, the n diagonal elements of the tridiagonal matrix
!>          A.
!>          On exit, if INFO = 0, the eigenvalues in ascending order.
!>
[in,out]E
!>          E is REAL array, dimension (N-1)
!>          On entry, the (n-1) subdiagonal elements of the tridiagonal
!>          matrix A, stored in elements 1 to N-1 of E.
!>          On exit, the contents of E are destroyed.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
!>          eigenvectors of the matrix A, with the i-th column of Z
!>          holding the eigenvector associated with D(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>
[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,2*N-2))
!>          If JOBZ = 'N', WORK 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, the algorithm failed to converge; i
!>                off-diagonal elements of E did not converge to zero.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 115 of file sstev.f.

116*
117* -- LAPACK driver routine --
118* -- LAPACK is a software package provided by Univ. of Tennessee, --
119* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
120*
121* .. Scalar Arguments ..
122 CHARACTER JOBZ
123 INTEGER INFO, LDZ, N
124* ..
125* .. Array Arguments ..
126 REAL D( * ), E( * ), WORK( * ), Z( LDZ, * )
127* ..
128*
129* =====================================================================
130*
131* .. Parameters ..
132 REAL ZERO, ONE
133 parameter( zero = 0.0e0, one = 1.0e0 )
134* ..
135* .. Local Scalars ..
136 LOGICAL WANTZ
137 INTEGER IMAX, ISCALE
138 REAL BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
139 $ TNRM
140* ..
141* .. External Functions ..
142 LOGICAL LSAME
143 REAL SLAMCH, SLANST
144 EXTERNAL lsame, slamch, slanst
145* ..
146* .. External Subroutines ..
147 EXTERNAL sscal, ssteqr, ssterf, xerbla
148* ..
149* .. Intrinsic Functions ..
150 INTRINSIC sqrt
151* ..
152* .. Executable Statements ..
153*
154* Test the input parameters.
155*
156 wantz = lsame( jobz, 'V' )
157*
158 info = 0
159 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
160 info = -1
161 ELSE IF( n.LT.0 ) THEN
162 info = -2
163 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
164 info = -6
165 END IF
166*
167 IF( info.NE.0 ) THEN
168 CALL xerbla( 'SSTEV ', -info )
169 RETURN
170 END IF
171*
172* Quick return if possible
173*
174 IF( n.EQ.0 )
175 $ RETURN
176*
177 IF( n.EQ.1 ) THEN
178 IF( wantz )
179 $ z( 1, 1 ) = one
180 RETURN
181 END IF
182*
183* Get machine constants.
184*
185 safmin = slamch( 'Safe minimum' )
186 eps = slamch( 'Precision' )
187 smlnum = safmin / eps
188 bignum = one / smlnum
189 rmin = sqrt( smlnum )
190 rmax = sqrt( bignum )
191*
192* Scale matrix to allowable range, if necessary.
193*
194 iscale = 0
195 tnrm = slanst( 'M', n, d, e )
196 IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
197 iscale = 1
198 sigma = rmin / tnrm
199 ELSE IF( tnrm.GT.rmax ) THEN
200 iscale = 1
201 sigma = rmax / tnrm
202 END IF
203 IF( iscale.EQ.1 ) THEN
204 CALL sscal( n, sigma, d, 1 )
205 CALL sscal( n-1, sigma, e( 1 ), 1 )
206 END IF
207*
208* For eigenvalues only, call SSTERF. For eigenvalues and
209* eigenvectors, call SSTEQR.
210*
211 IF( .NOT.wantz ) THEN
212 CALL ssterf( n, d, e, info )
213 ELSE
214 CALL ssteqr( 'I', n, d, e, z, ldz, work, info )
215 END IF
216*
217* If matrix was scaled, then rescale eigenvalues appropriately.
218*
219 IF( iscale.EQ.1 ) THEN
220 IF( info.EQ.0 ) THEN
221 imax = n
222 ELSE
223 imax = info - 1
224 END IF
225 CALL sscal( imax, one / sigma, d, 1 )
226 END IF
227*
228 RETURN
229*
230* End of SSTEV
231*
slanst
real function slanst(norm, n, d, e)
SLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slanst.f:100

◆ sstevd()

subroutine sstevd ( character jobz,
integer n,
real, dimension( * ) d,
real, dimension( * ) e,
real, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

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

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

Purpose:
!>
!> SSTEVD computes all eigenvalues and, optionally, eigenvectors of a
!> real symmetric tridiagonal matrix. If eigenvectors are desired, it
!> uses a divide and conquer algorithm.
!>
!> The divide and conquer algorithm makes very mild assumptions about
!> floating point arithmetic. It will work on machines with a guard
!> digit in add/subtract, or on those binary machines without guard
!> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
!> Cray-2. It could conceivably fail on hexadecimal or decimal machines
!> without guard digits, but we know of none.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix.  N >= 0.
!>
[in,out]D
!>          D is REAL array, dimension (N)
!>          On entry, the n diagonal elements of the tridiagonal matrix
!>          A.
!>          On exit, if INFO = 0, the eigenvalues in ascending order.
!>
[in,out]E
!>          E is REAL array, dimension (N-1)
!>          On entry, the (n-1) subdiagonal elements of the tridiagonal
!>          matrix A, stored in elements 1 to N-1 of E.
!>          On exit, the contents of E are destroyed.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, N)
!>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
!>          eigenvectors of the matrix A, with the i-th column of Z
!>          holding the eigenvector associated with D(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>
[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)
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          If JOBZ  = 'N' or N <= 1 then LWORK must be at least 1.
!>          If JOBZ  = 'V' and N > 1 then LWORK must be at least
!>                         ( 1 + 4*N + N**2 ).
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal sizes of the WORK and IWORK
!>          arrays, returns these values as the first entries of the WORK
!>          and IWORK arrays, and no error message related to LWORK or
!>          LIWORK is issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.
!>          If JOBZ  = 'N' or N <= 1 then LIWORK must be at least 1.
!>          If JOBZ  = 'V' and N > 1 then LIWORK must be at least 3+5*N.
!>
!>          If LIWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal sizes of the WORK and
!>          IWORK arrays, returns these values as the first entries of
!>          the WORK and IWORK arrays, and no error message related to
!>          LWORK or LIWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, the algorithm failed to converge; i
!>                off-diagonal elements of E did not converge to zero.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 161 of file sstevd.f.

163*
164* -- LAPACK driver routine --
165* -- LAPACK is a software package provided by Univ. of Tennessee, --
166* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
167*
168* .. Scalar Arguments ..
169 CHARACTER JOBZ
170 INTEGER INFO, LDZ, LIWORK, LWORK, N
171* ..
172* .. Array Arguments ..
173 INTEGER IWORK( * )
174 REAL D( * ), E( * ), WORK( * ), Z( LDZ, * )
175* ..
176*
177* =====================================================================
178*
179* .. Parameters ..
180 REAL ZERO, ONE
181 parameter( zero = 0.0e0, one = 1.0e0 )
182* ..
183* .. Local Scalars ..
184 LOGICAL LQUERY, WANTZ
185 INTEGER ISCALE, LIWMIN, LWMIN
186 REAL BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
187 $ TNRM
188* ..
189* .. External Functions ..
190 LOGICAL LSAME
191 REAL SLAMCH, SLANST
192 EXTERNAL lsame, slamch, slanst
193* ..
194* .. External Subroutines ..
195 EXTERNAL sscal, sstedc, ssterf, xerbla
196* ..
197* .. Intrinsic Functions ..
198 INTRINSIC sqrt
199* ..
200* .. Executable Statements ..
201*
202* Test the input parameters.
203*
204 wantz = lsame( jobz, 'V' )
205 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
206*
207 info = 0
208 liwmin = 1
209 lwmin = 1
210 IF( n.GT.1 .AND. wantz ) THEN
211 lwmin = 1 + 4*n + n**2
212 liwmin = 3 + 5*n
213 END IF
214*
215 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
216 info = -1
217 ELSE IF( n.LT.0 ) THEN
218 info = -2
219 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
220 info = -6
221 END IF
222*
223 IF( info.EQ.0 ) THEN
224 work( 1 ) = lwmin
225 iwork( 1 ) = liwmin
226*
227 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
228 info = -8
229 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
230 info = -10
231 END IF
232 END IF
233*
234 IF( info.NE.0 ) THEN
235 CALL xerbla( 'SSTEVD', -info )
236 RETURN
237 ELSE IF( lquery ) THEN
238 RETURN
239 END IF
240*
241* Quick return if possible
242*
243 IF( n.EQ.0 )
244 $ RETURN
245*
246 IF( n.EQ.1 ) THEN
247 IF( wantz )
248 $ z( 1, 1 ) = one
249 RETURN
250 END IF
251*
252* Get machine constants.
253*
254 safmin = slamch( 'Safe minimum' )
255 eps = slamch( 'Precision' )
256 smlnum = safmin / eps
257 bignum = one / smlnum
258 rmin = sqrt( smlnum )
259 rmax = sqrt( bignum )
260*
261* Scale matrix to allowable range, if necessary.
262*
263 iscale = 0
264 tnrm = slanst( 'M', n, d, e )
265 IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
266 iscale = 1
267 sigma = rmin / tnrm
268 ELSE IF( tnrm.GT.rmax ) THEN
269 iscale = 1
270 sigma = rmax / tnrm
271 END IF
272 IF( iscale.EQ.1 ) THEN
273 CALL sscal( n, sigma, d, 1 )
274 CALL sscal( n-1, sigma, e( 1 ), 1 )
275 END IF
276*
277* For eigenvalues only, call SSTERF. For eigenvalues and
278* eigenvectors, call SSTEDC.
279*
280 IF( .NOT.wantz ) THEN
281 CALL ssterf( n, d, e, info )
282 ELSE
283 CALL sstedc( 'I', n, d, e, z, ldz, work, lwork, iwork, liwork,
284 $ info )
285 END IF
286*
287* If matrix was scaled, then rescale eigenvalues appropriately.
288*
289 IF( iscale.EQ.1 )
290 $ CALL sscal( n, one / sigma, d, 1 )
291*
292 work( 1 ) = lwmin
293 iwork( 1 ) = liwmin
294*
295 RETURN
296*
297* End of SSTEVD
298*

◆ sstevr()

subroutine sstevr ( character jobz,
character range,
integer n,
real, dimension( * ) d,
real, dimension( * ) e,
real vl,
real vu,
integer il,
integer iu,
real abstol,
integer m,
real, dimension( * ) w,
real, dimension( ldz, * ) z,
integer ldz,
integer, dimension( * ) isuppz,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

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

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

Purpose:
!>
!> SSTEVR computes selected eigenvalues and, optionally, eigenvectors
!> of a real symmetric tridiagonal matrix T.  Eigenvalues and
!> eigenvectors can be selected by specifying either a range of values
!> or a range of indices for the desired eigenvalues.
!>
!> Whenever possible, SSTEVR calls SSTEMR to compute the
!> eigenspectrum using Relatively Robust Representations.  SSTEMR
!> computes eigenvalues by the dqds algorithm, while orthogonal
!> eigenvectors are computed from various  L D L^T representations
!> (also known as Relatively Robust Representations). Gram-Schmidt
!> orthogonalization is avoided as far as possible. More specifically,
!> the various steps of the algorithm are as follows. For the i-th
!> unreduced block of T,
!>    (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T
!>         is a relatively robust representation,
!>    (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high
!>        relative accuracy by the dqds algorithm,
!>    (c) If there is a cluster of close eigenvalues,  sigma_i
!>        close to the cluster, and go to step (a),
!>    (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T,
!>        compute the corresponding eigenvector by forming a
!>        rank-revealing twisted factorization.
!> The desired accuracy of the output can be specified by the input
!> parameter ABSTOL.
!>
!> For more details, see , by Inderjit Dhillon,
!> Computer Science Division Technical Report No. UCB//CSD-97-971,
!> UC Berkeley, May 1997.
!>
!>
!> Note 1 : SSTEVR calls SSTEMR when the full spectrum is requested
!> on machines which conform to the ieee-754 floating point standard.
!> SSTEVR calls SSTEBZ and SSTEIN on non-ieee machines and
!> when partial spectrum requests are made.
!>
!> Normal execution of SSTEMR may create NaNs and infinities and
!> hence may abort due to a floating point exception in environments
!> which do not handle NaNs and infinities in the ieee standard default
!> manner.
!>
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.
!>          For RANGE = 'V' or 'I' and IU - IL < N - 1, SSTEBZ and
!>          SSTEIN are called
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix.  N >= 0.
!>
[in,out]D
!>          D is REAL 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 REAL 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 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.
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
!>          If high relative accuracy is important, set ABSTOL to
!>          SLAMCH( 'Safe minimum' ).  Doing so will guarantee that
!>          eigenvalues are computed to high relative accuracy when
!>          possible in future releases.  The current code does not
!>          make any guarantees about high relative accuracy, but
!>          future releases will. See J. Barlow and J. Demmel,
!>          , LAPACK Working Note #7, for a discussion
!>          of which matrices define their eigenvalues to high relative
!>          accuracy.
!>
[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).
!>          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]ISUPPZ
!>          ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
!>          The support of the eigenvectors in Z, i.e., the indices
!>          indicating the nonzero elements in Z. The i-th eigenvector
!>          is nonzero only in elements ISUPPZ( 2*i-1 ) through
!>          ISUPPZ( 2*i ).
!>          Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
!>
[out]WORK
!>          WORK is REAL array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal (and
!>          minimal) LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.  LWORK >= 20*N.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal sizes of the WORK and IWORK
!>          arrays, returns these values as the first entries of the WORK
!>          and IWORK arrays, and no error message related to LWORK or
!>          LIWORK is issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          On exit, if INFO = 0, IWORK(1) returns the optimal (and
!>          minimal) LIWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.  LIWORK >= 10*N.
!>
!>          If LIWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal sizes of the WORK and
!>          IWORK arrays, returns these values as the first entries of
!>          the WORK and IWORK arrays, and no error message related to
!>          LWORK or LIWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  Internal error
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Inderjit Dhillon, IBM Almaden, USA
Osni Marques, LBNL/NERSC, USA
Ken Stanley, Computer Science Division, University of California at Berkeley, USA
Jason Riedy, Computer Science Division, University of California at Berkeley, USA

Definition at line 303 of file sstevr.f.

306*
307* -- LAPACK driver routine --
308* -- LAPACK is a software package provided by Univ. of Tennessee, --
309* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
310*
311* .. Scalar Arguments ..
312 CHARACTER JOBZ, RANGE
313 INTEGER IL, INFO, IU, LDZ, LIWORK, LWORK, M, N
314 REAL ABSTOL, VL, VU
315* ..
316* .. Array Arguments ..
317 INTEGER ISUPPZ( * ), IWORK( * )
318 REAL D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
319* ..
320*
321* =====================================================================
322*
323* .. Parameters ..
324 REAL ZERO, ONE, TWO
325 parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
326* ..
327* .. Local Scalars ..
328 LOGICAL ALLEIG, INDEIG, TEST, LQUERY, VALEIG, WANTZ,
329 $ TRYRAC
330 CHARACTER ORDER
331 INTEGER I, IEEEOK, IMAX, INDIBL, INDIFL, INDISP,
332 $ INDIWO, ISCALE, J, JJ, LIWMIN, LWMIN, NSPLIT
333 REAL BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
334 $ TMP1, TNRM, VLL, VUU
335* ..
336* .. External Functions ..
337 LOGICAL LSAME
338 INTEGER ILAENV
339 REAL SLAMCH, SLANST
340 EXTERNAL lsame, ilaenv, slamch, slanst
341* ..
342* .. External Subroutines ..
343 EXTERNAL scopy, sscal, sstebz, sstemr, sstein, ssterf,
344 $ sswap, xerbla
345* ..
346* .. Intrinsic Functions ..
347 INTRINSIC max, min, sqrt
348* ..
349* .. Executable Statements ..
350*
351*
352* Test the input parameters.
353*
354 ieeeok = ilaenv( 10, 'SSTEVR', 'N', 1, 2, 3, 4 )
355*
356 wantz = lsame( jobz, 'V' )
357 alleig = lsame( range, 'A' )
358 valeig = lsame( range, 'V' )
359 indeig = lsame( range, 'I' )
360*
361 lquery = ( ( lwork.EQ.-1 ) .OR. ( liwork.EQ.-1 ) )
362 lwmin = max( 1, 20*n )
363 liwmin = max(1, 10*n )
364*
365*
366 info = 0
367 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
368 info = -1
369 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
370 info = -2
371 ELSE IF( n.LT.0 ) THEN
372 info = -3
373 ELSE
374 IF( valeig ) THEN
375 IF( n.GT.0 .AND. vu.LE.vl )
376 $ info = -7
377 ELSE IF( indeig ) THEN
378 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
379 info = -8
380 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
381 info = -9
382 END IF
383 END IF
384 END IF
385 IF( info.EQ.0 ) THEN
386 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
387 info = -14
388 END IF
389 END IF
390*
391 IF( info.EQ.0 ) THEN
392 work( 1 ) = lwmin
393 iwork( 1 ) = liwmin
394*
395 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
396 info = -17
397 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
398 info = -19
399 END IF
400 END IF
401*
402 IF( info.NE.0 ) THEN
403 CALL xerbla( 'SSTEVR', -info )
404 RETURN
405 ELSE IF( lquery ) THEN
406 RETURN
407 END IF
408*
409* Quick return if possible
410*
411 m = 0
412 IF( n.EQ.0 )
413 $ RETURN
414*
415 IF( n.EQ.1 ) THEN
416 IF( alleig .OR. indeig ) THEN
417 m = 1
418 w( 1 ) = d( 1 )
419 ELSE
420 IF( vl.LT.d( 1 ) .AND. vu.GE.d( 1 ) ) THEN
421 m = 1
422 w( 1 ) = d( 1 )
423 END IF
424 END IF
425 IF( wantz )
426 $ z( 1, 1 ) = one
427 RETURN
428 END IF
429*
430* Get machine constants.
431*
432 safmin = slamch( 'Safe minimum' )
433 eps = slamch( 'Precision' )
434 smlnum = safmin / eps
435 bignum = one / smlnum
436 rmin = sqrt( smlnum )
437 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
438*
439*
440* Scale matrix to allowable range, if necessary.
441*
442 iscale = 0
443 IF( valeig ) THEN
444 vll = vl
445 vuu = vu
446 END IF
447*
448 tnrm = slanst( 'M', n, d, e )
449 IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
450 iscale = 1
451 sigma = rmin / tnrm
452 ELSE IF( tnrm.GT.rmax ) THEN
453 iscale = 1
454 sigma = rmax / tnrm
455 END IF
456 IF( iscale.EQ.1 ) THEN
457 CALL sscal( n, sigma, d, 1 )
458 CALL sscal( n-1, sigma, e( 1 ), 1 )
459 IF( valeig ) THEN
460 vll = vl*sigma
461 vuu = vu*sigma
462 END IF
463 END IF
464
465* Initialize indices into workspaces. Note: These indices are used only
466* if SSTERF or SSTEMR fail.
467
468* IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in SSTEBZ and
469* stores the block indices of each of the M<=N eigenvalues.
470 indibl = 1
471* IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in SSTEBZ and
472* stores the starting and finishing indices of each block.
473 indisp = indibl + n
474* IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
475* that corresponding to eigenvectors that fail to converge in
476* SSTEIN. This information is discarded; if any fail, the driver
477* returns INFO > 0.
478 indifl = indisp + n
479* INDIWO is the offset of the remaining integer workspace.
480 indiwo = indisp + n
481*
482* If all eigenvalues are desired, then
483* call SSTERF or SSTEMR. If this fails for some eigenvalue, then
484* try SSTEBZ.
485*
486*
487 test = .false.
488 IF( indeig ) THEN
489 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
490 test = .true.
491 END IF
492 END IF
493 IF( ( alleig .OR. test ) .AND. ieeeok.EQ.1 ) THEN
494 CALL scopy( n-1, e( 1 ), 1, work( 1 ), 1 )
495 IF( .NOT.wantz ) THEN
496 CALL scopy( n, d, 1, w, 1 )
497 CALL ssterf( n, w, work, info )
498 ELSE
499 CALL scopy( n, d, 1, work( n+1 ), 1 )
500 IF (abstol .LE. two*n*eps) THEN
501 tryrac = .true.
502 ELSE
503 tryrac = .false.
504 END IF
505 CALL sstemr( jobz, 'A', n, work( n+1 ), work, vl, vu, il,
506 $ iu, m, w, z, ldz, n, isuppz, tryrac,
507 $ work( 2*n+1 ), lwork-2*n, iwork, liwork, info )
508*
509 END IF
510 IF( info.EQ.0 ) THEN
511 m = n
512 GO TO 10
513 END IF
514 info = 0
515 END IF
516*
517* Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
518*
519 IF( wantz ) THEN
520 order = 'B'
521 ELSE
522 order = 'E'
523 END IF
524
525 CALL sstebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,
526 $ nsplit, w, iwork( indibl ), iwork( indisp ), work,
527 $ iwork( indiwo ), info )
528*
529 IF( wantz ) THEN
530 CALL sstein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),
531 $ z, ldz, work, iwork( indiwo ), iwork( indifl ),
532 $ info )
533 END IF
534*
535* If matrix was scaled, then rescale eigenvalues appropriately.
536*
537 10 CONTINUE
538 IF( iscale.EQ.1 ) THEN
539 IF( info.EQ.0 ) THEN
540 imax = m
541 ELSE
542 imax = info - 1
543 END IF
544 CALL sscal( imax, one / sigma, w, 1 )
545 END IF
546*
547* If eigenvalues are not in order, then sort them, along with
548* eigenvectors.
549*
550 IF( wantz ) THEN
551 DO 30 j = 1, m - 1
552 i = 0
553 tmp1 = w( j )
554 DO 20 jj = j + 1, m
555 IF( w( jj ).LT.tmp1 ) THEN
556 i = jj
557 tmp1 = w( jj )
558 END IF
559 20 CONTINUE
560*
561 IF( i.NE.0 ) THEN
562 w( i ) = w( j )
563 w( j ) = tmp1
564 CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
565 END IF
566 30 CONTINUE
567 END IF
568*
569* Causes problems with tests 19 & 20:
570* IF (wantz .and. INDEIG ) Z( 1,1) = Z(1,1) / 1.002 + .002
571*
572*
573 work( 1 ) = lwmin
574 iwork( 1 ) = liwmin
575 RETURN
576*
577* End of SSTEVR
578*
sstemr
subroutine sstemr(jobz, range, n, d, e, vl, vu, il, iu, m, w, z, ldz, nzc, isuppz, tryrac, work, lwork, iwork, liwork, info)
SSTEMR
Definition sstemr.f:321

◆ sstevx()

subroutine sstevx ( character jobz,
character range,
integer n,
real, dimension( * ) d,
real, dimension( * ) e,
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, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

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

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

Purpose:
!>
!> SSTEVX 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 REAL 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 REAL 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 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.
!>
!>          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  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 (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 REAL 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 225 of file sstevx.f.

227*
228* -- LAPACK driver routine --
229* -- LAPACK is a software package provided by Univ. of Tennessee, --
230* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
231*
232* .. Scalar Arguments ..
233 CHARACTER JOBZ, RANGE
234 INTEGER IL, INFO, IU, LDZ, M, N
235 REAL ABSTOL, VL, VU
236* ..
237* .. Array Arguments ..
238 INTEGER IFAIL( * ), IWORK( * )
239 REAL D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
240* ..
241*
242* =====================================================================
243*
244* .. Parameters ..
245 REAL ZERO, ONE
246 parameter( zero = 0.0e0, one = 1.0e0 )
247* ..
248* .. Local Scalars ..
249 LOGICAL ALLEIG, INDEIG, TEST, VALEIG, WANTZ
250 CHARACTER ORDER
251 INTEGER I, IMAX, INDIBL, INDISP, INDIWO, INDWRK,
252 $ ISCALE, ITMP1, J, JJ, NSPLIT
253 REAL BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
254 $ TMP1, TNRM, VLL, VUU
255* ..
256* .. External Functions ..
257 LOGICAL LSAME
258 REAL SLAMCH, SLANST
259 EXTERNAL lsame, slamch, slanst
260* ..
261* .. External Subroutines ..
262 EXTERNAL scopy, sscal, sstebz, sstein, ssteqr, ssterf,
263 $ sswap, 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( 'SSTEVX', -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 = slamch( 'Safe minimum' )
330 eps = slamch( '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 ENDIF
346 tnrm = slanst( '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 sscal( n, sigma, d, 1 )
356 CALL sscal( 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 SSTERF or SSTEQR. If this fails for some eigenvalue, then
365* try SSTEBZ.
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 scopy( n, d, 1, w, 1 )
375 CALL scopy( n-1, e( 1 ), 1, work( 1 ), 1 )
376 indwrk = n + 1
377 IF( .NOT.wantz ) THEN
378 CALL ssterf( n, w, work, info )
379 ELSE
380 CALL ssteqr( 'I', n, w, work, z, ldz, work( indwrk ), info )
381 IF( info.EQ.0 ) THEN
382 DO 10 i = 1, n
383 ifail( i ) = 0
384 10 CONTINUE
385 END IF
386 END IF
387 IF( info.EQ.0 ) THEN
388 m = n
389 GO TO 20
390 END IF
391 info = 0
392 END IF
393*
394* Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
395*
396 IF( wantz ) THEN
397 order = 'B'
398 ELSE
399 order = 'E'
400 END IF
401 indwrk = 1
402 indibl = 1
403 indisp = indibl + n
404 indiwo = indisp + n
405 CALL sstebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,
406 $ nsplit, w, iwork( indibl ), iwork( indisp ),
407 $ work( indwrk ), iwork( indiwo ), info )
408*
409 IF( wantz ) THEN
410 CALL sstein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),
411 $ z, ldz, work( indwrk ), iwork( indiwo ), ifail,
412 $ info )
413 END IF
414*
415* If matrix was scaled, then rescale eigenvalues appropriately.
416*
417 20 CONTINUE
418 IF( iscale.EQ.1 ) THEN
419 IF( info.EQ.0 ) THEN
420 imax = m
421 ELSE
422 imax = info - 1
423 END IF
424 CALL sscal( imax, one / sigma, w, 1 )
425 END IF
426*
427* If eigenvalues are not in order, then sort them, along with
428* eigenvectors.
429*
430 IF( wantz ) THEN
431 DO 40 j = 1, m - 1
432 i = 0
433 tmp1 = w( j )
434 DO 30 jj = j + 1, m
435 IF( w( jj ).LT.tmp1 ) THEN
436 i = jj
437 tmp1 = w( jj )
438 END IF
439 30 CONTINUE
440*
441 IF( i.NE.0 ) THEN
442 itmp1 = iwork( indibl+i-1 )
443 w( i ) = w( j )
444 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
445 w( j ) = tmp1
446 iwork( indibl+j-1 ) = itmp1
447 CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
448 IF( info.NE.0 ) THEN
449 itmp1 = ifail( i )
450 ifail( i ) = ifail( j )
451 ifail( j ) = itmp1
452 END IF
453 END IF
454 40 CONTINUE
455 END IF
456*
457 RETURN
458*
459* End of SSTEVX
460*