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

Functions

subroutine dbdsvdx (uplo, jobz, range, n, d, e, vl, vu, il, iu, ns, s, z, ldz, work, iwork, info)
 DBDSVDX
subroutine dggglm (n, m, p, a, lda, b, ldb, d, x, y, work, lwork, info)
 DGGGLM
subroutine dsbev (jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, info)
  DSBEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine dsbev_2stage (jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, lwork, info)
  DSBEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine dsbevd (jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, lwork, iwork, liwork, info)
  DSBEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine dsbevd_2stage (jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, lwork, iwork, liwork, info)
  DSBEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine dsbevx (jobz, range, uplo, n, kd, ab, ldab, q, ldq, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info)
  DSBEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine dsbevx_2stage (jobz, range, uplo, n, kd, ab, ldab, q, ldq, vl, vu, il, iu, abstol, m, w, z, ldz, work, lwork, iwork, ifail, info)
  DSBEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine dsbgv (jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, w, z, ldz, work, info)
 DSBGV
subroutine dsbgvd (jobz, uplo, n, ka, kb, ab, ldab, bb, ldbb, w, z, ldz, work, lwork, iwork, liwork, info)
 DSBGVD
subroutine dsbgvx (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)
 DSBGVX
subroutine dspev (jobz, uplo, n, ap, w, z, ldz, work, info)
  DSPEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine dspevd (jobz, uplo, n, ap, w, z, ldz, work, lwork, iwork, liwork, info)
  DSPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine dspevx (jobz, range, uplo, n, ap, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info)
  DSPEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine dspgv (itype, jobz, uplo, n, ap, bp, w, z, ldz, work, info)
 DSPGV
subroutine dspgvd (itype, jobz, uplo, n, ap, bp, w, z, ldz, work, lwork, iwork, liwork, info)
 DSPGVD
subroutine dspgvx (itype, jobz, range, uplo, n, ap, bp, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info)
 DSPGVX
subroutine dstev (jobz, n, d, e, z, ldz, work, info)
  DSTEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine dstevd (jobz, n, d, e, z, ldz, work, lwork, iwork, liwork, info)
  DSTEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine dstevr (jobz, range, n, d, e, vl, vu, il, iu, abstol, m, w, z, ldz, isuppz, work, lwork, iwork, liwork, info)
  DSTEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
subroutine dstevx (jobz, range, n, d, e, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info)
  DSTEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices

Detailed Description

This is the group of double Other Eigenvalue routines

Function Documentation

◆ dbdsvdx()

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

DBDSVDX

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

Purpose:
!>
!>  DBDSVDX 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 ], DBDSVDX 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 DGESVD/DGESDD) or b) compute s,v and reorder
!>  the values (and corresponding vectors). DBDSVDX implements a) by
!>  calling DSTEVX (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 DOUBLE PRECISION array, dimension (N)
!>          The n diagonal elements of the bidiagonal matrix B.
!>
[in]E
!>          E is DOUBLE PRECISION 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 DOUBLE PRECISION
!>          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 DOUBLE PRECISION
!>          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 DOUBLE PRECISION array, dimension (N)
!>          The first NS elements contain the selected singular values in
!>          ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DSTEVX. The indices of the eigenvectors
!>                   (as returned by DSTEVX) 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 dbdsvdx.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 DOUBLE PRECISION VL, VU
235* ..
236* .. Array Arguments ..
237 INTEGER IWORK( * )
238 DOUBLE PRECISION D( * ), E( * ), S( * ), WORK( * ),
239 $ Z( LDZ, * )
240* ..
241*
242* =====================================================================
243*
244* .. Parameters ..
245 DOUBLE PRECISION ZERO, ONE, TEN, HNDRD, MEIGTH
246 parameter( zero = 0.0d0, one = 1.0d0, ten = 10.0d0,
247 $ hndrd = 100.0d0, meigth = -0.1250d0 )
248 DOUBLE PRECISION FUDGE
249 parameter( fudge = 2.0d0 )
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 DOUBLE PRECISION 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 IDAMAX
265 DOUBLE PRECISION DDOT, DLAMCH, DNRM2
266 EXTERNAL idamax, lsame, daxpy, ddot, dlamch, dnrm2
267* ..
268* .. External Subroutines ..
269 EXTERNAL dstevx, dcopy, dlaset, dscal, dswap, xerbla
270* ..
271* .. Intrinsic Functions ..
272 INTRINSIC abs, dble, 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( 'DBDSVDX', -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*dlamch( 'Safe Minimum' )
340 ulp = dlamch( 'Precision' )
341 eps = dlamch( 'Epsilon' )
342 sqrt2 = sqrt( 2.0d0 )
343 ortol = sqrt( ulp )
344*
345* Criterion for splitting is taken from DBDSQR 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 = idamax( n, d, 1 )
355 smax = abs( d( i ) )
356 i = idamax( 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( dble( 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 DSTEVX.
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 DSTEVX 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 dlaset( '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 dcopy( n, d, 1, work( ietgk ), 2 )
418 CALL dcopy( n-1, e, 1, work( ietgk+1 ), 2 )
419 CALL dstevx( '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 dlaset( '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 DSTEBZ, 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 dcopy( n, d, 1, work( ietgk ), 2 )
442 CALL dcopy( n-1, e, 1, work( ietgk+1 ), 2 )
443 CALL dstevx( '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 dcopy( n, d, 1, work( ietgk ), 2 )
450 CALL dcopy( n-1, e, 1, work( ietgk+1 ), 2 )
451 CALL dstevx( '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, DSTEVX 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 dlaset( '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 dcopy( n, d, 1, work( ietgk ), 2 )
491 CALL dcopy( 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 DSTEVX:
587* WORK( ITEMP: ): 2*5*NTGK
588* IWORK( 1: ): 2*6*NTGK
589*
590 CALL dstevx( 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 DSTEVX.
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 = dnrm2( nru, z( irowu, icolz+i ), 2 )
635 IF( nrmu.EQ.zero ) THEN
636 info = n*2 + 1
637 RETURN
638 END IF
639 CALL dscal( 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 = -ddot( nru, z( irowu, icolz+j ),
646 $ 2, z( irowu, icolz+i ), 2 )
647 CALL daxpy( nru, zjtji,
648 $ z( irowu, icolz+j ), 2,
649 $ z( irowu, icolz+i ), 2 )
650 END DO
651 nrmu = dnrm2( nru, z( irowu, icolz+i ), 2 )
652 CALL dscal( 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 = dnrm2( nrv, z( irowv, icolz+i ), 2 )
658 IF( nrmv.EQ.zero ) THEN
659 info = n*2 + 1
660 RETURN
661 END IF
662 CALL dscal( 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 = -ddot( nrv, z( irowv, icolz+j ),
669 $ 2, z( irowv, icolz+i ), 2 )
670 CALL daxpy( nru, zjtji,
671 $ z( irowv, icolz+j ), 2,
672 $ z( irowv, icolz+i ), 2 )
673 END DO
674 nrmv = dnrm2( nrv, z( irowv, icolz+i ), 2 )
675 CALL dscal( 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 dswap( 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 dcopy( n*2, z( 1,i ), 1, work, 1 )
775 IF( lower ) THEN
776 CALL dcopy( n, work( 2 ), 2, z( n+1,i ), 1 )
777 CALL dcopy( n, work( 1 ), 2, z( 1 ,i ), 1 )
778 ELSE
779 CALL dcopy( n, work( 2 ), 2, z( 1 ,i ), 1 )
780 CALL dcopy( 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 DBDSVDX
788*
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
dstevx
subroutine dstevx(jobz, range, n, d, e, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info)
DSTEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition dstevx.f:227
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
ddot
double precision function ddot(n, dx, incx, dy, incy)
DDOT
Definition ddot.f:82
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dnrm2
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition dnrm2.f90:89
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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)

◆ dggglm()

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

DGGGLM

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

Purpose:
!>
!> DGGGLM 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
!>          On entry, D is the left hand side of the GLM equation.
!>          On exit, D is destroyed.
!>
[out]X
!>          X is DOUBLE PRECISION array, dimension (M)
!>
[out]Y
!>          Y is DOUBLE PRECISION array, dimension (P)
!>
!>          On exit, X and Y are the solutions of the GLM problem.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION 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
!>          DGEQRF, SGERQF, DORMQR 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 dggglm.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 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
195 $ X( * ), Y( * )
196* ..
197*
198* ===================================================================
199*
200* .. Parameters ..
201 DOUBLE PRECISION ZERO, ONE
202 parameter( zero = 0.0d+0, one = 1.0d+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 dcopy, dgemv, dggqrf, dormqr, dormrq, dtrtrs,
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, 'DGEQRF', ' ', n, m, -1, -1 )
247 nb2 = ilaenv( 1, 'DGERQF', ' ', n, m, -1, -1 )
248 nb3 = ilaenv( 1, 'DORMQR', ' ', n, m, p, -1 )
249 nb4 = ilaenv( 1, 'DORMRQ', ' ', 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( 'DGGGLM', -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 dggqrf( 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 dormqr( '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 dtrtrs( '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 dcopy( 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 dgemv( '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 dtrtrs( '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 dcopy( m, d, 1, x, 1 )
339 END IF
340*
341* Backward transformation y = Z**T *y
342*
343 CALL dormrq( '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 DGGGLM
351*
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
dggqrf
subroutine dggqrf(n, m, p, a, lda, taua, b, ldb, taub, work, lwork, info)
DGGQRF
Definition dggqrf.f:215
dormqr
subroutine dormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMQR
Definition dormqr.f:167
dtrtrs
subroutine dtrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
DTRTRS
Definition dtrtrs.f:140
dormrq
subroutine dormrq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMRQ
Definition dormrq.f:167
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:156

◆ dsbev()

subroutine dsbev ( character jobz,
character uplo,
integer n,
integer kd,
double precision, dimension( ldab, * ) ab,
integer ldab,
double precision, dimension( * ) w,
double precision, dimension( ldz, * ) z,
integer ldz,
double precision, dimension( * ) work,
integer info )

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

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

Purpose:
!>
!> DSBEV 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 DOUBLE PRECISION array, dimension (LDAB, N)
!>          On entry, the upper or lower triangle of the symmetric band
!>          matrix A, stored in the first KD+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
!>
!>          On exit, AB is overwritten by values generated during the
!>          reduction to tridiagonal form.  If UPLO = 'U', the first
!>          superdiagonal and the diagonal of the tridiagonal matrix T
!>          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
!>          the diagonal and first subdiagonal of T are returned in the
!>          first two rows of AB.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KD + 1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION 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 DOUBLE PRECISION 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 dsbev.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 DOUBLE PRECISION AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
157* ..
158*
159* =====================================================================
160*
161* .. Parameters ..
162 DOUBLE PRECISION ZERO, ONE
163 parameter( zero = 0.0d0, one = 1.0d0 )
164* ..
165* .. Local Scalars ..
166 LOGICAL LOWER, WANTZ
167 INTEGER IINFO, IMAX, INDE, INDWRK, ISCALE
168 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
169 $ SMLNUM
170* ..
171* .. External Functions ..
172 LOGICAL LSAME
173 DOUBLE PRECISION DLAMCH, DLANSB
174 EXTERNAL lsame, dlamch, dlansb
175* ..
176* .. External Subroutines ..
177 EXTERNAL dlascl, dsbtrd, dscal, dsteqr, dsterf, 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( 'DSBEV ', -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 = dlamch( 'Safe minimum' )
228 eps = dlamch( '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 = dlansb( '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 dlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
248 ELSE
249 CALL dlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
250 END IF
251 END IF
252*
253* Call DSBTRD to reduce symmetric band matrix to tridiagonal form.
254*
255 inde = 1
256 indwrk = inde + n
257 CALL dsbtrd( jobz, uplo, n, kd, ab, ldab, w, work( inde ), z, ldz,
258 $ work( indwrk ), iinfo )
259*
260* For eigenvalues only, call DSTERF. For eigenvectors, call SSTEQR.
261*
262 IF( .NOT.wantz ) THEN
263 CALL dsterf( n, w, work( inde ), info )
264 ELSE
265 CALL dsteqr( 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 dscal( imax, one / sigma, w, 1 )
278 END IF
279*
280 RETURN
281*
282* End of DSBEV
283*
dlascl
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:143
dsteqr
subroutine dsteqr(compz, n, d, e, z, ldz, work, info)
DSTEQR
Definition dsteqr.f:131
dsterf
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86
dlansb
double precision function dlansb(norm, uplo, n, k, ab, ldab, work)
DLANSB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition dlansb.f:129
dsbtrd
subroutine dsbtrd(vect, uplo, n, kd, ab, ldab, d, e, q, ldq, work, info)
DSBTRD
Definition dsbtrd.f:163

◆ dsbev_2stage()

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

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

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

Purpose:
!>
!> DSBEV_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 DOUBLE PRECISION array, dimension (LDAB, N)
!>          On entry, the upper or lower triangle of the symmetric band
!>          matrix A, stored in the first KD+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
!>
!>          On exit, AB is overwritten by values generated during the
!>          reduction to tridiagonal form.  If UPLO = 'U', the first
!>          superdiagonal and the diagonal of the tridiagonal matrix T
!>          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
!>          the diagonal and first subdiagonal of T are returned in the
!>          first two rows of AB.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KD + 1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION 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 DOUBLE PRECISION 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 dsbev_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 DOUBLE PRECISION AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
217* ..
218*
219* =====================================================================
220*
221* .. Parameters ..
222 DOUBLE PRECISION ZERO, ONE
223 parameter( zero = 0.0d0, one = 1.0d0 )
224* ..
225* .. Local Scalars ..
226 LOGICAL LOWER, WANTZ, LQUERY
227 INTEGER IINFO, IMAX, INDE, INDWRK, ISCALE,
228 $ LLWORK, LWMIN, LHTRD, LWTRD, IB, INDHOUS
229 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
230 $ SMLNUM
231* ..
232* .. External Functions ..
233 LOGICAL LSAME
234 INTEGER ILAENV2STAGE
235 DOUBLE PRECISION DLAMCH, DLANSB
236 EXTERNAL lsame, dlamch, dlansb, ilaenv2stage
237* ..
238* .. External Subroutines ..
239 EXTERNAL dlascl, dscal, dsteqr, dsterf, xerbla,
240 $ dsytrd_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, 'DSYTRD_SB2ST', jobz,
274 $ n, kd, -1, -1 )
275 lhtrd = ilaenv2stage( 3, 'DSYTRD_SB2ST', jobz,
276 $ n, kd, ib, -1 )
277 lwtrd = ilaenv2stage( 4, 'DSYTRD_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( 'DSBEV_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 = dlamch( 'Safe minimum' )
313 eps = dlamch( '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 = dlansb( '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 dlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
333 ELSE
334 CALL dlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
335 END IF
336 END IF
337*
338* Call DSYTRD_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 dsytrd_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 DSTERF. For eigenvectors, call SSTEQR.
350*
351 IF( .NOT.wantz ) THEN
352 CALL dsterf( n, w, work( inde ), info )
353 ELSE
354 CALL dsteqr( 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 dscal( 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 DSBEV_2STAGE
376*
dsytrd_sb2st
subroutine dsytrd_sb2st(stage1, vect, uplo, n, kd, ab, ldab, d, e, hous, lhous, work, lwork, info)
DSYTRD_SB2ST reduces a real symmetric band matrix A to real symmetric tridiagonal form T
Definition dsytrd_sb2st.F:230
ilaenv2stage
integer function ilaenv2stage(ispec, name, opts, n1, n2, n3, n4)
ILAENV2STAGE
Definition ilaenv2stage.f:149

◆ dsbevd()

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

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

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

Purpose:
!>
!> DSBEVD 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 DOUBLE PRECISION array, dimension (LDAB, N)
!>          On entry, the upper or lower triangle of the symmetric band
!>          matrix A, stored in the first KD+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
!>
!>          On exit, AB is overwritten by values generated during the
!>          reduction to tridiagonal form.  If UPLO = 'U', the first
!>          superdiagonal and the diagonal of the tridiagonal matrix T
!>          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
!>          the diagonal and first subdiagonal of T are returned in the
!>          first two rows of AB.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KD + 1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION 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 DOUBLE PRECISION 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 dsbevd.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 DOUBLE PRECISION AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
205* ..
206*
207* =====================================================================
208*
209* .. Parameters ..
210 DOUBLE PRECISION ZERO, ONE
211 parameter( zero = 0.0d+0, one = 1.0d+0 )
212* ..
213* .. Local Scalars ..
214 LOGICAL LOWER, LQUERY, WANTZ
215 INTEGER IINFO, INDE, INDWK2, INDWRK, ISCALE, LIWMIN,
216 $ LLWRK2, LWMIN
217 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
218 $ SMLNUM
219* ..
220* .. External Functions ..
221 LOGICAL LSAME
222 DOUBLE PRECISION DLAMCH, DLANSB
223 EXTERNAL lsame, dlamch, dlansb
224* ..
225* .. External Subroutines ..
226 EXTERNAL dgemm, dlacpy, dlascl, dsbtrd, dscal, dstedc,
227 $ dsterf, 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( 'DSBEVD', -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 = dlamch( 'Safe minimum' )
300 eps = dlamch( '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 = dlansb( '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 dlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
320 ELSE
321 CALL dlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
322 END IF
323 END IF
324*
325* Call DSBTRD 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 dsbtrd( jobz, uplo, n, kd, ab, ldab, w, work( inde ), z, ldz,
332 $ work( indwrk ), iinfo )
333*
334* For eigenvalues only, call DSTERF. For eigenvectors, call SSTEDC.
335*
336 IF( .NOT.wantz ) THEN
337 CALL dsterf( n, w, work( inde ), info )
338 ELSE
339 CALL dstedc( 'I', n, w, work( inde ), work( indwrk ), n,
340 $ work( indwk2 ), llwrk2, iwork, liwork, info )
341 CALL dgemm( 'N', 'N', n, n, n, one, z, ldz, work( indwrk ), n,
342 $ zero, work( indwk2 ), n )
343 CALL dlacpy( '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 dscal( n, one / sigma, w, 1 )
350*
351 work( 1 ) = lwmin
352 iwork( 1 ) = liwmin
353 RETURN
354*
355* End of DSBEVD
356*
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
dstedc
subroutine dstedc(compz, n, d, e, z, ldz, work, lwork, iwork, liwork, info)
DSTEDC
Definition dstedc.f:188
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:187

◆ dsbevd_2stage()

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

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

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

Purpose:
!>
!> DSBEVD_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 DOUBLE PRECISION array, dimension (LDAB, N)
!>          On entry, the upper or lower triangle of the symmetric band
!>          matrix A, stored in the first KD+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
!>
!>          On exit, AB is overwritten by values generated during the
!>          reduction to tridiagonal form.  If UPLO = 'U', the first
!>          superdiagonal and the diagonal of the tridiagonal matrix T
!>          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
!>          the diagonal and first subdiagonal of T are returned in the
!>          first two rows of AB.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KD + 1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION 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 DOUBLE PRECISION 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 dsbevd_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 DOUBLE PRECISION AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
248* ..
249*
250* =====================================================================
251*
252* .. Parameters ..
253 DOUBLE PRECISION ZERO, ONE
254 parameter( zero = 0.0d+0, one = 1.0d+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 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
262 $ SMLNUM
263* ..
264* .. External Functions ..
265 LOGICAL LSAME
266 INTEGER ILAENV2STAGE
267 DOUBLE PRECISION DLAMCH, DLANSB
268 EXTERNAL lsame, dlamch, dlansb, ilaenv2stage
269* ..
270* .. External Subroutines ..
271 EXTERNAL dgemm, dlacpy, dlascl, dscal, dstedc,
272 $ dsterf, xerbla, dsytrd_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, 'DSYTRD_SB2ST', jobz, n, kd, -1, -1 )
291 lhtrd = ilaenv2stage( 3, 'DSYTRD_SB2ST', jobz, n, kd, ib, -1 )
292 lwtrd = ilaenv2stage( 4, 'DSYTRD_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( 'DSBEVD_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 = dlamch( 'Safe minimum' )
348 eps = dlamch( '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 = dlansb( '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 dlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
368 ELSE
369 CALL dlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
370 END IF
371 END IF
372*
373* Call DSYTRD_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 dsytrd_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 DSTERF. For eigenvectors, call SSTEDC.
387*
388 IF( .NOT.wantz ) THEN
389 CALL dsterf( n, w, work( inde ), info )
390 ELSE
391 CALL dstedc( 'I', n, w, work( inde ), work( indwrk ), n,
392 $ work( indwk2 ), llwrk2, iwork, liwork, info )
393 CALL dgemm( 'N', 'N', n, n, n, one, z, ldz, work( indwrk ), n,
394 $ zero, work( indwk2 ), n )
395 CALL dlacpy( '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 dscal( n, one / sigma, w, 1 )
402*
403 work( 1 ) = lwmin
404 iwork( 1 ) = liwmin
405 RETURN
406*
407* End of DSBEVD_2STAGE
408*

◆ dsbevx()

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

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

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

Purpose:
!>
!> DSBEVX 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 DOUBLE PRECISION array, dimension (LDAB, N)
!>          On entry, the upper or lower triangle of the symmetric band
!>          matrix A, stored in the first KD+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
!>
!>          On exit, AB is overwritten by values generated during the
!>          reduction to tridiagonal form.  If UPLO = 'U', the first
!>          superdiagonal and the diagonal of the tridiagonal matrix T
!>          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
!>          the diagonal and first subdiagonal of T are returned in the
!>          first two rows of AB.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KD + 1.
!>
[out]Q
!>          Q is DOUBLE PRECISION array, dimension (LDQ, N)
!>          If JOBZ = 'V', the N-by-N orthogonal matrix used in the
!>                         reduction to tridiagonal form.
!>          If JOBZ = 'N', the array Q is not referenced.
!>
[in]LDQ
!>          LDQ is INTEGER
!>          The leading dimension of the array Q.  If JOBZ = 'V', then
!>          LDQ >= max(1,N).
!>
[in]VL
!>          VL is DOUBLE PRECISION
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing AB to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*DLAMCH('S').
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          The first M elements contain the selected eigenvalues in
!>          ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          If an eigenvector fails to converge, then that column of Z
!>          contains the latest approximation to the eigenvector, and the
!>          index of the eigenvector is returned in IFAIL.
!>          If JOBZ = 'N', then Z is not referenced.
!>          Note: the user must ensure that at least max(1,M) columns are
!>          supplied in the array Z; if RANGE = 'V', the exact value of M
!>          is not known in advance and an upper bound must be used.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (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 dsbevx.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 DOUBLE PRECISION ABSTOL, VL, VU
274* ..
275* .. Array Arguments ..
276 INTEGER IFAIL( * ), IWORK( * )
277 DOUBLE PRECISION AB( LDAB, * ), Q( LDQ, * ), W( * ), WORK( * ),
278 $ Z( LDZ, * )
279* ..
280*
281* =====================================================================
282*
283* .. Parameters ..
284 DOUBLE PRECISION ZERO, ONE
285 parameter( zero = 0.0d0, one = 1.0d0 )
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 DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
294 $ SIGMA, SMLNUM, TMP1, VLL, VUU
295* ..
296* .. External Functions ..
297 LOGICAL LSAME
298 DOUBLE PRECISION DLAMCH, DLANSB
299 EXTERNAL lsame, dlamch, dlansb
300* ..
301* .. External Subroutines ..
302 EXTERNAL dcopy, dgemv, dlacpy, dlascl, dsbtrd, dscal,
303 $ dstebz, dstein, dsteqr, dsterf, dswap, 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( 'DSBEVX', -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 = dlamch( 'Safe minimum' )
383 eps = dlamch( '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 END IF
400 anrm = dlansb( '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 dlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
411 ELSE
412 CALL dlascl( '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 DSBTRD to reduce symmetric band matrix to tridiagonal form.
423*
424 indd = 1
425 inde = indd + n
426 indwrk = inde + n
427 CALL dsbtrd( 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 DSTERF or SSTEQR. If this fails for some
432* eigenvalue, then try DSTEBZ.
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 dcopy( n, work( indd ), 1, w, 1 )
442 indee = indwrk + 2*n
443 IF( .NOT.wantz ) THEN
444 CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
445 CALL dsterf( n, w, work( indee ), info )
446 ELSE
447 CALL dlacpy( 'A', n, n, q, ldq, z, ldz )
448 CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
449 CALL dsteqr( 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 DSTEBZ 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 dstebz( 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 dstein( 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 DSTEIN.
486*
487 DO 20 j = 1, m
488 CALL dcopy( n, z( 1, j ), 1, work( 1 ), 1 )
489 CALL dgemv( '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 dscal( 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 dswap( 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 DSBEVX
539*
dstebz
subroutine dstebz(range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, iwork, info)
DSTEBZ
Definition dstebz.f:273
dstein
subroutine dstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
DSTEIN
Definition dstein.f:174

◆ dsbevx_2stage()

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

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

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

Purpose:
!>
!> DSBEVX_2STAGE computes selected eigenvalues and, optionally, eigenvectors
!> of a real symmetric band matrix A using the 2stage technique for
!> the reduction to tridiagonal. Eigenvalues and eigenvectors can
!> be selected by specifying either a range of values or a range of
!> indices for the desired eigenvalues.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>                  Not available in this release.
!>
[in]RANGE
!>          RANGE is CHARACTER*1
!>          = 'A': all eigenvalues will be found;
!>          = 'V': all eigenvalues in the half-open interval (VL,VU]
!>                 will be found;
!>          = 'I': the IL-th through IU-th eigenvalues will be found.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!>
[in]KD
!>          KD is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
!>
[in,out]AB
!>          AB is DOUBLE PRECISION array, dimension (LDAB, N)
!>          On entry, the upper or lower triangle of the symmetric band
!>          matrix A, stored in the first KD+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
!>
!>          On exit, AB is overwritten by values generated during the
!>          reduction to tridiagonal form.  If UPLO = 'U', the first
!>          superdiagonal and the diagonal of the tridiagonal matrix T
!>          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
!>          the diagonal and first subdiagonal of T are returned in the
!>          first two rows of AB.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KD + 1.
!>
[out]Q
!>          Q is DOUBLE PRECISION array, dimension (LDQ, N)
!>          If JOBZ = 'V', the N-by-N orthogonal matrix used in the
!>                         reduction to tridiagonal form.
!>          If JOBZ = 'N', the array Q is not referenced.
!>
[in]LDQ
!>          LDQ is INTEGER
!>          The leading dimension of the array Q.  If JOBZ = 'V', then
!>          LDQ >= max(1,N).
!>
[in]VL
!>          VL is DOUBLE PRECISION
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing AB to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*DLAMCH('S').
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          The first M elements contain the selected eigenvalues in
!>          ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          If an eigenvector fails to converge, then that column of Z
!>          contains the latest approximation to the eigenvector, and the
!>          index of the eigenvector is returned in IFAIL.
!>          If JOBZ = 'N', then Z is not referenced.
!>          Note: the user must ensure that at least max(1,M) columns are
!>          supplied in the array Z; if RANGE = 'V', the exact value of M
!>          is not known in advance and an upper bound must be used.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (LWORK)
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The length of the array WORK. LWORK >= 1, when N <= 1;
!>          otherwise  
!>          If JOBZ = 'N' and N > 1, LWORK must be queried.
!>                                   LWORK = MAX(1, 7*N, dimension) where
!>                                   dimension = (2KD+1)*N + KD*NTHREADS + 2*N
!>                                   where KD is the size of the band.
!>                                   NTHREADS is the number of threads used when
!>                                   openMP compilation is enabled, otherwise =1.
!>          If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (5*N)
!>
[out]IFAIL
!>          IFAIL is INTEGER array, dimension (N)
!>          If JOBZ = 'V', then if INFO = 0, the first M elements of
!>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
!>          indices of the eigenvectors that failed to converge.
!>          If JOBZ = 'N', then IFAIL is not referenced.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  if INFO = i, then i eigenvectors failed to converge.
!>                Their indices are stored in array IFAIL.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  All details about the 2stage techniques are available in:
!>
!>  Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
!>  Parallel reduction to condensed forms for symmetric eigenvalue problems
!>  using aggregated fine-grained and memory-aware kernels. In Proceedings
!>  of 2011 International Conference for High Performance Computing,
!>  Networking, Storage and Analysis (SC '11), New York, NY, USA,
!>  Article 8 , 11 pages.
!>  http://doi.acm.org/10.1145/2063384.2063394
!>
!>  A. Haidar, J. Kurzak, P. Luszczek, 2013.
!>  An improved parallel singular value algorithm and its implementation 
!>  for multicore hardware, In Proceedings of 2013 International Conference
!>  for High Performance Computing, Networking, Storage and Analysis (SC '13).
!>  Denver, Colorado, USA, 2013.
!>  Article 90, 12 pages.
!>  http://doi.acm.org/10.1145/2503210.2503292
!>
!>  A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
!>  A novel hybrid CPU-GPU generalized eigensolver for electronic structure 
!>  calculations based on fine-grained memory aware tasks.
!>  International Journal of High Performance Computing Applications.
!>  Volume 28 Issue 2, Pages 196-209, May 2014.
!>  http://hpc.sagepub.com/content/28/2/196 
!>
!>

Definition at line 319 of file dsbevx_2stage.f.

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

◆ dsbgv()

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

DSBGV

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

Purpose:
!>
!> DSBGV 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DPBSTF.
!>
[in]LDBB
!>          LDBB is INTEGER
!>          The leading dimension of the array BB.  LDBB >= KB+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DPBSTF
!>                    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 dsbgv.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 DOUBLE PRECISION 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 dpbstf, dsbgst, dsbtrd, dsteqr, dsterf, 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( 'DSBGV ', -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 dpbstf( 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 dsbgst( 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 dsbtrd( vect, uplo, n, ka, ab, ldab, w, work( inde ), z, ldz,
263 $ work( indwrk ), iinfo )
264*
265* For eigenvalues only, call DSTERF. For eigenvectors, call SSTEQR.
266*
267 IF( .NOT.wantz ) THEN
268 CALL dsterf( n, w, work( inde ), info )
269 ELSE
270 CALL dsteqr( jobz, n, w, work( inde ), z, ldz, work( indwrk ),
271 $ info )
272 END IF
273 RETURN
274*
275* End of DSBGV
276*
dsbgst
subroutine dsbgst(vect, uplo, n, ka, kb, ab, ldab, bb, ldbb, x, ldx, work, info)
DSBGST
Definition dsbgst.f:159
dpbstf
subroutine dpbstf(uplo, n, kd, ab, ldab, info)
DPBSTF
Definition dpbstf.f:152

◆ dsbgvd()

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

DSBGVD

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

Purpose:
!>
!> DSBGVD 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DPBSTF.
!>
[in]LDBB
!>          LDBB is INTEGER
!>          The leading dimension of the array BB.  LDBB >= KB+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION 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 DOUBLE PRECISION 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 >= 2*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 DPBSTF
!>                    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 dsbgvd.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 DOUBLE PRECISION AB( LDAB, * ), BB( LDBB, * ), W( * ),
239 $ WORK( * ), Z( LDZ, * )
240* ..
241*
242* =====================================================================
243*
244* .. Parameters ..
245 DOUBLE PRECISION ONE, ZERO
246 parameter( one = 1.0d+0, zero = 0.0d+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 dgemm, dlacpy, dpbstf, dsbgst, dsbtrd, dstedc,
260 $ dsterf, 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( 'DSBGVD', -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 dpbstf( 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 dsbgst( 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 dsbtrd( vect, uplo, n, ka, ab, ldab, w, work( inde ), z, ldz,
348 $ work( indwrk ), iinfo )
349*
350* For eigenvalues only, call DSTERF. For eigenvectors, call SSTEDC.
351*
352 IF( .NOT.wantz ) THEN
353 CALL dsterf( n, w, work( inde ), info )
354 ELSE
355 CALL dstedc( 'I', n, w, work( inde ), work( indwrk ), n,
356 $ work( indwk2 ), llwrk2, iwork, liwork, info )
357 CALL dgemm( 'N', 'N', n, n, n, one, z, ldz, work( indwrk ), n,
358 $ zero, work( indwk2 ), n )
359 CALL dlacpy( '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 DSBGVD
368*

◆ dsbgvx()

subroutine dsbgvx ( character jobz,
character range,
character uplo,
integer n,
integer ka,
integer kb,
double precision, dimension( ldab, * ) ab,
integer ldab,
double precision, dimension( ldbb, * ) bb,
integer ldbb,
double precision, dimension( ldq, * ) q,
integer ldq,
double precision vl,
double precision vu,
integer il,
integer iu,
double precision abstol,
integer m,
double precision, dimension( * ) w,
double precision, dimension( ldz, * ) z,
integer ldz,
double precision, dimension( * ) work,
integer, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

DSBGVX

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

Purpose:
!>
!> DSBGVX 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DPBSTF.
!>
[in]LDBB
!>          LDBB is INTEGER
!>          The leading dimension of the array BB.  LDBB >= KB+1.
!>
[out]Q
!>          Q is DOUBLE PRECISION 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 DOUBLE PRECISION
!>
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing A to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*DLAMCH('S').
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION 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 DOUBLE PRECISION 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:  DPBSTF 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 dsbgvx.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 DOUBLE PRECISION ABSTOL, VL, VU
304* ..
305* .. Array Arguments ..
306 INTEGER IFAIL( * ), IWORK( * )
307 DOUBLE PRECISION AB( LDAB, * ), BB( LDBB, * ), Q( LDQ, * ),
308 $ W( * ), WORK( * ), Z( LDZ, * )
309* ..
310*
311* =====================================================================
312*
313* .. Parameters ..
314 DOUBLE PRECISION ZERO, ONE
315 parameter( zero = 0.0d+0, one = 1.0d+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 DOUBLE PRECISION TMP1
323* ..
324* .. External Functions ..
325 LOGICAL LSAME
326 EXTERNAL lsame
327* ..
328* .. External Subroutines ..
329 EXTERNAL dcopy, dgemv, dlacpy, dpbstf, dsbgst, dsbtrd,
330 $ dstebz, dstein, dsteqr, dsterf, dswap, 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( 'DSBGVX', -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 dpbstf( 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 dsbgst( 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 dsbtrd( 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 DSTERF or SSTEQR. If this fails for some
421* eigenvalue, then try DSTEBZ.
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 dcopy( n, work( indd ), 1, w, 1 )
431 indee = indwrk + 2*n
432 CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
433 IF( .NOT.wantz ) THEN
434 CALL dsterf( n, w, work( indee ), info )
435 ELSE
436 CALL dlacpy( 'A', n, n, q, ldq, z, ldz )
437 CALL dsteqr( 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 DSTEBZ and, if eigenvectors are desired,
453* call DSTEIN.
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 dstebz( 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 dstein( 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 DSTEIN.
475*
476 DO 20 j = 1, m
477 CALL dcopy( n, z( 1, j ), 1, work( 1 ), 1 )
478 CALL dgemv( '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 dswap( 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 DSBGVX
518*

◆ dspev()

subroutine dspev ( character jobz,
character uplo,
integer n,
double precision, dimension( * ) ap,
double precision, dimension( * ) w,
double precision, dimension( ldz, * ) z,
integer ldz,
double precision, dimension( * ) work,
integer info )

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

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

Purpose:
!>
!> DSPEV 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION 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 DOUBLE PRECISION 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 dspev.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 DOUBLE PRECISION AP( * ), W( * ), WORK( * ), Z( LDZ, * )
141* ..
142*
143* =====================================================================
144*
145* .. Parameters ..
146 DOUBLE PRECISION ZERO, ONE
147 parameter( zero = 0.0d0, one = 1.0d0 )
148* ..
149* .. Local Scalars ..
150 LOGICAL WANTZ
151 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE
152 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
153 $ SMLNUM
154* ..
155* .. External Functions ..
156 LOGICAL LSAME
157 DOUBLE PRECISION DLAMCH, DLANSP
158 EXTERNAL lsame, dlamch, dlansp
159* ..
160* .. External Subroutines ..
161 EXTERNAL dopgtr, dscal, dsptrd, dsteqr, dsterf, 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( 'DSPEV ', -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 = dlamch( 'Safe minimum' )
204 eps = dlamch( '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 = dlansp( '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 dscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
223 END IF
224*
225* Call DSPTRD to reduce symmetric packed matrix to tridiagonal form.
226*
227 inde = 1
228 indtau = inde + n
229 CALL dsptrd( uplo, n, ap, w, work( inde ), work( indtau ), iinfo )
230*
231* For eigenvalues only, call DSTERF. For eigenvectors, first call
232* DOPGTR to generate the orthogonal matrix, then call DSTEQR.
233*
234 IF( .NOT.wantz ) THEN
235 CALL dsterf( n, w, work( inde ), info )
236 ELSE
237 indwrk = indtau + n
238 CALL dopgtr( uplo, n, ap, work( indtau ), z, ldz,
239 $ work( indwrk ), iinfo )
240 CALL dsteqr( 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 dscal( imax, one / sigma, w, 1 )
253 END IF
254*
255 RETURN
256*
257* End of DSPEV
258*
dlansp
double precision function dlansp(norm, uplo, n, ap, work)
DLANSP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition dlansp.f:114
dsptrd
subroutine dsptrd(uplo, n, ap, d, e, tau, info)
DSPTRD
Definition dsptrd.f:150
dopgtr
subroutine dopgtr(uplo, n, ap, tau, q, ldq, work, info)
DOPGTR
Definition dopgtr.f:114

◆ dspevd()

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

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

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

Purpose:
!>
!> DSPEVD 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION 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 DOUBLE PRECISION 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 dspevd.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 DOUBLE PRECISION AP( * ), W( * ), WORK( * ), Z( LDZ, * )
190* ..
191*
192* =====================================================================
193*
194* .. Parameters ..
195 DOUBLE PRECISION ZERO, ONE
196 parameter( zero = 0.0d+0, one = 1.0d+0 )
197* ..
198* .. Local Scalars ..
199 LOGICAL LQUERY, WANTZ
200 INTEGER IINFO, INDE, INDTAU, INDWRK, ISCALE, LIWMIN,
201 $ LLWORK, LWMIN
202 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
203 $ SMLNUM
204* ..
205* .. External Functions ..
206 LOGICAL LSAME
207 DOUBLE PRECISION DLAMCH, DLANSP
208 EXTERNAL lsame, dlamch, dlansp
209* ..
210* .. External Subroutines ..
211 EXTERNAL dopmtr, dscal, dsptrd, dstedc, dsterf, 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( 'DSPEVD', -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 = dlamch( 'Safe minimum' )
280 eps = dlamch( '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 = dlansp( '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 dscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
299 END IF
300*
301* Call DSPTRD to reduce symmetric packed matrix to tridiagonal form.
302*
303 inde = 1
304 indtau = inde + n
305 CALL dsptrd( uplo, n, ap, w, work( inde ), work( indtau ), iinfo )
306*
307* For eigenvalues only, call DSTERF. For eigenvectors, first call
308* DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
309* tridiagonal matrix, then call DOPMTR to multiply it by the
310* Householder transformations represented in AP.
311*
312 IF( .NOT.wantz ) THEN
313 CALL dsterf( n, w, work( inde ), info )
314 ELSE
315 indwrk = indtau + n
316 llwork = lwork - indwrk + 1
317 CALL dstedc( 'I', n, w, work( inde ), z, ldz, work( indwrk ),
318 $ llwork, iwork, liwork, info )
319 CALL dopmtr( '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 dscal( n, one / sigma, w, 1 )
327*
328 work( 1 ) = lwmin
329 iwork( 1 ) = liwmin
330 RETURN
331*
332* End of DSPEVD
333*
dopmtr
subroutine dopmtr(side, uplo, trans, m, n, ap, tau, c, ldc, work, info)
DOPMTR
Definition dopmtr.f:150

◆ dspevx()

subroutine dspevx ( character jobz,
character range,
character uplo,
integer n,
double precision, dimension( * ) ap,
double precision vl,
double precision vu,
integer il,
integer iu,
double precision abstol,
integer m,
double precision, dimension( * ) w,
double precision, dimension( ldz, * ) z,
integer ldz,
double precision, dimension( * ) work,
integer, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

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

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

Purpose:
!>
!> DSPEVX 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 DOUBLE PRECISION 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 DOUBLE PRECISION
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing AP to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*DLAMCH('S').
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the selected eigenvalues in ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          If an eigenvector fails to converge, then that column of Z
!>          contains the latest approximation to the eigenvector, and the
!>          index of the eigenvector is returned in IFAIL.
!>          If JOBZ = 'N', then Z is not referenced.
!>          Note: the user must ensure that at least max(1,M) columns are
!>          supplied in the array Z; if RANGE = 'V', the exact value of M
!>          is not known in advance and an upper bound must be used.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (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 dspevx.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 DOUBLE PRECISION ABSTOL, VL, VU
243* ..
244* .. Array Arguments ..
245 INTEGER IFAIL( * ), IWORK( * )
246 DOUBLE PRECISION AP( * ), W( * ), WORK( * ), Z( LDZ, * )
247* ..
248*
249* =====================================================================
250*
251* .. Parameters ..
252 DOUBLE PRECISION ZERO, ONE
253 parameter( zero = 0.0d0, one = 1.0d0 )
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 DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
262 $ SIGMA, SMLNUM, TMP1, VLL, VUU
263* ..
264* .. External Functions ..
265 LOGICAL LSAME
266 DOUBLE PRECISION DLAMCH, DLANSP
267 EXTERNAL lsame, dlamch, dlansp
268* ..
269* .. External Subroutines ..
270 EXTERNAL dcopy, dopgtr, dopmtr, dscal, dsptrd, dstebz,
271 $ dstein, dsteqr, dsterf, dswap, 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( 'DSPEVX', -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 = dlamch( 'Safe minimum' )
341 eps = dlamch( '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 END IF
358 anrm = dlansp( '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 dscal( ( 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 DSPTRD 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 dsptrd( 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 DSTERF or DOPGTR and SSTEQR. If this fails
387* for some eigenvalue, then try DSTEBZ.
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 dcopy( n, work( indd ), 1, w, 1 )
397 indee = indwrk + 2*n
398 IF( .NOT.wantz ) THEN
399 CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
400 CALL dsterf( n, w, work( indee ), info )
401 ELSE
402 CALL dopgtr( uplo, n, ap, work( indtau ), z, ldz,
403 $ work( indwrk ), iinfo )
404 CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
405 CALL dsteqr( 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 DSTEBZ 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 dstebz( 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 dstein( 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 DSTEIN.
442*
443 CALL dopmtr( '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 dscal( 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 dswap( 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 DSPEVX
492*

◆ dspgv()

subroutine dspgv ( integer itype,
character jobz,
character uplo,
integer n,
double precision, dimension( * ) ap,
double precision, dimension( * ) bp,
double precision, dimension( * ) w,
double precision, dimension( ldz, * ) z,
integer ldz,
double precision, dimension( * ) work,
integer info )

DSPGV

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

Purpose:
!>
!> DSPGV 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION 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 DOUBLE PRECISION 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:  DPPTRF or DSPEV returned an error code:
!>             <= N:  if INFO = i, DSPEV 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 dspgv.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 DOUBLE PRECISION 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 dpptrf, dspev, dspgst, dtpmv, dtpsv, 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( 'DSPGV ', -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 dpptrf( 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 dspgst( itype, uplo, n, ap, bp, info )
228 CALL dspev( 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 dtpsv( 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 dtpmv( 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 DSPGV
273*
dspgst
subroutine dspgst(itype, uplo, n, ap, bp, info)
DSPGST
Definition dspgst.f:113
dpptrf
subroutine dpptrf(uplo, n, ap, info)
DPPTRF
Definition dpptrf.f:119
dspev
subroutine dspev(jobz, uplo, n, ap, w, z, ldz, work, info)
DSPEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices
Definition dspev.f:130
dtpsv
subroutine dtpsv(uplo, trans, diag, n, ap, x, incx)
DTPSV
Definition dtpsv.f:144
dtpmv
subroutine dtpmv(uplo, trans, diag, n, ap, x, incx)
DTPMV
Definition dtpmv.f:142

◆ dspgvd()

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

DSPGVD

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

Purpose:
!>
!> DSPGVD 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION 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 DOUBLE PRECISION 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:  DPPTRF or DSPEVD returned an error code:
!>             <= N:  if INFO = i, DSPEVD 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 dspgvd.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 DOUBLE PRECISION 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 dpptrf, dspevd, dspgst, dtpmv, dtpsv, xerbla
238* ..
239* .. Intrinsic Functions ..
240 INTRINSIC dble, max
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( 'DSPGVD', -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 dpptrf( 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 dspgst( itype, uplo, n, ap, bp, info )
308 CALL dspevd( jobz, uplo, n, ap, w, z, ldz, work, lwork, iwork,
309 $ liwork, info )
310 lwmin = max( dble( lwmin ), dble( work( 1 ) ) )
311 liwmin = max( dble( liwmin ), dble( 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 dtpsv( 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 dtpmv( 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 DSPGVD
360*
dspevd
subroutine dspevd(jobz, uplo, n, ap, w, z, ldz, work, lwork, iwork, liwork, info)
DSPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition dspevd.f:178

◆ dspgvx()

subroutine dspgvx ( integer itype,
character jobz,
character range,
character uplo,
integer n,
double precision, dimension( * ) ap,
double precision, dimension( * ) bp,
double precision vl,
double precision vu,
integer il,
integer iu,
double precision abstol,
integer m,
double precision, dimension( * ) w,
double precision, dimension( ldz, * ) z,
integer ldz,
double precision, dimension( * ) work,
integer, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

DSPGVX

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

Purpose:
!>
!> DSPGVX 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
!>
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing A to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*DLAMCH('S').
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          On normal exit, the first M elements contain the selected
!>          eigenvalues in ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
!>          If JOBZ = '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 DOUBLE PRECISION 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:  DPPTRF or DSPEVX returned an error code:
!>             <= N:  if INFO = i, DSPEVX 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 dspgvx.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 DOUBLE PRECISION ABSTOL, VL, VU
281* ..
282* .. Array Arguments ..
283 INTEGER IFAIL( * ), IWORK( * )
284 DOUBLE PRECISION 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 dpptrf, dspevx, dspgst, dtpmv, dtpsv, 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( 'DSPGVX', -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 dpptrf( 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 dspgst( itype, uplo, n, ap, bp, info )
367 CALL dspevx( 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 dtpsv( 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 dtpmv( 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 DSPGVX
413*
dspevx
subroutine dspevx(jobz, range, uplo, n, ap, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info)
DSPEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition dspevx.f:234

◆ dstev()

subroutine dstev ( character jobz,
integer n,
double precision, dimension( * ) d,
double precision, dimension( * ) e,
double precision, dimension( ldz, * ) z,
integer ldz,
double precision, dimension( * ) work,
integer info )

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

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

Purpose:
!>
!> DSTEV 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 dstev.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 DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
127* ..
128*
129* =====================================================================
130*
131* .. Parameters ..
132 DOUBLE PRECISION ZERO, ONE
133 parameter( zero = 0.0d0, one = 1.0d0 )
134* ..
135* .. Local Scalars ..
136 LOGICAL WANTZ
137 INTEGER IMAX, ISCALE
138 DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
139 $ TNRM
140* ..
141* .. External Functions ..
142 LOGICAL LSAME
143 DOUBLE PRECISION DLAMCH, DLANST
144 EXTERNAL lsame, dlamch, dlanst
145* ..
146* .. External Subroutines ..
147 EXTERNAL dscal, dsteqr, dsterf, 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( 'DSTEV ', -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 = dlamch( 'Safe minimum' )
186 eps = dlamch( '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 = dlanst( '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 dscal( n, sigma, d, 1 )
205 CALL dscal( n-1, sigma, e( 1 ), 1 )
206 END IF
207*
208* For eigenvalues only, call DSTERF. For eigenvalues and
209* eigenvectors, call DSTEQR.
210*
211 IF( .NOT.wantz ) THEN
212 CALL dsterf( n, d, e, info )
213 ELSE
214 CALL dsteqr( '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 dscal( imax, one / sigma, d, 1 )
226 END IF
227*
228 RETURN
229*
230* End of DSTEV
231*
dlanst
double precision function dlanst(norm, n, d, e)
DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition dlanst.f:100

◆ dstevd()

subroutine dstevd ( character jobz,
integer n,
double precision, dimension( * ) d,
double precision, dimension( * ) e,
double precision, dimension( ldz, * ) z,
integer ldz,
double precision, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

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

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

Purpose:
!>
!> DSTEVD 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 dstevd.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 DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
175* ..
176*
177* =====================================================================
178*
179* .. Parameters ..
180 DOUBLE PRECISION ZERO, ONE
181 parameter( zero = 0.0d0, one = 1.0d0 )
182* ..
183* .. Local Scalars ..
184 LOGICAL LQUERY, WANTZ
185 INTEGER ISCALE, LIWMIN, LWMIN
186 DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
187 $ TNRM
188* ..
189* .. External Functions ..
190 LOGICAL LSAME
191 DOUBLE PRECISION DLAMCH, DLANST
192 EXTERNAL lsame, dlamch, dlanst
193* ..
194* .. External Subroutines ..
195 EXTERNAL dscal, dstedc, dsterf, 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( 'DSTEVD', -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 = dlamch( 'Safe minimum' )
255 eps = dlamch( '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 = dlanst( '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 dscal( n, sigma, d, 1 )
274 CALL dscal( n-1, sigma, e( 1 ), 1 )
275 END IF
276*
277* For eigenvalues only, call DSTERF. For eigenvalues and
278* eigenvectors, call DSTEDC.
279*
280 IF( .NOT.wantz ) THEN
281 CALL dsterf( n, d, e, info )
282 ELSE
283 CALL dstedc( '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 dscal( n, one / sigma, d, 1 )
291*
292 work( 1 ) = lwmin
293 iwork( 1 ) = liwmin
294*
295 RETURN
296*
297* End of DSTEVD
298*

◆ dstevr()

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

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

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

Purpose:
!>
!> DSTEVR 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, DSTEVR calls DSTEMR to compute the
!> eigenspectrum using Relatively Robust Representations.  DSTEMR
!> 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 : DSTEVR calls DSTEMR when the full spectrum is requested
!> on machines which conform to the ieee-754 floating point standard.
!> DSTEVR calls DSTEBZ and DSTEIN on non-ieee machines and
!> when partial spectrum requests are made.
!>
!> Normal execution of DSTEMR 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, DSTEBZ and
!>          DSTEIN are called
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix.  N >= 0.
!>
[in,out]D
!>          D is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing A to tridiagonal form.
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
!>          If high relative accuracy is important, set ABSTOL to
!>          DLAMCH( '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 DOUBLE PRECISION array, dimension (N)
!>          The first M elements contain the selected eigenvalues in
!>          ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          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 DOUBLE PRECISION 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 >= max(1,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 >= max(1,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

Definition at line 301 of file dstevr.f.

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

◆ dstevx()

subroutine dstevx ( character jobz,
character range,
integer n,
double precision, dimension( * ) d,
double precision, dimension( * ) e,
double precision vl,
double precision vu,
integer il,
integer iu,
double precision abstol,
integer m,
double precision, dimension( * ) w,
double precision, dimension( ldz, * ) z,
integer ldz,
double precision, dimension( * ) work,
integer, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

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

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

Purpose:
!>
!> DSTEVX 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is DOUBLE PRECISION
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less
!>          than or equal to zero, then  EPS*|T|  will be used in
!>          its place, where |T| is the 1-norm of the tridiagonal
!>          matrix.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*DLAMCH('S').
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>          The first M elements contain the selected eigenvalues in
!>          ascending order.
!>
[out]Z
!>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          If an eigenvector fails to converge (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 DOUBLE PRECISION 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 dstevx.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 DOUBLE PRECISION ABSTOL, VL, VU
236* ..
237* .. Array Arguments ..
238 INTEGER IFAIL( * ), IWORK( * )
239 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
240* ..
241*
242* =====================================================================
243*
244* .. Parameters ..
245 DOUBLE PRECISION ZERO, ONE
246 parameter( zero = 0.0d0, one = 1.0d0 )
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 DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
254 $ TMP1, TNRM, VLL, VUU
255* ..
256* .. External Functions ..
257 LOGICAL LSAME
258 DOUBLE PRECISION DLAMCH, DLANST
259 EXTERNAL lsame, dlamch, dlanst
260* ..
261* .. External Subroutines ..
262 EXTERNAL dcopy, dscal, dstebz, dstein, dsteqr, dsterf,
263 $ dswap, 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( 'DSTEVX', -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 = dlamch( 'Safe minimum' )
330 eps = dlamch( '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 END IF
346 tnrm = dlanst( '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 dscal( n, sigma, d, 1 )
356 CALL dscal( 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 DSTERF or SSTEQR. If this fails for some eigenvalue, then
365* try DSTEBZ.
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 dcopy( n, d, 1, w, 1 )
375 CALL dcopy( n-1, e( 1 ), 1, work( 1 ), 1 )
376 indwrk = n + 1
377 IF( .NOT.wantz ) THEN
378 CALL dsterf( n, w, work, info )
379 ELSE
380 CALL dsteqr( '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 DSTEBZ 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 dstebz( 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 dstein( 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 dscal( 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 dswap( 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 DSTEVX
460*