LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
Other Auxiliary Routines
Collaboration diagram for Other Auxiliary Routines:

Functions

double precision function dlamch (CMACH)
 DLAMCH More...
 
double precision function dlamc3 (A, B)
 DLAMC3 More...
 
subroutine dlamc1 (BETA, T, RND, IEEE1)
 DLAMC1 More...
 
subroutine dlamc2 (BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX)
 DLAMC2 More...
 
program __dlamchtst.f__
 DLAMCHTST More...
 
double precision function dsecnd ()
 DSECND Using ETIME More...
 
program __dsecndtst.f__
 DSECNDTST More...
 
subroutine ilaver (VERS_MAJOR, VERS_MINOR, VERS_PATCH)
 ILAVER returns the LAPACK version. More...
 
program __lapack_version.f__
 LAPACK_VERSION More...
 
logical function lsame (CA, CB)
 LSAME More...
 
program __lsametst.f__
 LSAMETST More...
 
real function second ()
 SECOND Using ETIME More...
 
real function slamch (CMACH)
 SLAMCH More...
 
real function slamc3 (A, B)
 SLAMC3 More...
 
subroutine slamc1 (BETA, T, RND, IEEE1)
 SLAMC1 More...
 
subroutine slamc2 (BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX)
 SLAMC2 More...
 
program tstiee
 TSTIEE More...
 
logical function disnan (DIN)
 DISNAN tests input for NaN. More...
 
subroutine dlabad (SMALL, LARGE)
 DLABAD More...
 
subroutine dlacpy (UPLO, M, N, A, LDA, B, LDB)
 DLACPY copies all or part of one two-dimensional array to another. More...
 
subroutine dladiv (A, B, C, D, P, Q)
 DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow. More...
 
subroutine dlae2 (A, B, C, RT1, RT2)
 DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix. More...
 
subroutine dlaebz (IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, NAB, WORK, IWORK, INFO)
 DLAEBZ computes the number of eigenvalues of a real symmetric tridiagonal matrix which are less than or equal to a given value, and performs other tasks required by the routine sstebz. More...
 
subroutine dlaev2 (A, B, C, RT1, RT2, CS1, SN1)
 DLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix. More...
 
subroutine dlagts (JOB, N, A, B, C, D, IN, Y, TOL, INFO)
 DLAGTS solves the system of equations (T-λI)x = y or (T-λI)Tx = y,where T is a general tridiagonal matrix and λ a scalar, using the LU factorization computed by slagtf. More...
 
logical function dlaisnan (DIN1, DIN2)
 DLAISNAN tests input for NaN by comparing two arguments for inequality. More...
 
integer function dlaneg (N, D, LLD, SIGMA, PIVMIN, R)
 DLANEG computes the Sturm count. More...
 
double precision function dlanst (NORM, N, D, E)
 DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric tridiagonal matrix. More...
 
double precision function dlapy2 (X, Y)
 DLAPY2 returns sqrt(x2+y2). More...
 
double precision function dlapy3 (X, Y, Z)
 DLAPY3 returns sqrt(x2+y2+z2). More...
 
subroutine dlarnv (IDIST, ISEED, N, X)
 DLARNV returns a vector of random numbers from a uniform or normal distribution. More...
 
subroutine dlarra (N, D, E, E2, SPLTOL, TNRM, NSPLIT, ISPLIT, INFO)
 DLARRA computes the splitting points with the specified threshold. More...
 
subroutine dlarrb (N, D, LLD, IFIRST, ILAST, RTOL1, RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK, PIVMIN, SPDIAM, TWIST, INFO)
 DLARRB provides limited bisection to locate eigenvalues for more accuracy. More...
 
subroutine dlarrc (JOBT, N, VL, VU, D, E, PIVMIN, EIGCNT, LCNT, RCNT, INFO)
 DLARRC computes the number of eigenvalues of the symmetric tridiagonal matrix. More...
 
subroutine dlarrd (RANGE, ORDER, N, VL, VU, IL, IU, GERS, RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, M, W, WERR, WL, WU, IBLOCK, INDEXW, WORK, IWORK, INFO)
 DLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy. More...
 
subroutine dlarre (RANGE, N, VL, VU, IL, IU, D, E, E2, RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M, W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN, WORK, IWORK, INFO)
 DLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduced block Ti, finds base representations and eigenvalues. More...
 
subroutine dlarrf (N, D, L, LD, CLSTRT, CLEND, W, WGAP, WERR, SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA, DPLUS, LPLUS, WORK, INFO)
 DLARRF finds a new relatively robust representation such that at least one of the eigenvalues is relatively isolated. More...
 
subroutine dlarrj (N, D, E2, IFIRST, ILAST, RTOL, OFFSET, W, WERR, WORK, IWORK, PIVMIN, SPDIAM, INFO)
 DLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T. More...
 
subroutine dlarrk (N, IW, GL, GU, D, E2, PIVMIN, RELTOL, W, WERR, INFO)
 DLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy. More...
 
subroutine dlarrr (N, D, E, INFO)
 DLARRR performs tests to decide whether the symmetric tridiagonal matrix T warrants expensive computations which guarantee high relative accuracy in the eigenvalues. More...
 
subroutine dlartg (F, G, CS, SN, R)
 DLARTG generates a plane rotation with real cosine and real sine. More...
 
subroutine dlartgp (F, G, CS, SN, R)
 DLARTGP generates a plane rotation so that the diagonal is nonnegative. More...
 
subroutine dlaruv (ISEED, N, X)
 DLARUV returns a vector of n random real numbers from a uniform distribution. More...
 
subroutine dlas2 (F, G, H, SSMIN, SSMAX)
 DLAS2 computes singular values of a 2-by-2 triangular matrix. More...
 
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. More...
 
subroutine dlasd0 (N, SQRE, D, E, U, LDU, VT, LDVT, SMLSIZ, IWORK, WORK, INFO)
 DLASD0 computes the singular values of a real upper bidiagonal n-by-m matrix B with diagonal d and off-diagonal e. Used by sbdsdc. More...
 
subroutine dlasd1 (NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT, IDXQ, IWORK, WORK, INFO)
 DLASD1 computes the SVD of an upper bidiagonal matrix B of the specified size. Used by sbdsdc. More...
 
subroutine dlasd2 (NL, NR, SQRE, K, D, Z, ALPHA, BETA, U, LDU, VT, LDVT, DSIGMA, U2, LDU2, VT2, LDVT2, IDXP, IDX, IDXC, IDXQ, COLTYP, INFO)
 DLASD2 merges the two sets of singular values together into a single sorted set. Used by sbdsdc. More...
 
subroutine dlasd3 (NL, NR, SQRE, K, D, Q, LDQ, DSIGMA, U, LDU, U2, LDU2, VT, LDVT, VT2, LDVT2, IDXC, CTOT, Z, INFO)
 DLASD3 finds all square roots of the roots of the secular equation, as defined by the values in D and Z, and then updates the singular vectors by matrix multiplication. Used by sbdsdc. More...
 
subroutine dlasd4 (N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO)
 DLASD4 computes the square root of the i-th updated eigenvalue of a positive symmetric rank-one modification to a positive diagonal matrix. Used by dbdsdc. More...
 
subroutine dlasd5 (I, D, Z, DELTA, RHO, DSIGMA, WORK)
 DLASD5 computes the square root of the i-th eigenvalue of a positive symmetric rank-one modification of a 2-by-2 diagonal matrix. Used by sbdsdc. More...
 
subroutine dlasd6 (ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA, IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK, IWORK, INFO)
 DLASD6 computes the SVD of an updated upper bidiagonal matrix obtained by merging two smaller ones by appending a row. Used by sbdsdc. More...
 
subroutine dlasd7 (ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL, VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S, INFO)
 DLASD7 merges the two sets of singular values together into a single sorted set. Then it tries to deflate the size of the problem. Used by sbdsdc. More...
 
subroutine dlasd8 (ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR, DSIGMA, WORK, INFO)
 DLASD8 finds the square roots of the roots of the secular equation, and stores, for each element in D, the distance to its two nearest poles. Used by sbdsdc. More...
 
subroutine dlasda (ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, IWORK, INFO)
 DLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagonal d and off-diagonal e. Used by sbdsdc. More...
 
subroutine dlasdq (UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
 DLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e. Used by sbdsdc. More...
 
subroutine dlasdt (N, LVL, ND, INODE, NDIML, NDIMR, MSUB)
 DLASDT creates a tree of subproblems for bidiagonal divide and conquer. Used by sbdsdc. More...
 
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. More...
 
subroutine dlasr (SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
 DLASR applies a sequence of plane rotations to a general rectangular matrix. More...
 
subroutine dlassq (N, X, INCX, SCALE, SUMSQ)
 DLASSQ updates a sum of squares represented in scaled form. More...
 
subroutine dlasv2 (F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL)
 DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix. More...
 
integer function ieeeck (ISPEC, ZERO, ONE)
 IEEECK More...
 
integer function iladlc (M, N, A, LDA)
 ILADLC scans a matrix for its last non-zero column. More...
 
integer function iladlr (M, N, A, LDA)
 ILADLR scans a matrix for its last non-zero row. More...
 
integer function ilaenv (ISPEC, NAME, OPTS, N1, N2, N3, N4)
 ILAENV More...
 
integer function iparmq (ISPEC, NAME, OPTS, N, ILO, IHI, LWORK)
 IPARMQ More...
 
logical function lsamen (N, CA, CB)
 LSAMEN More...
 
logical function sisnan (SIN)
 SISNAN tests input for NaN. More...
 
subroutine slabad (SMALL, LARGE)
 SLABAD More...
 
subroutine slacpy (UPLO, M, N, A, LDA, B, LDB)
 SLACPY copies all or part of one two-dimensional array to another. More...
 
subroutine sladiv (A, B, C, D, P, Q)
 SLADIV performs complex division in real arithmetic, avoiding unnecessary overflow. More...
 
subroutine slae2 (A, B, C, RT1, RT2)
 SLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix. More...
 
subroutine slaebz (IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, NAB, WORK, IWORK, INFO)
 SLAEBZ computes the number of eigenvalues of a real symmetric tridiagonal matrix which are less than or equal to a given value, and performs other tasks required by the routine sstebz. More...
 
subroutine slaev2 (A, B, C, RT1, RT2, CS1, SN1)
 SLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix. More...
 
subroutine slag2d (M, N, SA, LDSA, A, LDA, INFO)
 SLAG2D converts a single precision matrix to a double precision matrix. More...
 
subroutine slagts (JOB, N, A, B, C, D, IN, Y, TOL, INFO)
 SLAGTS solves the system of equations (T-λI)x = y or (T-λI)Tx = y,where T is a general tridiagonal matrix and λ a scalar, using the LU factorization computed by slagtf. More...
 
logical function slaisnan (SIN1, SIN2)
 SLAISNAN tests input for NaN by comparing two arguments for inequality. More...
 
integer function slaneg (N, D, LLD, SIGMA, PIVMIN, R)
 SLANEG computes the Sturm count. More...
 
real function slanst (NORM, N, D, E)
 SLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric tridiagonal matrix. More...
 
real function slapy2 (X, Y)
 SLAPY2 returns sqrt(x2+y2). More...
 
real function slapy3 (X, Y, Z)
 SLAPY3 returns sqrt(x2+y2+z2). More...
 
subroutine slarnv (IDIST, ISEED, N, X)
 SLARNV returns a vector of random numbers from a uniform or normal distribution. More...
 
subroutine slarra (N, D, E, E2, SPLTOL, TNRM, NSPLIT, ISPLIT, INFO)
 SLARRA computes the splitting points with the specified threshold. More...
 
subroutine slarrb (N, D, LLD, IFIRST, ILAST, RTOL1, RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK, PIVMIN, SPDIAM, TWIST, INFO)
 SLARRB provides limited bisection to locate eigenvalues for more accuracy. More...
 
subroutine slarrc (JOBT, N, VL, VU, D, E, PIVMIN, EIGCNT, LCNT, RCNT, INFO)
 SLARRC computes the number of eigenvalues of the symmetric tridiagonal matrix. More...
 
subroutine slarrd (RANGE, ORDER, N, VL, VU, IL, IU, GERS, RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, M, W, WERR, WL, WU, IBLOCK, INDEXW, WORK, IWORK, INFO)
 SLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy. More...
 
subroutine slarre (RANGE, N, VL, VU, IL, IU, D, E, E2, RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M, W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN, WORK, IWORK, INFO)
 SLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduced block Ti, finds base representations and eigenvalues. More...
 
subroutine slarrf (N, D, L, LD, CLSTRT, CLEND, W, WGAP, WERR, SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA, DPLUS, LPLUS, WORK, INFO)
 SLARRF finds a new relatively robust representation such that at least one of the eigenvalues is relatively isolated. More...
 
subroutine slarrj (N, D, E2, IFIRST, ILAST, RTOL, OFFSET, W, WERR, WORK, IWORK, PIVMIN, SPDIAM, INFO)
 SLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T. More...
 
subroutine slarrk (N, IW, GL, GU, D, E2, PIVMIN, RELTOL, W, WERR, INFO)
 SLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy. More...
 
subroutine slarrr (N, D, E, INFO)
 SLARRR performs tests to decide whether the symmetric tridiagonal matrix T warrants expensive computations which guarantee high relative accuracy in the eigenvalues. More...
 
subroutine slartg (F, G, CS, SN, R)
 SLARTG generates a plane rotation with real cosine and real sine. More...
 
subroutine slartgp (F, G, CS, SN, R)
 SLARTGP generates a plane rotation so that the diagonal is nonnegative. More...
 
subroutine slaruv (ISEED, N, X)
 SLARUV returns a vector of n random real numbers from a uniform distribution. More...
 
subroutine slas2 (F, G, H, SSMIN, SSMAX)
 SLAS2 computes singular values of a 2-by-2 triangular matrix. More...
 
subroutine slascl (TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
 SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom. More...
 
subroutine slasd0 (N, SQRE, D, E, U, LDU, VT, LDVT, SMLSIZ, IWORK, WORK, INFO)
 SLASD0 computes the singular values of a real upper bidiagonal n-by-m matrix B with diagonal d and off-diagonal e. Used by sbdsdc. More...
 
subroutine slasd1 (NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT, IDXQ, IWORK, WORK, INFO)
 SLASD1 computes the SVD of an upper bidiagonal matrix B of the specified size. Used by sbdsdc. More...
 
subroutine slasd2 (NL, NR, SQRE, K, D, Z, ALPHA, BETA, U, LDU, VT, LDVT, DSIGMA, U2, LDU2, VT2, LDVT2, IDXP, IDX, IDXC, IDXQ, COLTYP, INFO)
 SLASD2 merges the two sets of singular values together into a single sorted set. Used by sbdsdc. More...
 
subroutine slasd3 (NL, NR, SQRE, K, D, Q, LDQ, DSIGMA, U, LDU, U2, LDU2, VT, LDVT, VT2, LDVT2, IDXC, CTOT, Z, INFO)
 SLASD3 finds all square roots of the roots of the secular equation, as defined by the values in D and Z, and then updates the singular vectors by matrix multiplication. Used by sbdsdc. More...
 
subroutine slasd4 (N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO)
 SLASD4 computes the square root of the i-th updated eigenvalue of a positive symmetric rank-one modification to a positive diagonal matrix. Used by sbdsdc. More...
 
subroutine slasd5 (I, D, Z, DELTA, RHO, DSIGMA, WORK)
 SLASD5 computes the square root of the i-th eigenvalue of a positive symmetric rank-one modification of a 2-by-2 diagonal matrix. Used by sbdsdc. More...
 
subroutine slasd6 (ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA, IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK, IWORK, INFO)
 SLASD6 computes the SVD of an updated upper bidiagonal matrix obtained by merging two smaller ones by appending a row. Used by sbdsdc. More...
 
subroutine slasd7 (ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL, VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S, INFO)
 SLASD7 merges the two sets of singular values together into a single sorted set. Then it tries to deflate the size of the problem. Used by sbdsdc. More...
 
subroutine slasd8 (ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR, DSIGMA, WORK, INFO)
 SLASD8 finds the square roots of the roots of the secular equation, and stores, for each element in D, the distance to its two nearest poles. Used by sbdsdc. More...
 
subroutine slasda (ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, IWORK, INFO)
 SLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagonal d and off-diagonal e. Used by sbdsdc. More...
 
subroutine slasdq (UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
 SLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e. Used by sbdsdc. More...
 
subroutine slasdt (N, LVL, ND, INODE, NDIML, NDIMR, MSUB)
 SLASDT creates a tree of subproblems for bidiagonal divide and conquer. Used by sbdsdc. More...
 
subroutine slaset (UPLO, M, N, ALPHA, BETA, A, LDA)
 SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values. More...
 
subroutine slasr (SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
 SLASR applies a sequence of plane rotations to a general rectangular matrix. More...
 
subroutine slassq (N, X, INCX, SCALE, SUMSQ)
 SLASSQ updates a sum of squares represented in scaled form. More...
 
subroutine slasv2 (F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL)
 SLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix. More...
 
subroutine xerbla (SRNAME, INFO)
 XERBLA More...
 
subroutine xerbla_array (SRNAME_ARRAY, SRNAME_LEN, INFO)
 XERBLA_ARRAY More...
 

Detailed Description

This is the group of Other Auxiliary routines

Function Documentation

program __dlamchtst.f__ ( )

DLAMCHTST

Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 35 of file dlamchtst.f.

35  DOUBLE PRECISION base, emax, emin, eps, prec, rmax, rmin, rnd,
36  $ sfmin, t

Here is the call graph for this function:

program __dsecndtst.f__ ( )

DSECNDTST

Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 37 of file dsecndtst.f.

37  INTEGER nmax, its

Here is the call graph for this function:

program __lapack_version.f__ ( )

LAPACK_VERSION

Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 32 of file LAPACK_version.f.

32  INTEGER major, minor, patch

Here is the call graph for this function:

program __lsametst.f__ ( )

LSAMETST

Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 36 of file lsametst.f.

36  INTEGER i1, i2

Here is the call graph for this function:

logical function disnan ( double precision  DIN)

DISNAN tests input for NaN.

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

Purpose:
 DISNAN returns .TRUE. if its argument is NaN, and .FALSE.
 otherwise.  To be replaced by the Fortran 2003 intrinsic in the
 future.
Parameters
[in]DIN
          DIN is DOUBLE PRECISION
          Input to test for NaN.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 61 of file disnan.f.

61 *
62 * -- LAPACK auxiliary routine (version 3.4.2) --
63 * -- LAPACK is a software package provided by Univ. of Tennessee, --
64 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
65 * September 2012
66 *
67 * .. Scalar Arguments ..
68  DOUBLE PRECISION din
69 * ..
70 *
71 * =====================================================================
72 *
73 * .. External Functions ..
74  LOGICAL dlaisnan
75  EXTERNAL dlaisnan
76 * ..
77 * .. Executable Statements ..
78  disnan = dlaisnan(din,din)
79  RETURN
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:61
logical function dlaisnan(DIN1, DIN2)
DLAISNAN tests input for NaN by comparing two arguments for inequality.
Definition: dlaisnan.f:76
subroutine dlabad ( double precision  SMALL,
double precision  LARGE 
)

DLABAD

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

Purpose:
 DLABAD takes as input the values computed by DLAMCH for underflow and
 overflow, and returns the square root of each of these values if the
 log of LARGE is sufficiently large.  This subroutine is intended to
 identify machines with a large exponent range, such as the Crays, and
 redefine the underflow and overflow limits to be the square roots of
 the values computed by DLAMCH.  This subroutine is needed because
 DLAMCH does not compensate for poor arithmetic in the upper half of
 the exponent range, as is found on a Cray.
Parameters
[in,out]SMALL
          SMALL is DOUBLE PRECISION
          On entry, the underflow threshold as computed by DLAMCH.
          On exit, if LOG10(LARGE) is sufficiently large, the square
          root of SMALL, otherwise unchanged.
[in,out]LARGE
          LARGE is DOUBLE PRECISION
          On entry, the overflow threshold as computed by DLAMCH.
          On exit, if LOG10(LARGE) is sufficiently large, the square
          root of LARGE, otherwise unchanged.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 76 of file dlabad.f.

76 *
77 * -- LAPACK auxiliary routine (version 3.4.0) --
78 * -- LAPACK is a software package provided by Univ. of Tennessee, --
79 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
80 * November 2011
81 *
82 * .. Scalar Arguments ..
83  DOUBLE PRECISION large, small
84 * ..
85 *
86 * =====================================================================
87 *
88 * .. Intrinsic Functions ..
89  INTRINSIC log10, sqrt
90 * ..
91 * .. Executable Statements ..
92 *
93 * If it looks like we're on a Cray, take the square root of
94 * SMALL and LARGE to avoid overflow and underflow problems.
95 *
96  IF( log10( large ).GT.2000.d0 ) THEN
97  small = sqrt( small )
98  large = sqrt( large )
99  END IF
100 *
101  RETURN
102 *
103 * End of DLABAD
104 *
subroutine dlacpy ( character  UPLO,
integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldb, * )  B,
integer  LDB 
)

DLACPY copies all or part of one two-dimensional array to another.

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

Purpose:
 DLACPY copies all or part of a two-dimensional matrix A to another
 matrix B.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies the part of the matrix A to be copied to B.
          = 'U':      Upper triangular part
          = 'L':      Lower triangular part
          Otherwise:  All of the matrix A
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in]A
          A is DOUBLE PRECISION array, dimension (LDA,N)
          The m by n matrix A.  If UPLO = 'U', only the upper triangle
          or trapezoid is accessed; if UPLO = 'L', only the lower
          triangle or trapezoid is accessed.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]B
          B is DOUBLE PRECISION array, dimension (LDB,N)
          On exit, B = A in the locations specified by UPLO.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,M).
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 105 of file dlacpy.f.

105 *
106 * -- LAPACK auxiliary routine (version 3.4.2) --
107 * -- LAPACK is a software package provided by Univ. of Tennessee, --
108 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
109 * September 2012
110 *
111 * .. Scalar Arguments ..
112  CHARACTER uplo
113  INTEGER lda, ldb, m, n
114 * ..
115 * .. Array Arguments ..
116  DOUBLE PRECISION a( lda, * ), b( ldb, * )
117 * ..
118 *
119 * =====================================================================
120 *
121 * .. Local Scalars ..
122  INTEGER i, j
123 * ..
124 * .. External Functions ..
125  LOGICAL lsame
126  EXTERNAL lsame
127 * ..
128 * .. Intrinsic Functions ..
129  INTRINSIC min
130 * ..
131 * .. Executable Statements ..
132 *
133  IF( lsame( uplo, 'U' ) ) THEN
134  DO 20 j = 1, n
135  DO 10 i = 1, min( j, m )
136  b( i, j ) = a( i, j )
137  10 CONTINUE
138  20 CONTINUE
139  ELSE IF( lsame( uplo, 'L' ) ) THEN
140  DO 40 j = 1, n
141  DO 30 i = j, m
142  b( i, j ) = a( i, j )
143  30 CONTINUE
144  40 CONTINUE
145  ELSE
146  DO 60 j = 1, n
147  DO 50 i = 1, m
148  b( i, j ) = a( i, j )
149  50 CONTINUE
150  60 CONTINUE
151  END IF
152  RETURN
153 *
154 * End of DLACPY
155 *
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dladiv ( double precision  A,
double precision  B,
double precision  C,
double precision  D,
double precision  P,
double precision  Q 
)

DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.

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

Purpose:
 DLADIV performs complex division in  real arithmetic

                       a + i*b
            p + i*q = ---------
                       c + i*d

 The algorithm is due to Michael Baudin and Robert L. Smith
 and can be found in the paper
 "A Robust Complex Division in Scilab"
Parameters
[in]A
          A is DOUBLE PRECISION
[in]B
          B is DOUBLE PRECISION
[in]C
          C is DOUBLE PRECISION
[in]D
          D is DOUBLE PRECISION
          The scalars a, b, c, and d in the above expression.
[out]P
          P is DOUBLE PRECISION
[out]Q
          Q is DOUBLE PRECISION
          The scalars p and q in the above expression.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
January 2013

Definition at line 93 of file dladiv.f.

93 *
94 * -- LAPACK auxiliary routine (version 3.6.0) --
95 * -- LAPACK is a software package provided by Univ. of Tennessee, --
96 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
97 * January 2013
98 *
99 * .. Scalar Arguments ..
100  DOUBLE PRECISION a, b, c, d, p, q
101 * ..
102 *
103 * =====================================================================
104 *
105 * .. Parameters ..
106  DOUBLE PRECISION bs
107  parameter( bs = 2.0d0 )
108  DOUBLE PRECISION half
109  parameter( half = 0.5d0 )
110  DOUBLE PRECISION two
111  parameter( two = 2.0d0 )
112 *
113 * .. Local Scalars ..
114  DOUBLE PRECISION aa, bb, cc, dd, ab, cd, s, ov, un, be, eps
115 * ..
116 * .. External Functions ..
117  DOUBLE PRECISION dlamch
118  EXTERNAL dlamch
119 * ..
120 * .. External Subroutines ..
121  EXTERNAL dladiv1
122 * ..
123 * .. Intrinsic Functions ..
124  INTRINSIC abs, max
125 * ..
126 * .. Executable Statements ..
127 *
128  aa = a
129  bb = b
130  cc = c
131  dd = d
132  ab = max( abs(a), abs(b) )
133  cd = max( abs(c), abs(d) )
134  s = 1.0d0
135 
136  ov = dlamch( 'Overflow threshold' )
137  un = dlamch( 'Safe minimum' )
138  eps = dlamch( 'Epsilon' )
139  be = bs / (eps*eps)
140 
141  IF( ab >= half*ov ) THEN
142  aa = half * aa
143  bb = half * bb
144  s = two * s
145  END IF
146  IF( cd >= half*ov ) THEN
147  cc = half * cc
148  dd = half * dd
149  s = half * s
150  END IF
151  IF( ab <= un*bs/eps ) THEN
152  aa = aa * be
153  bb = bb * be
154  s = s / be
155  END IF
156  IF( cd <= un*bs/eps ) THEN
157  cc = cc * be
158  dd = dd * be
159  s = s * be
160  END IF
161  IF( abs( d ).LE.abs( c ) ) THEN
162  CALL dladiv1(aa, bb, cc, dd, p, q)
163  ELSE
164  CALL dladiv1(bb, aa, dd, cc, p, q)
165  q = -q
166  END IF
167  p = p * s
168  q = q * s
169 *
170  RETURN
171 *
172 * End of DLADIV
173 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dladiv1(A, B, C, D, P, Q)
Definition: dladiv.f:179

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dlae2 ( double precision  A,
double precision  B,
double precision  C,
double precision  RT1,
double precision  RT2 
)

DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.

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

Purpose:
 DLAE2  computes the eigenvalues of a 2-by-2 symmetric matrix
    [  A   B  ]
    [  B   C  ].
 On return, RT1 is the eigenvalue of larger absolute value, and RT2
 is the eigenvalue of smaller absolute value.
Parameters
[in]A
          A is DOUBLE PRECISION
          The (1,1) element of the 2-by-2 matrix.
[in]B
          B is DOUBLE PRECISION
          The (1,2) and (2,1) elements of the 2-by-2 matrix.
[in]C
          C is DOUBLE PRECISION
          The (2,2) element of the 2-by-2 matrix.
[out]RT1
          RT1 is DOUBLE PRECISION
          The eigenvalue of larger absolute value.
[out]RT2
          RT2 is DOUBLE PRECISION
          The eigenvalue of smaller absolute value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
  RT1 is accurate to a few ulps barring over/underflow.

  RT2 may be inaccurate if there is massive cancellation in the
  determinant A*C-B*B; higher precision or correctly rounded or
  correctly truncated arithmetic would be needed to compute RT2
  accurately in all cases.

  Overflow is possible only if RT1 is within a factor of 5 of overflow.
  Underflow is harmless if the input data is 0 or exceeds
     underflow_threshold / macheps.

Definition at line 104 of file dlae2.f.

104 *
105 * -- LAPACK auxiliary routine (version 3.4.2) --
106 * -- LAPACK is a software package provided by Univ. of Tennessee, --
107 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
108 * September 2012
109 *
110 * .. Scalar Arguments ..
111  DOUBLE PRECISION a, b, c, rt1, rt2
112 * ..
113 *
114 * =====================================================================
115 *
116 * .. Parameters ..
117  DOUBLE PRECISION one
118  parameter( one = 1.0d0 )
119  DOUBLE PRECISION two
120  parameter( two = 2.0d0 )
121  DOUBLE PRECISION zero
122  parameter( zero = 0.0d0 )
123  DOUBLE PRECISION half
124  parameter( half = 0.5d0 )
125 * ..
126 * .. Local Scalars ..
127  DOUBLE PRECISION ab, acmn, acmx, adf, df, rt, sm, tb
128 * ..
129 * .. Intrinsic Functions ..
130  INTRINSIC abs, sqrt
131 * ..
132 * .. Executable Statements ..
133 *
134 * Compute the eigenvalues
135 *
136  sm = a + c
137  df = a - c
138  adf = abs( df )
139  tb = b + b
140  ab = abs( tb )
141  IF( abs( a ).GT.abs( c ) ) THEN
142  acmx = a
143  acmn = c
144  ELSE
145  acmx = c
146  acmn = a
147  END IF
148  IF( adf.GT.ab ) THEN
149  rt = adf*sqrt( one+( ab / adf )**2 )
150  ELSE IF( adf.LT.ab ) THEN
151  rt = ab*sqrt( one+( adf / ab )**2 )
152  ELSE
153 *
154 * Includes case AB=ADF=0
155 *
156  rt = ab*sqrt( two )
157  END IF
158  IF( sm.LT.zero ) THEN
159  rt1 = half*( sm-rt )
160 *
161 * Order of execution important.
162 * To get fully accurate smaller eigenvalue,
163 * next line needs to be executed in higher precision.
164 *
165  rt2 = ( acmx / rt1 )*acmn - ( b / rt1 )*b
166  ELSE IF( sm.GT.zero ) THEN
167  rt1 = half*( sm+rt )
168 *
169 * Order of execution important.
170 * To get fully accurate smaller eigenvalue,
171 * next line needs to be executed in higher precision.
172 *
173  rt2 = ( acmx / rt1 )*acmn - ( b / rt1 )*b
174  ELSE
175 *
176 * Includes case RT1 = RT2 = 0
177 *
178  rt1 = half*rt
179  rt2 = -half*rt
180  END IF
181  RETURN
182 *
183 * End of DLAE2
184 *

Here is the caller graph for this function:

subroutine dlaebz ( integer  IJOB,
integer  NITMAX,
integer  N,
integer  MMAX,
integer  MINP,
integer  NBMIN,
double precision  ABSTOL,
double precision  RELTOL,
double precision  PIVMIN,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( * )  E2,
integer, dimension( * )  NVAL,
double precision, dimension( mmax, * )  AB,
double precision, dimension( * )  C,
integer  MOUT,
integer, dimension( mmax, * )  NAB,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

DLAEBZ computes the number of eigenvalues of a real symmetric tridiagonal matrix which are less than or equal to a given value, and performs other tasks required by the routine sstebz.

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

Purpose:
 DLAEBZ contains the iteration loops which compute and use the
 function N(w), which is the count of eigenvalues of a symmetric
 tridiagonal matrix T less than or equal to its argument  w.  It
 performs a choice of two types of loops:

 IJOB=1, followed by
 IJOB=2: It takes as input a list of intervals and returns a list of
         sufficiently small intervals whose union contains the same
         eigenvalues as the union of the original intervals.
         The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
         The output interval (AB(j,1),AB(j,2)] will contain
         eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.

 IJOB=3: It performs a binary search in each input interval
         (AB(j,1),AB(j,2)] for a point  w(j)  such that
         N(w(j))=NVAL(j), and uses  C(j)  as the starting point of
         the search.  If such a w(j) is found, then on output
         AB(j,1)=AB(j,2)=w.  If no such w(j) is found, then on output
         (AB(j,1),AB(j,2)] will be a small interval containing the
         point where N(w) jumps through NVAL(j), unless that point
         lies outside the initial interval.

 Note that the intervals are in all cases half-open intervals,
 i.e., of the form  (a,b] , which includes  b  but not  a .

 To avoid underflow, the matrix should be scaled so that its largest
 element is no greater than  overflow**(1/2) * underflow**(1/4)
 in absolute value.  To assure the most accurate computation
 of small eigenvalues, the matrix should be scaled to be
 not much smaller than that, either.

 See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
 Matrix", Report CS41, Computer Science Dept., Stanford
 University, July 21, 1966

 Note: the arguments are, in general, *not* checked for unreasonable
 values.
Parameters
[in]IJOB
          IJOB is INTEGER
          Specifies what is to be done:
          = 1:  Compute NAB for the initial intervals.
          = 2:  Perform bisection iteration to find eigenvalues of T.
          = 3:  Perform bisection iteration to invert N(w), i.e.,
                to find a point which has a specified number of
                eigenvalues of T to its left.
          Other values will cause DLAEBZ to return with INFO=-1.
[in]NITMAX
          NITMAX is INTEGER
          The maximum number of "levels" of bisection to be
          performed, i.e., an interval of width W will not be made
          smaller than 2^(-NITMAX) * W.  If not all intervals
          have converged after NITMAX iterations, then INFO is set
          to the number of non-converged intervals.
[in]N
          N is INTEGER
          The dimension n of the tridiagonal matrix T.  It must be at
          least 1.
[in]MMAX
          MMAX is INTEGER
          The maximum number of intervals.  If more than MMAX intervals
          are generated, then DLAEBZ will quit with INFO=MMAX+1.
[in]MINP
          MINP is INTEGER
          The initial number of intervals.  It may not be greater than
          MMAX.
[in]NBMIN
          NBMIN is INTEGER
          The smallest number of intervals that should be processed
          using a vector loop.  If zero, then only the scalar loop
          will be used.
[in]ABSTOL
          ABSTOL is DOUBLE PRECISION
          The minimum (absolute) width of an interval.  When an
          interval is narrower than ABSTOL, or than RELTOL times the
          larger (in magnitude) endpoint, then it is considered to be
          sufficiently small, i.e., converged.  This must be at least
          zero.
[in]RELTOL
          RELTOL is DOUBLE PRECISION
          The minimum relative width of an interval.  When an interval
          is narrower than ABSTOL, or than RELTOL times the larger (in
          magnitude) endpoint, then it is considered to be
          sufficiently small, i.e., converged.  Note: this should
          always be at least radix*machine epsilon.
[in]PIVMIN
          PIVMIN is DOUBLE PRECISION
          The minimum absolute value of a "pivot" in the Sturm
          sequence loop.
          This must be at least  max |e(j)**2|*safe_min  and at
          least safe_min, where safe_min is at least
          the smallest number that can divide one without overflow.
[in]D
          D is DOUBLE PRECISION array, dimension (N)
          The diagonal elements of the tridiagonal matrix T.
[in]E
          E is DOUBLE PRECISION array, dimension (N)
          The offdiagonal elements of the tridiagonal matrix T in
          positions 1 through N-1.  E(N) is arbitrary.
[in]E2
          E2 is DOUBLE PRECISION array, dimension (N)
          The squares of the offdiagonal elements of the tridiagonal
          matrix T.  E2(N) is ignored.
[in,out]NVAL
          NVAL is INTEGER array, dimension (MINP)
          If IJOB=1 or 2, not referenced.
          If IJOB=3, the desired values of N(w).  The elements of NVAL
          will be reordered to correspond with the intervals in AB.
          Thus, NVAL(j) on output will not, in general be the same as
          NVAL(j) on input, but it will correspond with the interval
          (AB(j,1),AB(j,2)] on output.
[in,out]AB
          AB is DOUBLE PRECISION array, dimension (MMAX,2)
          The endpoints of the intervals.  AB(j,1) is  a(j), the left
          endpoint of the j-th interval, and AB(j,2) is b(j), the
          right endpoint of the j-th interval.  The input intervals
          will, in general, be modified, split, and reordered by the
          calculation.
[in,out]C
          C is DOUBLE PRECISION array, dimension (MMAX)
          If IJOB=1, ignored.
          If IJOB=2, workspace.
          If IJOB=3, then on input C(j) should be initialized to the
          first search point in the binary search.
[out]MOUT
          MOUT is INTEGER
          If IJOB=1, the number of eigenvalues in the intervals.
          If IJOB=2 or 3, the number of intervals output.
          If IJOB=3, MOUT will equal MINP.
[in,out]NAB
          NAB is INTEGER array, dimension (MMAX,2)
          If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)).
          If IJOB=2, then on input, NAB(i,j) should be set.  It must
             satisfy the condition:
             N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)),
             which means that in interval i only eigenvalues
             NAB(i,1)+1,...,NAB(i,2) will be considered.  Usually,
             NAB(i,j)=N(AB(i,j)), from a previous call to DLAEBZ with
             IJOB=1.
             On output, NAB(i,j) will contain
             max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of
             the input interval that the output interval
             (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the
             the input values of NAB(k,1) and NAB(k,2).
          If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)),
             unless N(w) > NVAL(i) for all search points  w , in which
             case NAB(i,1) will not be modified, i.e., the output
             value will be the same as the input value (modulo
             reorderings -- see NVAL and AB), or unless N(w) < NVAL(i)
             for all search points  w , in which case NAB(i,2) will
             not be modified.  Normally, NAB should be set to some
             distinctive value(s) before DLAEBZ is called.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (MMAX)
          Workspace.
[out]IWORK
          IWORK is INTEGER array, dimension (MMAX)
          Workspace.
[out]INFO
          INFO is INTEGER
          = 0:       All intervals converged.
          = 1--MMAX: The last INFO intervals did not converge.
          = MMAX+1:  More than MMAX intervals were generated.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
      This routine is intended to be called only by other LAPACK
  routines, thus the interface is less user-friendly.  It is intended
  for two purposes:

  (a) finding eigenvalues.  In this case, DLAEBZ should have one or
      more initial intervals set up in AB, and DLAEBZ should be called
      with IJOB=1.  This sets up NAB, and also counts the eigenvalues.
      Intervals with no eigenvalues would usually be thrown out at
      this point.  Also, if not all the eigenvalues in an interval i
      are desired, NAB(i,1) can be increased or NAB(i,2) decreased.
      For example, set NAB(i,1)=NAB(i,2)-1 to get the largest
      eigenvalue.  DLAEBZ is then called with IJOB=2 and MMAX
      no smaller than the value of MOUT returned by the call with
      IJOB=1.  After this (IJOB=2) call, eigenvalues NAB(i,1)+1
      through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the
      tolerance specified by ABSTOL and RELTOL.

  (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l).
      In this case, start with a Gershgorin interval  (a,b).  Set up
      AB to contain 2 search intervals, both initially (a,b).  One
      NVAL element should contain  f-1  and the other should contain  l
      , while C should contain a and b, resp.  NAB(i,1) should be -1
      and NAB(i,2) should be N+1, to flag an error if the desired
      interval does not lie in (a,b).  DLAEBZ is then called with
      IJOB=3.  On exit, if w(f-1) < w(f), then one of the intervals --
      j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while
      if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r
      >= 0, then the interval will have  N(AB(j,1))=NAB(j,1)=f-k and
      N(AB(j,2))=NAB(j,2)=f+r.  The cases w(l) < w(l+1) and
      w(l-r)=...=w(l+k) are handled similarly.

Definition at line 321 of file dlaebz.f.

321 *
322 * -- LAPACK auxiliary routine (version 3.4.2) --
323 * -- LAPACK is a software package provided by Univ. of Tennessee, --
324 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
325 * September 2012
326 *
327 * .. Scalar Arguments ..
328  INTEGER ijob, info, minp, mmax, mout, n, nbmin, nitmax
329  DOUBLE PRECISION abstol, pivmin, reltol
330 * ..
331 * .. Array Arguments ..
332  INTEGER iwork( * ), nab( mmax, * ), nval( * )
333  DOUBLE PRECISION ab( mmax, * ), c( * ), d( * ), e( * ), e2( * ),
334  $ work( * )
335 * ..
336 *
337 * =====================================================================
338 *
339 * .. Parameters ..
340  DOUBLE PRECISION zero, two, half
341  parameter( zero = 0.0d0, two = 2.0d0,
342  $ half = 1.0d0 / two )
343 * ..
344 * .. Local Scalars ..
345  INTEGER itmp1, itmp2, j, ji, jit, jp, kf, kfnew, kl,
346  $ klnew
347  DOUBLE PRECISION tmp1, tmp2
348 * ..
349 * .. Intrinsic Functions ..
350  INTRINSIC abs, max, min
351 * ..
352 * .. Executable Statements ..
353 *
354 * Check for Errors
355 *
356  info = 0
357  IF( ijob.LT.1 .OR. ijob.GT.3 ) THEN
358  info = -1
359  RETURN
360  END IF
361 *
362 * Initialize NAB
363 *
364  IF( ijob.EQ.1 ) THEN
365 *
366 * Compute the number of eigenvalues in the initial intervals.
367 *
368  mout = 0
369  DO 30 ji = 1, minp
370  DO 20 jp = 1, 2
371  tmp1 = d( 1 ) - ab( ji, jp )
372  IF( abs( tmp1 ).LT.pivmin )
373  $ tmp1 = -pivmin
374  nab( ji, jp ) = 0
375  IF( tmp1.LE.zero )
376  $ nab( ji, jp ) = 1
377 *
378  DO 10 j = 2, n
379  tmp1 = d( j ) - e2( j-1 ) / tmp1 - ab( ji, jp )
380  IF( abs( tmp1 ).LT.pivmin )
381  $ tmp1 = -pivmin
382  IF( tmp1.LE.zero )
383  $ nab( ji, jp ) = nab( ji, jp ) + 1
384  10 CONTINUE
385  20 CONTINUE
386  mout = mout + nab( ji, 2 ) - nab( ji, 1 )
387  30 CONTINUE
388  RETURN
389  END IF
390 *
391 * Initialize for loop
392 *
393 * KF and KL have the following meaning:
394 * Intervals 1,...,KF-1 have converged.
395 * Intervals KF,...,KL still need to be refined.
396 *
397  kf = 1
398  kl = minp
399 *
400 * If IJOB=2, initialize C.
401 * If IJOB=3, use the user-supplied starting point.
402 *
403  IF( ijob.EQ.2 ) THEN
404  DO 40 ji = 1, minp
405  c( ji ) = half*( ab( ji, 1 )+ab( ji, 2 ) )
406  40 CONTINUE
407  END IF
408 *
409 * Iteration loop
410 *
411  DO 130 jit = 1, nitmax
412 *
413 * Loop over intervals
414 *
415  IF( kl-kf+1.GE.nbmin .AND. nbmin.GT.0 ) THEN
416 *
417 * Begin of Parallel Version of the loop
418 *
419  DO 60 ji = kf, kl
420 *
421 * Compute N(c), the number of eigenvalues less than c
422 *
423  work( ji ) = d( 1 ) - c( ji )
424  iwork( ji ) = 0
425  IF( work( ji ).LE.pivmin ) THEN
426  iwork( ji ) = 1
427  work( ji ) = min( work( ji ), -pivmin )
428  END IF
429 *
430  DO 50 j = 2, n
431  work( ji ) = d( j ) - e2( j-1 ) / work( ji ) - c( ji )
432  IF( work( ji ).LE.pivmin ) THEN
433  iwork( ji ) = iwork( ji ) + 1
434  work( ji ) = min( work( ji ), -pivmin )
435  END IF
436  50 CONTINUE
437  60 CONTINUE
438 *
439  IF( ijob.LE.2 ) THEN
440 *
441 * IJOB=2: Choose all intervals containing eigenvalues.
442 *
443  klnew = kl
444  DO 70 ji = kf, kl
445 *
446 * Insure that N(w) is monotone
447 *
448  iwork( ji ) = min( nab( ji, 2 ),
449  $ max( nab( ji, 1 ), iwork( ji ) ) )
450 *
451 * Update the Queue -- add intervals if both halves
452 * contain eigenvalues.
453 *
454  IF( iwork( ji ).EQ.nab( ji, 2 ) ) THEN
455 *
456 * No eigenvalue in the upper interval:
457 * just use the lower interval.
458 *
459  ab( ji, 2 ) = c( ji )
460 *
461  ELSE IF( iwork( ji ).EQ.nab( ji, 1 ) ) THEN
462 *
463 * No eigenvalue in the lower interval:
464 * just use the upper interval.
465 *
466  ab( ji, 1 ) = c( ji )
467  ELSE
468  klnew = klnew + 1
469  IF( klnew.LE.mmax ) THEN
470 *
471 * Eigenvalue in both intervals -- add upper to
472 * queue.
473 *
474  ab( klnew, 2 ) = ab( ji, 2 )
475  nab( klnew, 2 ) = nab( ji, 2 )
476  ab( klnew, 1 ) = c( ji )
477  nab( klnew, 1 ) = iwork( ji )
478  ab( ji, 2 ) = c( ji )
479  nab( ji, 2 ) = iwork( ji )
480  ELSE
481  info = mmax + 1
482  END IF
483  END IF
484  70 CONTINUE
485  IF( info.NE.0 )
486  $ RETURN
487  kl = klnew
488  ELSE
489 *
490 * IJOB=3: Binary search. Keep only the interval containing
491 * w s.t. N(w) = NVAL
492 *
493  DO 80 ji = kf, kl
494  IF( iwork( ji ).LE.nval( ji ) ) THEN
495  ab( ji, 1 ) = c( ji )
496  nab( ji, 1 ) = iwork( ji )
497  END IF
498  IF( iwork( ji ).GE.nval( ji ) ) THEN
499  ab( ji, 2 ) = c( ji )
500  nab( ji, 2 ) = iwork( ji )
501  END IF
502  80 CONTINUE
503  END IF
504 *
505  ELSE
506 *
507 * End of Parallel Version of the loop
508 *
509 * Begin of Serial Version of the loop
510 *
511  klnew = kl
512  DO 100 ji = kf, kl
513 *
514 * Compute N(w), the number of eigenvalues less than w
515 *
516  tmp1 = c( ji )
517  tmp2 = d( 1 ) - tmp1
518  itmp1 = 0
519  IF( tmp2.LE.pivmin ) THEN
520  itmp1 = 1
521  tmp2 = min( tmp2, -pivmin )
522  END IF
523 *
524  DO 90 j = 2, n
525  tmp2 = d( j ) - e2( j-1 ) / tmp2 - tmp1
526  IF( tmp2.LE.pivmin ) THEN
527  itmp1 = itmp1 + 1
528  tmp2 = min( tmp2, -pivmin )
529  END IF
530  90 CONTINUE
531 *
532  IF( ijob.LE.2 ) THEN
533 *
534 * IJOB=2: Choose all intervals containing eigenvalues.
535 *
536 * Insure that N(w) is monotone
537 *
538  itmp1 = min( nab( ji, 2 ),
539  $ max( nab( ji, 1 ), itmp1 ) )
540 *
541 * Update the Queue -- add intervals if both halves
542 * contain eigenvalues.
543 *
544  IF( itmp1.EQ.nab( ji, 2 ) ) THEN
545 *
546 * No eigenvalue in the upper interval:
547 * just use the lower interval.
548 *
549  ab( ji, 2 ) = tmp1
550 *
551  ELSE IF( itmp1.EQ.nab( ji, 1 ) ) THEN
552 *
553 * No eigenvalue in the lower interval:
554 * just use the upper interval.
555 *
556  ab( ji, 1 ) = tmp1
557  ELSE IF( klnew.LT.mmax ) THEN
558 *
559 * Eigenvalue in both intervals -- add upper to queue.
560 *
561  klnew = klnew + 1
562  ab( klnew, 2 ) = ab( ji, 2 )
563  nab( klnew, 2 ) = nab( ji, 2 )
564  ab( klnew, 1 ) = tmp1
565  nab( klnew, 1 ) = itmp1
566  ab( ji, 2 ) = tmp1
567  nab( ji, 2 ) = itmp1
568  ELSE
569  info = mmax + 1
570  RETURN
571  END IF
572  ELSE
573 *
574 * IJOB=3: Binary search. Keep only the interval
575 * containing w s.t. N(w) = NVAL
576 *
577  IF( itmp1.LE.nval( ji ) ) THEN
578  ab( ji, 1 ) = tmp1
579  nab( ji, 1 ) = itmp1
580  END IF
581  IF( itmp1.GE.nval( ji ) ) THEN
582  ab( ji, 2 ) = tmp1
583  nab( ji, 2 ) = itmp1
584  END IF
585  END IF
586  100 CONTINUE
587  kl = klnew
588 *
589  END IF
590 *
591 * Check for convergence
592 *
593  kfnew = kf
594  DO 110 ji = kf, kl
595  tmp1 = abs( ab( ji, 2 )-ab( ji, 1 ) )
596  tmp2 = max( abs( ab( ji, 2 ) ), abs( ab( ji, 1 ) ) )
597  IF( tmp1.LT.max( abstol, pivmin, reltol*tmp2 ) .OR.
598  $ nab( ji, 1 ).GE.nab( ji, 2 ) ) THEN
599 *
600 * Converged -- Swap with position KFNEW,
601 * then increment KFNEW
602 *
603  IF( ji.GT.kfnew ) THEN
604  tmp1 = ab( ji, 1 )
605  tmp2 = ab( ji, 2 )
606  itmp1 = nab( ji, 1 )
607  itmp2 = nab( ji, 2 )
608  ab( ji, 1 ) = ab( kfnew, 1 )
609  ab( ji, 2 ) = ab( kfnew, 2 )
610  nab( ji, 1 ) = nab( kfnew, 1 )
611  nab( ji, 2 ) = nab( kfnew, 2 )
612  ab( kfnew, 1 ) = tmp1
613  ab( kfnew, 2 ) = tmp2
614  nab( kfnew, 1 ) = itmp1
615  nab( kfnew, 2 ) = itmp2
616  IF( ijob.EQ.3 ) THEN
617  itmp1 = nval( ji )
618  nval( ji ) = nval( kfnew )
619  nval( kfnew ) = itmp1
620  END IF
621  END IF
622  kfnew = kfnew + 1
623  END IF
624  110 CONTINUE
625  kf = kfnew
626 *
627 * Choose Midpoints
628 *
629  DO 120 ji = kf, kl
630  c( ji ) = half*( ab( ji, 1 )+ab( ji, 2 ) )
631  120 CONTINUE
632 *
633 * If no more intervals to refine, quit.
634 *
635  IF( kf.GT.kl )
636  $ GO TO 140
637  130 CONTINUE
638 *
639 * Converged
640 *
641  140 CONTINUE
642  info = max( kl+1-kf, 0 )
643  mout = kl
644 *
645  RETURN
646 *
647 * End of DLAEBZ
648 *

Here is the caller graph for this function:

subroutine dlaev2 ( double precision  A,
double precision  B,
double precision  C,
double precision  RT1,
double precision  RT2,
double precision  CS1,
double precision  SN1 
)

DLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.

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

Purpose:
 DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
    [  A   B  ]
    [  B   C  ].
 On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
 eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
 eigenvector for RT1, giving the decomposition

    [ CS1  SN1 ] [  A   B  ] [ CS1 -SN1 ]  =  [ RT1  0  ]
    [-SN1  CS1 ] [  B   C  ] [ SN1  CS1 ]     [  0  RT2 ].
Parameters
[in]A
          A is DOUBLE PRECISION
          The (1,1) element of the 2-by-2 matrix.
[in]B
          B is DOUBLE PRECISION
          The (1,2) element and the conjugate of the (2,1) element of
          the 2-by-2 matrix.
[in]C
          C is DOUBLE PRECISION
          The (2,2) element of the 2-by-2 matrix.
[out]RT1
          RT1 is DOUBLE PRECISION
          The eigenvalue of larger absolute value.
[out]RT2
          RT2 is DOUBLE PRECISION
          The eigenvalue of smaller absolute value.
[out]CS1
          CS1 is DOUBLE PRECISION
[out]SN1
          SN1 is DOUBLE PRECISION
          The vector (CS1, SN1) is a unit right eigenvector for RT1.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
  RT1 is accurate to a few ulps barring over/underflow.

  RT2 may be inaccurate if there is massive cancellation in the
  determinant A*C-B*B; higher precision or correctly rounded or
  correctly truncated arithmetic would be needed to compute RT2
  accurately in all cases.

  CS1 and SN1 are accurate to a few ulps barring over/underflow.

  Overflow is possible only if RT1 is within a factor of 5 of overflow.
  Underflow is harmless if the input data is 0 or exceeds
     underflow_threshold / macheps.

Definition at line 122 of file dlaev2.f.

122 *
123 * -- LAPACK auxiliary routine (version 3.4.2) --
124 * -- LAPACK is a software package provided by Univ. of Tennessee, --
125 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
126 * September 2012
127 *
128 * .. Scalar Arguments ..
129  DOUBLE PRECISION a, b, c, cs1, rt1, rt2, sn1
130 * ..
131 *
132 * =====================================================================
133 *
134 * .. Parameters ..
135  DOUBLE PRECISION one
136  parameter( one = 1.0d0 )
137  DOUBLE PRECISION two
138  parameter( two = 2.0d0 )
139  DOUBLE PRECISION zero
140  parameter( zero = 0.0d0 )
141  DOUBLE PRECISION half
142  parameter( half = 0.5d0 )
143 * ..
144 * .. Local Scalars ..
145  INTEGER sgn1, sgn2
146  DOUBLE PRECISION ab, acmn, acmx, acs, adf, cs, ct, df, rt, sm,
147  $ tb, tn
148 * ..
149 * .. Intrinsic Functions ..
150  INTRINSIC abs, sqrt
151 * ..
152 * .. Executable Statements ..
153 *
154 * Compute the eigenvalues
155 *
156  sm = a + c
157  df = a - c
158  adf = abs( df )
159  tb = b + b
160  ab = abs( tb )
161  IF( abs( a ).GT.abs( c ) ) THEN
162  acmx = a
163  acmn = c
164  ELSE
165  acmx = c
166  acmn = a
167  END IF
168  IF( adf.GT.ab ) THEN
169  rt = adf*sqrt( one+( ab / adf )**2 )
170  ELSE IF( adf.LT.ab ) THEN
171  rt = ab*sqrt( one+( adf / ab )**2 )
172  ELSE
173 *
174 * Includes case AB=ADF=0
175 *
176  rt = ab*sqrt( two )
177  END IF
178  IF( sm.LT.zero ) THEN
179  rt1 = half*( sm-rt )
180  sgn1 = -1
181 *
182 * Order of execution important.
183 * To get fully accurate smaller eigenvalue,
184 * next line needs to be executed in higher precision.
185 *
186  rt2 = ( acmx / rt1 )*acmn - ( b / rt1 )*b
187  ELSE IF( sm.GT.zero ) THEN
188  rt1 = half*( sm+rt )
189  sgn1 = 1
190 *
191 * Order of execution important.
192 * To get fully accurate smaller eigenvalue,
193 * next line needs to be executed in higher precision.
194 *
195  rt2 = ( acmx / rt1 )*acmn - ( b / rt1 )*b
196  ELSE
197 *
198 * Includes case RT1 = RT2 = 0
199 *
200  rt1 = half*rt
201  rt2 = -half*rt
202  sgn1 = 1
203  END IF
204 *
205 * Compute the eigenvector
206 *
207  IF( df.GE.zero ) THEN
208  cs = df + rt
209  sgn2 = 1
210  ELSE
211  cs = df - rt
212  sgn2 = -1
213  END IF
214  acs = abs( cs )
215  IF( acs.GT.ab ) THEN
216  ct = -tb / cs
217  sn1 = one / sqrt( one+ct*ct )
218  cs1 = ct*sn1
219  ELSE
220  IF( ab.EQ.zero ) THEN
221  cs1 = one
222  sn1 = zero
223  ELSE
224  tn = -cs / tb
225  cs1 = one / sqrt( one+tn*tn )
226  sn1 = tn*cs1
227  END IF
228  END IF
229  IF( sgn1.EQ.sgn2 ) THEN
230  tn = cs1
231  cs1 = -sn1
232  sn1 = tn
233  END IF
234  RETURN
235 *
236 * End of DLAEV2
237 *

Here is the caller graph for this function:

subroutine dlagts ( integer  JOB,
integer  N,
double precision, dimension( * )  A,
double precision, dimension( * )  B,
double precision, dimension( * )  C,
double precision, dimension( * )  D,
integer, dimension( * )  IN,
double precision, dimension( * )  Y,
double precision  TOL,
integer  INFO 
)

DLAGTS solves the system of equations (T-λI)x = y or (T-λI)Tx = y,where T is a general tridiagonal matrix and λ a scalar, using the LU factorization computed by slagtf.

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

Purpose:
 DLAGTS may be used to solve one of the systems of equations

    (T - lambda*I)*x = y   or   (T - lambda*I)**T*x = y,

 where T is an n by n tridiagonal matrix, for x, following the
 factorization of (T - lambda*I) as

    (T - lambda*I) = P*L*U ,

 by routine DLAGTF. The choice of equation to be solved is
 controlled by the argument JOB, and in each case there is an option
 to perturb zero or very small diagonal elements of U, this option
 being intended for use in applications such as inverse iteration.
Parameters
[in]JOB
          JOB is INTEGER
          Specifies the job to be performed by DLAGTS as follows:
          =  1: The equations  (T - lambda*I)x = y  are to be solved,
                but diagonal elements of U are not to be perturbed.
          = -1: The equations  (T - lambda*I)x = y  are to be solved
                and, if overflow would otherwise occur, the diagonal
                elements of U are to be perturbed. See argument TOL
                below.
          =  2: The equations  (T - lambda*I)**Tx = y  are to be solved,
                but diagonal elements of U are not to be perturbed.
          = -2: The equations  (T - lambda*I)**Tx = y  are to be solved
                and, if overflow would otherwise occur, the diagonal
                elements of U are to be perturbed. See argument TOL
                below.
[in]N
          N is INTEGER
          The order of the matrix T.
[in]A
          A is DOUBLE PRECISION array, dimension (N)
          On entry, A must contain the diagonal elements of U as
          returned from DLAGTF.
[in]B
          B is DOUBLE PRECISION array, dimension (N-1)
          On entry, B must contain the first super-diagonal elements of
          U as returned from DLAGTF.
[in]C
          C is DOUBLE PRECISION array, dimension (N-1)
          On entry, C must contain the sub-diagonal elements of L as
          returned from DLAGTF.
[in]D
          D is DOUBLE PRECISION array, dimension (N-2)
          On entry, D must contain the second super-diagonal elements
          of U as returned from DLAGTF.
[in]IN
          IN is INTEGER array, dimension (N)
          On entry, IN must contain details of the matrix P as returned
          from DLAGTF.
[in,out]Y
          Y is DOUBLE PRECISION array, dimension (N)
          On entry, the right hand side vector y.
          On exit, Y is overwritten by the solution vector x.
[in,out]TOL
          TOL is DOUBLE PRECISION
          On entry, with  JOB .lt. 0, TOL should be the minimum
          perturbation to be made to very small diagonal elements of U.
          TOL should normally be chosen as about eps*norm(U), where eps
          is the relative machine precision, but if TOL is supplied as
          non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
          If  JOB .gt. 0  then TOL is not referenced.

          On exit, TOL is changed as described above, only if TOL is
          non-positive on entry. Otherwise TOL is unchanged.
[out]INFO
          INFO is INTEGER
          = 0   : successful exit
          .lt. 0: if INFO = -i, the i-th argument had an illegal value
          .gt. 0: overflow would occur when computing the INFO(th)
                  element of the solution vector x. This can only occur
                  when JOB is supplied as positive and either means
                  that a diagonal element of U is very small, or that
                  the elements of the right-hand side vector y are very
                  large.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 163 of file dlagts.f.

163 *
164 * -- LAPACK auxiliary routine (version 3.4.2) --
165 * -- LAPACK is a software package provided by Univ. of Tennessee, --
166 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
167 * September 2012
168 *
169 * .. Scalar Arguments ..
170  INTEGER info, job, n
171  DOUBLE PRECISION tol
172 * ..
173 * .. Array Arguments ..
174  INTEGER in( * )
175  DOUBLE PRECISION a( * ), b( * ), c( * ), d( * ), y( * )
176 * ..
177 *
178 * =====================================================================
179 *
180 * .. Parameters ..
181  DOUBLE PRECISION one, zero
182  parameter( one = 1.0d+0, zero = 0.0d+0 )
183 * ..
184 * .. Local Scalars ..
185  INTEGER k
186  DOUBLE PRECISION absak, ak, bignum, eps, pert, sfmin, temp
187 * ..
188 * .. Intrinsic Functions ..
189  INTRINSIC abs, max, sign
190 * ..
191 * .. External Functions ..
192  DOUBLE PRECISION dlamch
193  EXTERNAL dlamch
194 * ..
195 * .. External Subroutines ..
196  EXTERNAL xerbla
197 * ..
198 * .. Executable Statements ..
199 *
200  info = 0
201  IF( ( abs( job ).GT.2 ) .OR. ( job.EQ.0 ) ) THEN
202  info = -1
203  ELSE IF( n.LT.0 ) THEN
204  info = -2
205  END IF
206  IF( info.NE.0 ) THEN
207  CALL xerbla( 'DLAGTS', -info )
208  RETURN
209  END IF
210 *
211  IF( n.EQ.0 )
212  $ RETURN
213 *
214  eps = dlamch( 'Epsilon' )
215  sfmin = dlamch( 'Safe minimum' )
216  bignum = one / sfmin
217 *
218  IF( job.LT.0 ) THEN
219  IF( tol.LE.zero ) THEN
220  tol = abs( a( 1 ) )
221  IF( n.GT.1 )
222  $ tol = max( tol, abs( a( 2 ) ), abs( b( 1 ) ) )
223  DO 10 k = 3, n
224  tol = max( tol, abs( a( k ) ), abs( b( k-1 ) ),
225  $ abs( d( k-2 ) ) )
226  10 CONTINUE
227  tol = tol*eps
228  IF( tol.EQ.zero )
229  $ tol = eps
230  END IF
231  END IF
232 *
233  IF( abs( job ).EQ.1 ) THEN
234  DO 20 k = 2, n
235  IF( in( k-1 ).EQ.0 ) THEN
236  y( k ) = y( k ) - c( k-1 )*y( k-1 )
237  ELSE
238  temp = y( k-1 )
239  y( k-1 ) = y( k )
240  y( k ) = temp - c( k-1 )*y( k )
241  END IF
242  20 CONTINUE
243  IF( job.EQ.1 ) THEN
244  DO 30 k = n, 1, -1
245  IF( k.LE.n-2 ) THEN
246  temp = y( k ) - b( k )*y( k+1 ) - d( k )*y( k+2 )
247  ELSE IF( k.EQ.n-1 ) THEN
248  temp = y( k ) - b( k )*y( k+1 )
249  ELSE
250  temp = y( k )
251  END IF
252  ak = a( k )
253  absak = abs( ak )
254  IF( absak.LT.one ) THEN
255  IF( absak.LT.sfmin ) THEN
256  IF( absak.EQ.zero .OR. abs( temp )*sfmin.GT.absak )
257  $ THEN
258  info = k
259  RETURN
260  ELSE
261  temp = temp*bignum
262  ak = ak*bignum
263  END IF
264  ELSE IF( abs( temp ).GT.absak*bignum ) THEN
265  info = k
266  RETURN
267  END IF
268  END IF
269  y( k ) = temp / ak
270  30 CONTINUE
271  ELSE
272  DO 50 k = n, 1, -1
273  IF( k.LE.n-2 ) THEN
274  temp = y( k ) - b( k )*y( k+1 ) - d( k )*y( k+2 )
275  ELSE IF( k.EQ.n-1 ) THEN
276  temp = y( k ) - b( k )*y( k+1 )
277  ELSE
278  temp = y( k )
279  END IF
280  ak = a( k )
281  pert = sign( tol, ak )
282  40 CONTINUE
283  absak = abs( ak )
284  IF( absak.LT.one ) THEN
285  IF( absak.LT.sfmin ) THEN
286  IF( absak.EQ.zero .OR. abs( temp )*sfmin.GT.absak )
287  $ THEN
288  ak = ak + pert
289  pert = 2*pert
290  GO TO 40
291  ELSE
292  temp = temp*bignum
293  ak = ak*bignum
294  END IF
295  ELSE IF( abs( temp ).GT.absak*bignum ) THEN
296  ak = ak + pert
297  pert = 2*pert
298  GO TO 40
299  END IF
300  END IF
301  y( k ) = temp / ak
302  50 CONTINUE
303  END IF
304  ELSE
305 *
306 * Come to here if JOB = 2 or -2
307 *
308  IF( job.EQ.2 ) THEN
309  DO 60 k = 1, n
310  IF( k.GE.3 ) THEN
311  temp = y( k ) - b( k-1 )*y( k-1 ) - d( k-2 )*y( k-2 )
312  ELSE IF( k.EQ.2 ) THEN
313  temp = y( k ) - b( k-1 )*y( k-1 )
314  ELSE
315  temp = y( k )
316  END IF
317  ak = a( k )
318  absak = abs( ak )
319  IF( absak.LT.one ) THEN
320  IF( absak.LT.sfmin ) THEN
321  IF( absak.EQ.zero .OR. abs( temp )*sfmin.GT.absak )
322  $ THEN
323  info = k
324  RETURN
325  ELSE
326  temp = temp*bignum
327  ak = ak*bignum
328  END IF
329  ELSE IF( abs( temp ).GT.absak*bignum ) THEN
330  info = k
331  RETURN
332  END IF
333  END IF
334  y( k ) = temp / ak
335  60 CONTINUE
336  ELSE
337  DO 80 k = 1, n
338  IF( k.GE.3 ) THEN
339  temp = y( k ) - b( k-1 )*y( k-1 ) - d( k-2 )*y( k-2 )
340  ELSE IF( k.EQ.2 ) THEN
341  temp = y( k ) - b( k-1 )*y( k-1 )
342  ELSE
343  temp = y( k )
344  END IF
345  ak = a( k )
346  pert = sign( tol, ak )
347  70 CONTINUE
348  absak = abs( ak )
349  IF( absak.LT.one ) THEN
350  IF( absak.LT.sfmin ) THEN
351  IF( absak.EQ.zero .OR. abs( temp )*sfmin.GT.absak )
352  $ THEN
353  ak = ak + pert
354  pert = 2*pert
355  GO TO 70
356  ELSE
357  temp = temp*bignum
358  ak = ak*bignum
359  END IF
360  ELSE IF( abs( temp ).GT.absak*bignum ) THEN
361  ak = ak + pert
362  pert = 2*pert
363  GO TO 70
364  END IF
365  END IF
366  y( k ) = temp / ak
367  80 CONTINUE
368  END IF
369 *
370  DO 90 k = n, 2, -1
371  IF( in( k-1 ).EQ.0 ) THEN
372  y( k-1 ) = y( k-1 ) - c( k-1 )*y( k )
373  ELSE
374  temp = y( k-1 )
375  y( k-1 ) = y( k )
376  y( k ) = temp - c( k-1 )*y( k )
377  END IF
378  90 CONTINUE
379  END IF
380 *
381 * End of DLAGTS
382 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65

Here is the call graph for this function:

Here is the caller graph for this function:

logical function dlaisnan ( double precision  DIN1,
double precision  DIN2 
)

DLAISNAN tests input for NaN by comparing two arguments for inequality.

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

Purpose:
 This routine is not for general use.  It exists solely to avoid
 over-optimization in DISNAN.

 DLAISNAN checks for NaNs by comparing its two arguments for
 inequality.  NaN is the only floating-point value where NaN != NaN
 returns .TRUE.  To check for NaNs, pass the same variable as both
 arguments.

 A compiler must assume that the two arguments are
 not the same variable, and the test will not be optimized away.
 Interprocedural or whole-program optimization may delete this
 test.  The ISNAN functions will be replaced by the correct
 Fortran 03 intrinsic once the intrinsic is widely available.
Parameters
[in]DIN1
          DIN1 is DOUBLE PRECISION
[in]DIN2
          DIN2 is DOUBLE PRECISION
          Two numbers to compare for inequality.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 76 of file dlaisnan.f.

76 *
77 * -- LAPACK auxiliary routine (version 3.4.2) --
78 * -- LAPACK is a software package provided by Univ. of Tennessee, --
79 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
80 * September 2012
81 *
82 * .. Scalar Arguments ..
83  DOUBLE PRECISION din1, din2
84 * ..
85 *
86 * =====================================================================
87 *
88 * .. Executable Statements ..
89  dlaisnan = (din1.NE.din2)
90  RETURN
logical function dlaisnan(DIN1, DIN2)
DLAISNAN tests input for NaN by comparing two arguments for inequality.
Definition: dlaisnan.f:76
subroutine dlamc1 ( integer  BETA,
integer  T,
logical  RND,
logical  IEEE1 
)

DLAMC1

Purpose:

 DLAMC1 determines the machine parameters given by BETA, T, RND, and
 IEEE1.
Parameters
[out]BETA
          The base of the machine.
[out]T
          The number of ( BETA ) digits in the mantissa.
[out]RND
          Specifies whether proper rounding  ( RND = .TRUE. )  or
          chopping  ( RND = .FALSE. )  occurs in addition. This may not
          be a reliable guide to the way in which the machine performs
          its arithmetic.
[out]IEEE1
          Specifies whether rounding appears to be done in the IEEE
          'round to nearest' style.
Author
LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
Date
April 2012

Further Details

  The routine is based on the routine  ENVRON  by Malcolm and
  incorporates suggestions by Gentleman and Marovich. See

     Malcolm M. A. (1972) Algorithms to reveal properties of
        floating-point arithmetic. Comms. of the ACM, 15, 949-951.

     Gentleman W. M. and Marovich S. B. (1974) More on algorithms
        that reveal properties of floating point arithmetic units.
        Comms. of the ACM, 17, 276-277.

Definition at line 207 of file dlamchf77.f.

207 *
208 * -- LAPACK auxiliary routine (version 3.4.1) --
209 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
210 * November 2010
211 *
212 * .. Scalar Arguments ..
213  LOGICAL ieee1, rnd
214  INTEGER beta, t
215 * ..
216 * =====================================================================
217 *
218 * .. Local Scalars ..
219  LOGICAL first, lieee1, lrnd
220  INTEGER lbeta, lt
221  DOUBLE PRECISION a, b, c, f, one, qtr, savec, t1, t2
222 * ..
223 * .. External Functions ..
224  DOUBLE PRECISION dlamc3
225  EXTERNAL dlamc3
226 * ..
227 * .. Save statement ..
228  SAVE first, lieee1, lbeta, lrnd, lt
229 * ..
230 * .. Data statements ..
231  DATA first / .true. /
232 * ..
233 * .. Executable Statements ..
234 *
235  IF( first ) THEN
236  one = 1
237 *
238 * LBETA, LIEEE1, LT and LRND are the local values of BETA,
239 * IEEE1, T and RND.
240 *
241 * Throughout this routine we use the function DLAMC3 to ensure
242 * that relevant values are stored and not held in registers, or
243 * are not affected by optimizers.
244 *
245 * Compute a = 2.0**m with the smallest positive integer m such
246 * that
247 *
248 * fl( a + 1.0 ) = a.
249 *
250  a = 1
251  c = 1
252 *
253 *+ WHILE( C.EQ.ONE )LOOP
254  10 CONTINUE
255  IF( c.EQ.one ) THEN
256  a = 2*a
257  c = dlamc3( a, one )
258  c = dlamc3( c, -a )
259  GO TO 10
260  END IF
261 *+ END WHILE
262 *
263 * Now compute b = 2.0**m with the smallest positive integer m
264 * such that
265 *
266 * fl( a + b ) .gt. a.
267 *
268  b = 1
269  c = dlamc3( a, b )
270 *
271 *+ WHILE( C.EQ.A )LOOP
272  20 CONTINUE
273  IF( c.EQ.a ) THEN
274  b = 2*b
275  c = dlamc3( a, b )
276  GO TO 20
277  END IF
278 *+ END WHILE
279 *
280 * Now compute the base. a and c are neighbouring floating point
281 * numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
282 * their difference is beta. Adding 0.25 to c is to ensure that it
283 * is truncated to beta and not ( beta - 1 ).
284 *
285  qtr = one / 4
286  savec = c
287  c = dlamc3( c, -a )
288  lbeta = c + qtr
289 *
290 * Now determine whether rounding or chopping occurs, by adding a
291 * bit less than beta/2 and a bit more than beta/2 to a.
292 *
293  b = lbeta
294  f = dlamc3( b / 2, -b / 100 )
295  c = dlamc3( f, a )
296  IF( c.EQ.a ) THEN
297  lrnd = .true.
298  ELSE
299  lrnd = .false.
300  END IF
301  f = dlamc3( b / 2, b / 100 )
302  c = dlamc3( f, a )
303  IF( ( lrnd ) .AND. ( c.EQ.a ) )
304  $ lrnd = .false.
305 *
306 * Try and decide whether rounding is done in the IEEE 'round to
307 * nearest' style. B/2 is half a unit in the last place of the two
308 * numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
309 * zero, and SAVEC is odd. Thus adding B/2 to A should not change
310 * A, but adding B/2 to SAVEC should change SAVEC.
311 *
312  t1 = dlamc3( b / 2, a )
313  t2 = dlamc3( b / 2, savec )
314  lieee1 = ( t1.EQ.a ) .AND. ( t2.GT.savec ) .AND. lrnd
315 *
316 * Now find the mantissa, t. It should be the integer part of
317 * log to the base beta of a, however it is safer to determine t
318 * by powering. So we find t as the smallest positive integer for
319 * which
320 *
321 * fl( beta**t + 1.0 ) = 1.0.
322 *
323  lt = 0
324  a = 1
325  c = 1
326 *
327 *+ WHILE( C.EQ.ONE )LOOP
328  30 CONTINUE
329  IF( c.EQ.one ) THEN
330  lt = lt + 1
331  a = a*lbeta
332  c = dlamc3( a, one )
333  c = dlamc3( c, -a )
334  GO TO 30
335  END IF
336 *+ END WHILE
337 *
338  END IF
339 *
340  beta = lbeta
341  t = lt
342  rnd = lrnd
343  ieee1 = lieee1
344  first = .false.
345  RETURN
346 *
347 * End of DLAMC1
348 *
double precision function dlamc3(A, B)
DLAMC3
Definition: dlamch.f:169

Here is the caller graph for this function:

subroutine dlamc2 ( integer  BETA,
integer  T,
logical  RND,
double precision  EPS,
integer  EMIN,
double precision  RMIN,
integer  EMAX,
double precision  RMAX 
)

DLAMC2

Purpose:

 DLAMC2 determines the machine parameters specified in its argument
 list.
Author
LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
Date
April 2012
Parameters
[out]BETA
          The base of the machine.
[out]T
          The number of ( BETA ) digits in the mantissa.
[out]RND
          Specifies whether proper rounding  ( RND = .TRUE. )  or
          chopping  ( RND = .FALSE. )  occurs in addition. This may not
          be a reliable guide to the way in which the machine performs
          its arithmetic.
[out]EPS
          The smallest positive number such that
             fl( 1.0 - EPS ) .LT. 1.0,
          where fl denotes the computed value.
[out]EMIN
          The minimum exponent before (gradual) underflow occurs.
[out]RMIN
          The smallest normalized number for the machine, given by
          BASE**( EMIN - 1 ), where  BASE  is the floating point value
          of BETA.
[out]EMAX
          The maximum exponent before overflow occurs.
[out]RMAX
          The largest positive number for the machine, given by
          BASE**EMAX * ( 1 - EPS ), where  BASE  is the floating point
          value of BETA.

Further Details

  The computation of  EPS  is based on a routine PARANOIA by
  W. Kahan of the University of California at Berkeley.

Definition at line 420 of file dlamchf77.f.

420 *
421 * -- LAPACK auxiliary routine (version 3.4.1) --
422 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
423 * November 2010
424 *
425 * .. Scalar Arguments ..
426  LOGICAL rnd
427  INTEGER beta, emax, emin, t
428  DOUBLE PRECISION eps, rmax, rmin
429 * ..
430 * =====================================================================
431 *
432 * .. Local Scalars ..
433  LOGICAL first, ieee, iwarn, lieee1, lrnd
434  INTEGER gnmin, gpmin, i, lbeta, lemax, lemin, lt,
435  $ ngnmin, ngpmin
436  DOUBLE PRECISION a, b, c, half, leps, lrmax, lrmin, one, rbase,
437  $ sixth, small, third, two, zero
438 * ..
439 * .. External Functions ..
440  DOUBLE PRECISION dlamc3
441  EXTERNAL dlamc3
442 * ..
443 * .. External Subroutines ..
444  EXTERNAL dlamc1, dlamc4, dlamc5
445 * ..
446 * .. Intrinsic Functions ..
447  INTRINSIC abs, max, min
448 * ..
449 * .. Save statement ..
450  SAVE first, iwarn, lbeta, lemax, lemin, leps, lrmax,
451  $ lrmin, lt
452 * ..
453 * .. Data statements ..
454  DATA first / .true. / , iwarn / .false. /
455 * ..
456 * .. Executable Statements ..
457 *
458  IF( first ) THEN
459  zero = 0
460  one = 1
461  two = 2
462 *
463 * LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
464 * BETA, T, RND, EPS, EMIN and RMIN.
465 *
466 * Throughout this routine we use the function DLAMC3 to ensure
467 * that relevant values are stored and not held in registers, or
468 * are not affected by optimizers.
469 *
470 * DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
471 *
472  CALL dlamc1( lbeta, lt, lrnd, lieee1 )
473 *
474 * Start to find EPS.
475 *
476  b = lbeta
477  a = b**( -lt )
478  leps = a
479 *
480 * Try some tricks to see whether or not this is the correct EPS.
481 *
482  b = two / 3
483  half = one / 2
484  sixth = dlamc3( b, -half )
485  third = dlamc3( sixth, sixth )
486  b = dlamc3( third, -half )
487  b = dlamc3( b, sixth )
488  b = abs( b )
489  IF( b.LT.leps )
490  $ b = leps
491 *
492  leps = 1
493 *
494 *+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
495  10 CONTINUE
496  IF( ( leps.GT.b ) .AND. ( b.GT.zero ) ) THEN
497  leps = b
498  c = dlamc3( half*leps, ( two**5 )*( leps**2 ) )
499  c = dlamc3( half, -c )
500  b = dlamc3( half, c )
501  c = dlamc3( half, -b )
502  b = dlamc3( half, c )
503  GO TO 10
504  END IF
505 *+ END WHILE
506 *
507  IF( a.LT.leps )
508  $ leps = a
509 *
510 * Computation of EPS complete.
511 *
512 * Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
513 * Keep dividing A by BETA until (gradual) underflow occurs. This
514 * is detected when we cannot recover the previous A.
515 *
516  rbase = one / lbeta
517  small = one
518  DO 20 i = 1, 3
519  small = dlamc3( small*rbase, zero )
520  20 CONTINUE
521  a = dlamc3( one, small )
522  CALL dlamc4( ngpmin, one, lbeta )
523  CALL dlamc4( ngnmin, -one, lbeta )
524  CALL dlamc4( gpmin, a, lbeta )
525  CALL dlamc4( gnmin, -a, lbeta )
526  ieee = .false.
527 *
528  IF( ( ngpmin.EQ.ngnmin ) .AND. ( gpmin.EQ.gnmin ) ) THEN
529  IF( ngpmin.EQ.gpmin ) THEN
530  lemin = ngpmin
531 * ( Non twos-complement machines, no gradual underflow;
532 * e.g., VAX )
533  ELSE IF( ( gpmin-ngpmin ).EQ.3 ) THEN
534  lemin = ngpmin - 1 + lt
535  ieee = .true.
536 * ( Non twos-complement machines, with gradual underflow;
537 * e.g., IEEE standard followers )
538  ELSE
539  lemin = min( ngpmin, gpmin )
540 * ( A guess; no known machine )
541  iwarn = .true.
542  END IF
543 *
544  ELSE IF( ( ngpmin.EQ.gpmin ) .AND. ( ngnmin.EQ.gnmin ) ) THEN
545  IF( abs( ngpmin-ngnmin ).EQ.1 ) THEN
546  lemin = max( ngpmin, ngnmin )
547 * ( Twos-complement machines, no gradual underflow;
548 * e.g., CYBER 205 )
549  ELSE
550  lemin = min( ngpmin, ngnmin )
551 * ( A guess; no known machine )
552  iwarn = .true.
553  END IF
554 *
555  ELSE IF( ( abs( ngpmin-ngnmin ).EQ.1 ) .AND.
556  $ ( gpmin.EQ.gnmin ) ) THEN
557  IF( ( gpmin-min( ngpmin, ngnmin ) ).EQ.3 ) THEN
558  lemin = max( ngpmin, ngnmin ) - 1 + lt
559 * ( Twos-complement machines with gradual underflow;
560 * no known machine )
561  ELSE
562  lemin = min( ngpmin, ngnmin )
563 * ( A guess; no known machine )
564  iwarn = .true.
565  END IF
566 *
567  ELSE
568  lemin = min( ngpmin, ngnmin, gpmin, gnmin )
569 * ( A guess; no known machine )
570  iwarn = .true.
571  END IF
572  first = .false.
573 ***
574 * Comment out this if block if EMIN is ok
575  IF( iwarn ) THEN
576  first = .true.
577  WRITE( 6, fmt = 9999 )lemin
578  END IF
579 ***
580 *
581 * Assume IEEE arithmetic if we found denormalised numbers above,
582 * or if arithmetic seems to round in the IEEE style, determined
583 * in routine DLAMC1. A true IEEE machine should have both things
584 * true; however, faulty machines may have one or the other.
585 *
586  ieee = ieee .OR. lieee1
587 *
588 * Compute RMIN by successive division by BETA. We could compute
589 * RMIN as BASE**( EMIN - 1 ), but some machines underflow during
590 * this computation.
591 *
592  lrmin = 1
593  DO 30 i = 1, 1 - lemin
594  lrmin = dlamc3( lrmin*rbase, zero )
595  30 CONTINUE
596 *
597 * Finally, call DLAMC5 to compute EMAX and RMAX.
598 *
599  CALL dlamc5( lbeta, lt, lemin, ieee, lemax, lrmax )
600  END IF
601 *
602  beta = lbeta
603  t = lt
604  rnd = lrnd
605  eps = leps
606  emin = lemin
607  rmin = lrmin
608  emax = lemax
609  rmax = lrmax
610 *
611  RETURN
612 *
613  9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
614  $ ' EMIN = ', i8, /
615  $ ' If, after inspection, the value EMIN looks',
616  $ ' acceptable please comment out ',
617  $ / ' the IF block as marked within the code of routine',
618  $ ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / )
619 *
620 * End of DLAMC2
621 *
subroutine dlamc5(BETA, P, EMIN, IEEE, EMAX, RMAX)
DLAMC5
Definition: dlamchf77.f:797
subroutine dlamc4(EMIN, START, BASE)
DLAMC4
Definition: dlamchf77.f:690
double precision function dlamc3(A, B)
DLAMC3
Definition: dlamch.f:169
subroutine dlamc1(BETA, T, RND, IEEE1)
DLAMC1
Definition: dlamchf77.f:207

Here is the call graph for this function:

double precision function dlamc3 ( double precision  A,
double precision  B 
)

DLAMC3

Purpose:

 DLAMC3  is intended to force  A  and  B  to be stored prior to doing
 the addition of  A  and  B ,  for use in situations where optimizers
 might hold one of these in a register.
Author
LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
Date
November 2015
Parameters
[in]A
          A is a DOUBLE PRECISION
[in]B
          B is a DOUBLE PRECISION
          The values A and B.

Definition at line 169 of file dlamch.f.

169 *
170 * -- LAPACK auxiliary routine (version 3.6.0) --
171 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
172 * November 2010
173 *
174 * .. Scalar Arguments ..
175  DOUBLE PRECISION a, b
176 * ..
177 * =====================================================================
178 *
179 * .. Executable Statements ..
180 *
181  dlamc3 = a + b
182 *
183  RETURN
184 *
185 * End of DLAMC3
186 *
double precision function dlamc3(A, B)
DLAMC3
Definition: dlamch.f:169

Here is the caller graph for this function:

double precision function dlamch ( character  CMACH)

DLAMCH

DLAMCHF77 deprecated

Purpose:
 DLAMCH determines double precision machine parameters.
Parameters
[in]CMACH
          Specifies the value to be returned by DLAMCH:
          = 'E' or 'e',   DLAMCH := eps
          = 'S' or 's ,   DLAMCH := sfmin
          = 'B' or 'b',   DLAMCH := base
          = 'P' or 'p',   DLAMCH := eps*base
          = 'N' or 'n',   DLAMCH := t
          = 'R' or 'r',   DLAMCH := rnd
          = 'M' or 'm',   DLAMCH := emin
          = 'U' or 'u',   DLAMCH := rmin
          = 'L' or 'l',   DLAMCH := emax
          = 'O' or 'o',   DLAMCH := rmax
          where
          eps   = relative machine precision
          sfmin = safe minimum, such that 1/sfmin does not overflow
          base  = base of the machine
          prec  = eps*base
          t     = number of (base) digits in the mantissa
          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise
          emin  = minimum exponent before (gradual) underflow
          rmin  = underflow threshold - base**(emin-1)
          emax  = largest exponent before overflow
          rmax  = overflow threshold  - (base**emax)*(1-eps)
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015
Purpose:
 DLAMCHF77 determines double precision machine parameters.
Parameters
[in]CMACH
          Specifies the value to be returned by DLAMCH:
          = 'E' or 'e',   DLAMCH := eps
          = 'S' or 's ,   DLAMCH := sfmin
          = 'B' or 'b',   DLAMCH := base
          = 'P' or 'p',   DLAMCH := eps*base
          = 'N' or 'n',   DLAMCH := t
          = 'R' or 'r',   DLAMCH := rnd
          = 'M' or 'm',   DLAMCH := emin
          = 'U' or 'u',   DLAMCH := rmin
          = 'L' or 'l',   DLAMCH := emax
          = 'O' or 'o',   DLAMCH := rmax
          where
          eps   = relative machine precision
          sfmin = safe minimum, such that 1/sfmin does not overflow
          base  = base of the machine
          prec  = eps*base
          t     = number of (base) digits in the mantissa
          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise
          emin  = minimum exponent before (gradual) underflow
          rmin  = underflow threshold - base**(emin-1)
          emax  = largest exponent before overflow
          rmax  = overflow threshold  - (base**emax)*(1-eps)
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
April 2012

Definition at line 65 of file dlamch.f.

65 *
66 * -- LAPACK auxiliary routine (version 3.6.0) --
67 * -- LAPACK is a software package provided by Univ. of Tennessee, --
68 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
69 * November 2015
70 *
71 * .. Scalar Arguments ..
72  CHARACTER cmach
73 * ..
74 *
75 * =====================================================================
76 *
77 * .. Parameters ..
78  DOUBLE PRECISION one, zero
79  parameter( one = 1.0d+0, zero = 0.0d+0 )
80 * ..
81 * .. Local Scalars ..
82  DOUBLE PRECISION rnd, eps, sfmin, small, rmach
83 * ..
84 * .. External Functions ..
85  LOGICAL lsame
86  EXTERNAL lsame
87 * ..
88 * .. Intrinsic Functions ..
89  INTRINSIC digits, epsilon, huge, maxexponent,
90  $ minexponent, radix, tiny
91 * ..
92 * .. Executable Statements ..
93 *
94 *
95 * Assume rounding, not chopping. Always.
96 *
97  rnd = one
98 *
99  IF( one.EQ.rnd ) THEN
100  eps = epsilon(zero) * 0.5
101  ELSE
102  eps = epsilon(zero)
103  END IF
104 *
105  IF( lsame( cmach, 'E' ) ) THEN
106  rmach = eps
107  ELSE IF( lsame( cmach, 'S' ) ) THEN
108  sfmin = tiny(zero)
109  small = one / huge(zero)
110  IF( small.GE.sfmin ) THEN
111 *
112 * Use SMALL plus a bit, to avoid the possibility of rounding
113 * causing overflow when computing 1/sfmin.
114 *
115  sfmin = small*( one+eps )
116  END IF
117  rmach = sfmin
118  ELSE IF( lsame( cmach, 'B' ) ) THEN
119  rmach = radix(zero)
120  ELSE IF( lsame( cmach, 'P' ) ) THEN
121  rmach = eps * radix(zero)
122  ELSE IF( lsame( cmach, 'N' ) ) THEN
123  rmach = digits(zero)
124  ELSE IF( lsame( cmach, 'R' ) ) THEN
125  rmach = rnd
126  ELSE IF( lsame( cmach, 'M' ) ) THEN
127  rmach = minexponent(zero)
128  ELSE IF( lsame( cmach, 'U' ) ) THEN
129  rmach = tiny(zero)
130  ELSE IF( lsame( cmach, 'L' ) ) THEN
131  rmach = maxexponent(zero)
132  ELSE IF( lsame( cmach, 'O' ) ) THEN
133  rmach = huge(zero)
134  ELSE
135  rmach = zero
136  END IF
137 *
138  dlamch = rmach
139  RETURN
140 *
141 * End of DLAMCH
142 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the caller graph for this function:

integer function dlaneg ( integer  N,
double precision, dimension( * )  D,
double precision, dimension( * )  LLD,
double precision  SIGMA,
double precision  PIVMIN,
integer  R 
)

DLANEG computes the Sturm count.

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

Purpose:
 DLANEG computes the Sturm count, the number of negative pivots
 encountered while factoring tridiagonal T - sigma I = L D L^T.
 This implementation works directly on the factors without forming
 the tridiagonal matrix T.  The Sturm count is also the number of
 eigenvalues of T less than sigma.

 This routine is called from DLARRB.

 The current routine does not use the PIVMIN parameter but rather
 requires IEEE-754 propagation of Infinities and NaNs.  This
 routine also has no input range restrictions but does require
 default exception handling such that x/0 produces Inf when x is
 non-zero, and Inf/Inf produces NaN.  For more information, see:

   Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in
   Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on
   Scientific Computing, v28, n5, 2006.  DOI 10.1137/050641624
   (Tech report version in LAWN 172 with the same title.)
Parameters
[in]N
          N is INTEGER
          The order of the matrix.
[in]D
          D is DOUBLE PRECISION array, dimension (N)
          The N diagonal elements of the diagonal matrix D.
[in]LLD
          LLD is DOUBLE PRECISION array, dimension (N-1)
          The (N-1) elements L(i)*L(i)*D(i).
[in]SIGMA
          SIGMA is DOUBLE PRECISION
          Shift amount in T - sigma I = L D L^T.
[in]PIVMIN
          PIVMIN is DOUBLE PRECISION
          The minimum pivot in the Sturm sequence.  May be used
          when zero pivots are encountered on non-IEEE-754
          architectures.
[in]R
          R is INTEGER
          The twist index for the twisted factorization that is used
          for the negcount.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Contributors:
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA
Jason Riedy, University of California, Berkeley, USA

Definition at line 120 of file dlaneg.f.

120 *
121 * -- LAPACK auxiliary routine (version 3.4.2) --
122 * -- LAPACK is a software package provided by Univ. of Tennessee, --
123 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
124 * September 2012
125 *
126 * .. Scalar Arguments ..
127  INTEGER n, r
128  DOUBLE PRECISION pivmin, sigma
129 * ..
130 * .. Array Arguments ..
131  DOUBLE PRECISION d( * ), lld( * )
132 * ..
133 *
134 * =====================================================================
135 *
136 * .. Parameters ..
137  DOUBLE PRECISION zero, one
138  parameter( zero = 0.0d0, one = 1.0d0 )
139 * Some architectures propagate Infinities and NaNs very slowly, so
140 * the code computes counts in BLKLEN chunks. Then a NaN can
141 * propagate at most BLKLEN columns before being detected. This is
142 * not a general tuning parameter; it needs only to be just large
143 * enough that the overhead is tiny in common cases.
144  INTEGER blklen
145  parameter( blklen = 128 )
146 * ..
147 * .. Local Scalars ..
148  INTEGER bj, j, neg1, neg2, negcnt
149  DOUBLE PRECISION bsav, dminus, dplus, gamma, p, t, tmp
150  LOGICAL sawnan
151 * ..
152 * .. Intrinsic Functions ..
153  INTRINSIC min, max
154 * ..
155 * .. External Functions ..
156  LOGICAL disnan
157  EXTERNAL disnan
158 * ..
159 * .. Executable Statements ..
160 
161  negcnt = 0
162 
163 * I) upper part: L D L^T - SIGMA I = L+ D+ L+^T
164  t = -sigma
165  DO 210 bj = 1, r-1, blklen
166  neg1 = 0
167  bsav = t
168  DO 21 j = bj, min(bj+blklen-1, r-1)
169  dplus = d( j ) + t
170  IF( dplus.LT.zero ) neg1 = neg1 + 1
171  tmp = t / dplus
172  t = tmp * lld( j ) - sigma
173  21 CONTINUE
174  sawnan = disnan( t )
175 * Run a slower version of the above loop if a NaN is detected.
176 * A NaN should occur only with a zero pivot after an infinite
177 * pivot. In that case, substituting 1 for T/DPLUS is the
178 * correct limit.
179  IF( sawnan ) THEN
180  neg1 = 0
181  t = bsav
182  DO 22 j = bj, min(bj+blklen-1, r-1)
183  dplus = d( j ) + t
184  IF( dplus.LT.zero ) neg1 = neg1 + 1
185  tmp = t / dplus
186  IF (disnan(tmp)) tmp = one
187  t = tmp * lld(j) - sigma
188  22 CONTINUE
189  END IF
190  negcnt = negcnt + neg1
191  210 CONTINUE
192 *
193 * II) lower part: L D L^T - SIGMA I = U- D- U-^T
194  p = d( n ) - sigma
195  DO 230 bj = n-1, r, -blklen
196  neg2 = 0
197  bsav = p
198  DO 23 j = bj, max(bj-blklen+1, r), -1
199  dminus = lld( j ) + p
200  IF( dminus.LT.zero ) neg2 = neg2 + 1
201  tmp = p / dminus
202  p = tmp * d( j ) - sigma
203  23 CONTINUE
204  sawnan = disnan( p )
205 * As above, run a slower version that substitutes 1 for Inf/Inf.
206 *
207  IF( sawnan ) THEN
208  neg2 = 0
209  p = bsav
210  DO 24 j = bj, max(bj-blklen+1, r), -1
211  dminus = lld( j ) + p
212  IF( dminus.LT.zero ) neg2 = neg2 + 1
213  tmp = p / dminus
214  IF (disnan(tmp)) tmp = one
215  p = tmp * d(j) - sigma
216  24 CONTINUE
217  END IF
218  negcnt = negcnt + neg2
219  230 CONTINUE
220 *
221 * III) Twist index
222 * T was shifted by SIGMA initially.
223  gamma = (t + sigma) + p
224  IF( gamma.LT.zero ) negcnt = negcnt+1
225 
226  dlaneg = negcnt
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:61
integer function dlaneg(N, D, LLD, SIGMA, PIVMIN, R)
DLANEG computes the Sturm count.
Definition: dlaneg.f:120
double precision function dlanst ( character  NORM,
integer  N,
double precision, dimension( * )  D,
double precision, dimension( * )  E 
)

DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric tridiagonal matrix.

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

Purpose:
 DLANST  returns the value of the one norm,  or the Frobenius norm, or
 the  infinity norm,  or the  element of  largest absolute value  of a
 real symmetric tridiagonal matrix A.
Returns
DLANST
    DLANST = ( max(abs(A(i,j))), NORM = 'M' or 'm'
             (
             ( norm1(A),         NORM = '1', 'O' or 'o'
             (
             ( normI(A),         NORM = 'I' or 'i'
             (
             ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

 where  norm1  denotes the  one norm of a matrix (maximum column sum),
 normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
 normF  denotes the  Frobenius norm of a matrix (square root of sum of
 squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.
Parameters
[in]NORM
          NORM is CHARACTER*1
          Specifies the value to be returned in DLANST as described
          above.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.  When N = 0, DLANST is
          set to zero.
[in]D
          D is DOUBLE PRECISION array, dimension (N)
          The diagonal elements of A.
[in]E
          E is DOUBLE PRECISION array, dimension (N-1)
          The (n-1) sub-diagonal or super-diagonal elements of A.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 102 of file dlanst.f.

102 *
103 * -- LAPACK auxiliary routine (version 3.4.2) --
104 * -- LAPACK is a software package provided by Univ. of Tennessee, --
105 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
106 * September 2012
107 *
108 * .. Scalar Arguments ..
109  CHARACTER norm
110  INTEGER n
111 * ..
112 * .. Array Arguments ..
113  DOUBLE PRECISION d( * ), e( * )
114 * ..
115 *
116 * =====================================================================
117 *
118 * .. Parameters ..
119  DOUBLE PRECISION one, zero
120  parameter( one = 1.0d+0, zero = 0.0d+0 )
121 * ..
122 * .. Local Scalars ..
123  INTEGER i
124  DOUBLE PRECISION anorm, scale, sum
125 * ..
126 * .. External Functions ..
127  LOGICAL lsame, disnan
128  EXTERNAL lsame, disnan
129 * ..
130 * .. External Subroutines ..
131  EXTERNAL dlassq
132 * ..
133 * .. Intrinsic Functions ..
134  INTRINSIC abs, sqrt
135 * ..
136 * .. Executable Statements ..
137 *
138  IF( n.LE.0 ) THEN
139  anorm = zero
140  ELSE IF( lsame( norm, 'M' ) ) THEN
141 *
142 * Find max(abs(A(i,j))).
143 *
144  anorm = abs( d( n ) )
145  DO 10 i = 1, n - 1
146  sum = abs( d( i ) )
147  IF( anorm .LT. sum .OR. disnan( sum ) ) anorm = sum
148  sum = abs( e( i ) )
149  IF( anorm .LT. sum .OR. disnan( sum ) ) anorm = sum
150  10 CONTINUE
151  ELSE IF( lsame( norm, 'O' ) .OR. norm.EQ.'1' .OR.
152  $ lsame( norm, 'I' ) ) THEN
153 *
154 * Find norm1(A).
155 *
156  IF( n.EQ.1 ) THEN
157  anorm = abs( d( 1 ) )
158  ELSE
159  anorm = abs( d( 1 ) )+abs( e( 1 ) )
160  sum = abs( e( n-1 ) )+abs( d( n ) )
161  IF( anorm .LT. sum .OR. disnan( sum ) ) anorm = sum
162  DO 20 i = 2, n - 1
163  sum = abs( d( i ) )+abs( e( i ) )+abs( e( i-1 ) )
164  IF( anorm .LT. sum .OR. disnan( sum ) ) anorm = sum
165  20 CONTINUE
166  END IF
167  ELSE IF( ( lsame( norm, 'F' ) ) .OR. ( lsame( norm, 'E' ) ) ) THEN
168 *
169 * Find normF(A).
170 *
171  scale = zero
172  sum = one
173  IF( n.GT.1 ) THEN
174  CALL dlassq( n-1, e, 1, scale, sum )
175  sum = 2*sum
176  END IF
177  CALL dlassq( n, d, 1, scale, sum )
178  anorm = scale*sqrt( sum )
179  END IF
180 *
181  dlanst = anorm
182  RETURN
183 *
184 * End of DLANST
185 *
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:61
double precision function dlanst(NORM, N, D, E)
DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric tridiagonal matrix.
Definition: dlanst.f:102
subroutine dlassq(N, X, INCX, SCALE, SUMSQ)
DLASSQ updates a sum of squares represented in scaled form.
Definition: dlassq.f:105
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

double precision function dlapy2 ( double precision  X,
double precision  Y 
)

DLAPY2 returns sqrt(x2+y2).

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

Purpose:
 DLAPY2 returns sqrt(x**2+y**2), taking care not to cause unnecessary
 overflow.
Parameters
[in]X
          X is DOUBLE PRECISION
[in]Y
          Y is DOUBLE PRECISION
          X and Y specify the values x and y.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 65 of file dlapy2.f.

65 *
66 * -- LAPACK auxiliary routine (version 3.4.2) --
67 * -- LAPACK is a software package provided by Univ. of Tennessee, --
68 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
69 * September 2012
70 *
71 * .. Scalar Arguments ..
72  DOUBLE PRECISION x, y
73 * ..
74 *
75 * =====================================================================
76 *
77 * .. Parameters ..
78  DOUBLE PRECISION zero
79  parameter( zero = 0.0d0 )
80  DOUBLE PRECISION one
81  parameter( one = 1.0d0 )
82 * ..
83 * .. Local Scalars ..
84  DOUBLE PRECISION w, xabs, yabs, z
85 * ..
86 * .. Intrinsic Functions ..
87  INTRINSIC abs, max, min, sqrt
88 * ..
89 * .. Executable Statements ..
90 *
91  xabs = abs( x )
92  yabs = abs( y )
93  w = max( xabs, yabs )
94  z = min( xabs, yabs )
95  IF( z.EQ.zero ) THEN
96  dlapy2 = w
97  ELSE
98  dlapy2 = w*sqrt( one+( z / w )**2 )
99  END IF
100  RETURN
101 *
102 * End of DLAPY2
103 *
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
Definition: dlapy2.f:65
double precision function dlapy3 ( double precision  X,
double precision  Y,
double precision  Z 
)

DLAPY3 returns sqrt(x2+y2+z2).

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

Purpose:
 DLAPY3 returns sqrt(x**2+y**2+z**2), taking care not to cause
 unnecessary overflow.
Parameters
[in]X
          X is DOUBLE PRECISION
[in]Y
          Y is DOUBLE PRECISION
[in]Z
          Z is DOUBLE PRECISION
          X, Y and Z specify the values x, y and z.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 70 of file dlapy3.f.

70 *
71 * -- LAPACK auxiliary routine (version 3.4.2) --
72 * -- LAPACK is a software package provided by Univ. of Tennessee, --
73 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
74 * September 2012
75 *
76 * .. Scalar Arguments ..
77  DOUBLE PRECISION x, y, z
78 * ..
79 *
80 * =====================================================================
81 *
82 * .. Parameters ..
83  DOUBLE PRECISION zero
84  parameter( zero = 0.0d0 )
85 * ..
86 * .. Local Scalars ..
87  DOUBLE PRECISION w, xabs, yabs, zabs
88 * ..
89 * .. Intrinsic Functions ..
90  INTRINSIC abs, max, sqrt
91 * ..
92 * .. Executable Statements ..
93 *
94  xabs = abs( x )
95  yabs = abs( y )
96  zabs = abs( z )
97  w = max( xabs, yabs, zabs )
98  IF( w.EQ.zero ) THEN
99 * W can be zero for max(0,nan,0)
100 * adding all three entries together will make sure
101 * NaN will not disappear.
102  dlapy3 = xabs + yabs + zabs
103  ELSE
104  dlapy3 = w*sqrt( ( xabs / w )**2+( yabs / w )**2+
105  $ ( zabs / w )**2 )
106  END IF
107  RETURN
108 *
109 * End of DLAPY3
110 *
double precision function dlapy3(X, Y, Z)
DLAPY3 returns sqrt(x2+y2+z2).
Definition: dlapy3.f:70
subroutine dlarnv ( integer  IDIST,
integer, dimension( 4 )  ISEED,
integer  N,
double precision, dimension( * )  X 
)

DLARNV returns a vector of random numbers from a uniform or normal distribution.

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

Purpose:
 DLARNV returns a vector of n random real numbers from a uniform or
 normal distribution.
Parameters
[in]IDIST
          IDIST is INTEGER
          Specifies the distribution of the random numbers:
          = 1:  uniform (0,1)
          = 2:  uniform (-1,1)
          = 3:  normal (0,1)
[in,out]ISEED
          ISEED is INTEGER array, dimension (4)
          On entry, the seed of the random number generator; the array
          elements must be between 0 and 4095, and ISEED(4) must be
          odd.
          On exit, the seed is updated.
[in]N
          N is INTEGER
          The number of random numbers to be generated.
[out]X
          X is DOUBLE PRECISION array, dimension (N)
          The generated random numbers.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
  This routine calls the auxiliary routine DLARUV to generate random
  real numbers from a uniform (0,1) distribution, in batches of up to
  128 using vectorisable code. The Box-Muller method is used to
  transform numbers from a uniform to a normal distribution.

Definition at line 99 of file dlarnv.f.

99 *
100 * -- LAPACK auxiliary routine (version 3.4.2) --
101 * -- LAPACK is a software package provided by Univ. of Tennessee, --
102 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
103 * September 2012
104 *
105 * .. Scalar Arguments ..
106  INTEGER idist, n
107 * ..
108 * .. Array Arguments ..
109  INTEGER iseed( 4 )
110  DOUBLE PRECISION x( * )
111 * ..
112 *
113 * =====================================================================
114 *
115 * .. Parameters ..
116  DOUBLE PRECISION one, two
117  parameter( one = 1.0d+0, two = 2.0d+0 )
118  INTEGER lv
119  parameter( lv = 128 )
120  DOUBLE PRECISION twopi
121  parameter( twopi = 6.2831853071795864769252867663d+0 )
122 * ..
123 * .. Local Scalars ..
124  INTEGER i, il, il2, iv
125 * ..
126 * .. Local Arrays ..
127  DOUBLE PRECISION u( lv )
128 * ..
129 * .. Intrinsic Functions ..
130  INTRINSIC cos, log, min, sqrt
131 * ..
132 * .. External Subroutines ..
133  EXTERNAL dlaruv
134 * ..
135 * .. Executable Statements ..
136 *
137  DO 40 iv = 1, n, lv / 2
138  il = min( lv / 2, n-iv+1 )
139  IF( idist.EQ.3 ) THEN
140  il2 = 2*il
141  ELSE
142  il2 = il
143  END IF
144 *
145 * Call DLARUV to generate IL2 numbers from a uniform (0,1)
146 * distribution (IL2 <= LV)
147 *
148  CALL dlaruv( iseed, il2, u )
149 *
150  IF( idist.EQ.1 ) THEN
151 *
152 * Copy generated numbers
153 *
154  DO 10 i = 1, il
155  x( iv+i-1 ) = u( i )
156  10 CONTINUE
157  ELSE IF( idist.EQ.2 ) THEN
158 *
159 * Convert generated numbers to uniform (-1,1) distribution
160 *
161  DO 20 i = 1, il
162  x( iv+i-1 ) = two*u( i ) - one
163  20 CONTINUE
164  ELSE IF( idist.EQ.3 ) THEN
165 *
166 * Convert generated numbers to normal (0,1) distribution
167 *
168  DO 30 i = 1, il
169  x( iv+i-1 ) = sqrt( -two*log( u( 2*i-1 ) ) )*
170  $ cos( twopi*u( 2*i ) )
171  30 CONTINUE
172  END IF
173  40 CONTINUE
174  RETURN
175 *
176 * End of DLARNV
177 *
subroutine dlaruv(ISEED, N, X)
DLARUV returns a vector of n random real numbers from a uniform distribution.
Definition: dlaruv.f:97

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dlarra ( integer  N,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( * )  E2,
double precision  SPLTOL,
double precision  TNRM,
integer  NSPLIT,
integer, dimension( * )  ISPLIT,
integer  INFO 
)

DLARRA computes the splitting points with the specified threshold.

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

Purpose:
 Compute the splitting points with threshold SPLTOL.
 DLARRA sets any "small" off-diagonal elements to zero.
Parameters
[in]N
          N is INTEGER
          The order of the matrix. N > 0.
[in]D
          D is DOUBLE PRECISION array, dimension (N)
          On entry, the N diagonal elements of the tridiagonal
          matrix T.
[in,out]E
          E is DOUBLE PRECISION array, dimension (N)
          On entry, the first (N-1) entries contain the subdiagonal
          elements of the tridiagonal matrix T; E(N) need not be set.
          On exit, the entries E( ISPLIT( I ) ), 1 <= I <= NSPLIT,
          are set to zero, the other entries of E are untouched.
[in,out]E2
          E2 is DOUBLE PRECISION array, dimension (N)
          On entry, the first (N-1) entries contain the SQUARES of the
          subdiagonal elements of the tridiagonal matrix T;
          E2(N) need not be set.
          On exit, the entries E2( ISPLIT( I ) ),
          1 <= I <= NSPLIT, have been set to zero
[in]SPLTOL
          SPLTOL is DOUBLE PRECISION
          The threshold for splitting. Two criteria can be used:
          SPLTOL<0 : criterion based on absolute off-diagonal value
          SPLTOL>0 : criterion that preserves relative accuracy
[in]TNRM
          TNRM is DOUBLE PRECISION
          The norm of the matrix.
[out]NSPLIT
          NSPLIT is INTEGER
          The number of blocks T splits into. 1 <= NSPLIT <= N.
[out]ISPLIT
          ISPLIT is INTEGER array, dimension (N)
          The splitting points, at which T breaks up into blocks.
          The first block consists of rows/columns 1 to ISPLIT(1),
          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
          etc., and the NSPLIT-th consists of rows/columns
          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Contributors:
Beresford Parlett, University of California, Berkeley, USA
Jim Demmel, University of California, Berkeley, USA
Inderjit Dhillon, University of Texas, Austin, USA
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA

Definition at line 138 of file dlarra.f.

138 *
139 * -- LAPACK auxiliary routine (version 3.4.2) --
140 * -- LAPACK is a software package provided by Univ. of Tennessee, --
141 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
142 * September 2012
143 *
144 * .. Scalar Arguments ..
145  INTEGER info, n, nsplit
146  DOUBLE PRECISION spltol, tnrm
147 * ..
148 * .. Array Arguments ..
149  INTEGER isplit( * )
150  DOUBLE PRECISION d( * ), e( * ), e2( * )
151 * ..
152 *
153 * =====================================================================
154 *
155 * .. Parameters ..
156  DOUBLE PRECISION zero
157  parameter( zero = 0.0d0 )
158 * ..
159 * .. Local Scalars ..
160  INTEGER i
161  DOUBLE PRECISION eabs, tmp1
162 
163 * ..
164 * .. Intrinsic Functions ..
165  INTRINSIC abs
166 * ..
167 * .. Executable Statements ..
168 *
169  info = 0
170 
171 * Compute splitting points
172  nsplit = 1
173  IF(spltol.LT.zero) THEN
174 * Criterion based on absolute off-diagonal value
175  tmp1 = abs(spltol)* tnrm
176  DO 9 i = 1, n-1
177  eabs = abs( e(i) )
178  IF( eabs .LE. tmp1) THEN
179  e(i) = zero
180  e2(i) = zero
181  isplit( nsplit ) = i
182  nsplit = nsplit + 1
183  END IF
184  9 CONTINUE
185  ELSE
186 * Criterion that guarantees relative accuracy
187  DO 10 i = 1, n-1
188  eabs = abs( e(i) )
189  IF( eabs .LE. spltol * sqrt(abs(d(i)))*sqrt(abs(d(i+1))) )
190  $ THEN
191  e(i) = zero
192  e2(i) = zero
193  isplit( nsplit ) = i
194  nsplit = nsplit + 1
195  END IF
196  10 CONTINUE
197  ENDIF
198  isplit( nsplit ) = n
199 
200  RETURN
201 *
202 * End of DLARRA
203 *

Here is the caller graph for this function:

subroutine dlarrb ( integer  N,
double precision, dimension( * )  D,
double precision, dimension( * )  LLD,
integer  IFIRST,
integer  ILAST,
double precision  RTOL1,
double precision  RTOL2,
integer  OFFSET,
double precision, dimension( * )  W,
double precision, dimension( * )  WGAP,
double precision, dimension( * )  WERR,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
double precision  PIVMIN,
double precision  SPDIAM,
integer  TWIST,
integer  INFO 
)

DLARRB provides limited bisection to locate eigenvalues for more accuracy.

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

Purpose:
 Given the relatively robust representation(RRR) L D L^T, DLARRB
 does "limited" bisection to refine the eigenvalues of L D L^T,
 W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
 guesses for these eigenvalues are input in W, the corresponding estimate
 of the error in these guesses and their gaps are input in WERR
 and WGAP, respectively. During bisection, intervals
 [left, right] are maintained by storing their mid-points and
 semi-widths in the arrays W and WERR respectively.
Parameters
[in]N
          N is INTEGER
          The order of the matrix.
[in]D
          D is DOUBLE PRECISION array, dimension (N)
          The N diagonal elements of the diagonal matrix D.
[in]LLD
          LLD is DOUBLE PRECISION array, dimension (N-1)
          The (N-1) elements L(i)*L(i)*D(i).
[in]IFIRST
          IFIRST is INTEGER
          The index of the first eigenvalue to be computed.
[in]ILAST
          ILAST is INTEGER
          The index of the last eigenvalue to be computed.
[in]RTOL1
          RTOL1 is DOUBLE PRECISION
[in]RTOL2
          RTOL2 is DOUBLE PRECISION
          Tolerance for the convergence of the bisection intervals.
          An interval [LEFT,RIGHT] has converged if
          RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
          where GAP is the (estimated) distance to the nearest
          eigenvalue.
[in]OFFSET
          OFFSET is INTEGER
          Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET
          through ILAST-OFFSET elements of these arrays are to be used.
[in,out]W
          W is DOUBLE PRECISION array, dimension (N)
          On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
          estimates of the eigenvalues of L D L^T indexed IFIRST throug
          ILAST.
          On output, these estimates are refined.
[in,out]WGAP
          WGAP is DOUBLE PRECISION array, dimension (N-1)
          On input, the (estimated) gaps between consecutive
          eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between
          eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST
          then WGAP(IFIRST-OFFSET) must be set to ZERO.
          On output, these gaps are refined.
[in,out]WERR
          WERR is DOUBLE PRECISION array, dimension (N)
          On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
          the errors in the estimates of the corresponding elements in W.
          On output, these errors are refined.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (2*N)
          Workspace.
[out]IWORK
          IWORK is INTEGER array, dimension (2*N)
          Workspace.
[in]PIVMIN
          PIVMIN is DOUBLE PRECISION
          The minimum pivot in the Sturm sequence.
[in]SPDIAM
          SPDIAM is DOUBLE PRECISION
          The spectral diameter of the matrix.
[in]TWIST
          TWIST is INTEGER
          The twist index for the twisted factorization that is used
          for the negcount.
          TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T
          TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T
          TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r)
[out]INFO
          INFO is INTEGER
          Error flag.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Contributors:
Beresford Parlett, University of California, Berkeley, USA
Jim Demmel, University of California, Berkeley, USA
Inderjit Dhillon, University of Texas, Austin, USA
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA

Definition at line 198 of file dlarrb.f.

198 *
199 * -- LAPACK auxiliary routine (version 3.4.2) --
200 * -- LAPACK is a software package provided by Univ. of Tennessee, --
201 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
202 * September 2012
203 *
204 * .. Scalar Arguments ..
205  INTEGER ifirst, ilast, info, n, offset, twist
206  DOUBLE PRECISION pivmin, rtol1, rtol2, spdiam
207 * ..
208 * .. Array Arguments ..
209  INTEGER iwork( * )
210  DOUBLE PRECISION d( * ), lld( * ), w( * ),
211  $ werr( * ), wgap( * ), work( * )
212 * ..
213 *
214 * =====================================================================
215 *
216 * .. Parameters ..
217  DOUBLE PRECISION zero, two, half
218  parameter( zero = 0.0d0, two = 2.0d0,
219  $ half = 0.5d0 )
220  INTEGER maxitr
221 * ..
222 * .. Local Scalars ..
223  INTEGER i, i1, ii, ip, iter, k, negcnt, next, nint,
224  $ olnint, prev, r
225  DOUBLE PRECISION back, cvrgd, gap, left, lgap, mid, mnwdth,
226  $ rgap, right, tmp, width
227 * ..
228 * .. External Functions ..
229  INTEGER dlaneg
230  EXTERNAL dlaneg
231 *
232 * ..
233 * .. Intrinsic Functions ..
234  INTRINSIC abs, max, min
235 * ..
236 * .. Executable Statements ..
237 *
238  info = 0
239 *
240  maxitr = int( ( log( spdiam+pivmin )-log( pivmin ) ) /
241  $ log( two ) ) + 2
242  mnwdth = two * pivmin
243 *
244  r = twist
245  IF((r.LT.1).OR.(r.GT.n)) r = n
246 *
247 * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
248 * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
249 * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
250 * for an unconverged interval is set to the index of the next unconverged
251 * interval, and is -1 or 0 for a converged interval. Thus a linked
252 * list of unconverged intervals is set up.
253 *
254  i1 = ifirst
255 * The number of unconverged intervals
256  nint = 0
257 * The last unconverged interval found
258  prev = 0
259 
260  rgap = wgap( i1-offset )
261  DO 75 i = i1, ilast
262  k = 2*i
263  ii = i - offset
264  left = w( ii ) - werr( ii )
265  right = w( ii ) + werr( ii )
266  lgap = rgap
267  rgap = wgap( ii )
268  gap = min( lgap, rgap )
269 
270 * Make sure that [LEFT,RIGHT] contains the desired eigenvalue
271 * Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT
272 *
273 * Do while( NEGCNT(LEFT).GT.I-1 )
274 *
275  back = werr( ii )
276  20 CONTINUE
277  negcnt = dlaneg( n, d, lld, left, pivmin, r )
278  IF( negcnt.GT.i-1 ) THEN
279  left = left - back
280  back = two*back
281  GO TO 20
282  END IF
283 *
284 * Do while( NEGCNT(RIGHT).LT.I )
285 * Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT
286 *
287  back = werr( ii )
288  50 CONTINUE
289 
290  negcnt = dlaneg( n, d, lld, right, pivmin, r )
291  IF( negcnt.LT.i ) THEN
292  right = right + back
293  back = two*back
294  GO TO 50
295  END IF
296  width = half*abs( left - right )
297  tmp = max( abs( left ), abs( right ) )
298  cvrgd = max(rtol1*gap,rtol2*tmp)
299  IF( width.LE.cvrgd .OR. width.LE.mnwdth ) THEN
300 * This interval has already converged and does not need refinement.
301 * (Note that the gaps might change through refining the
302 * eigenvalues, however, they can only get bigger.)
303 * Remove it from the list.
304  iwork( k-1 ) = -1
305 * Make sure that I1 always points to the first unconverged interval
306  IF((i.EQ.i1).AND.(i.LT.ilast)) i1 = i + 1
307  IF((prev.GE.i1).AND.(i.LE.ilast)) iwork( 2*prev-1 ) = i + 1
308  ELSE
309 * unconverged interval found
310  prev = i
311  nint = nint + 1
312  iwork( k-1 ) = i + 1
313  iwork( k ) = negcnt
314  END IF
315  work( k-1 ) = left
316  work( k ) = right
317  75 CONTINUE
318 
319 *
320 * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
321 * and while (ITER.LT.MAXITR)
322 *
323  iter = 0
324  80 CONTINUE
325  prev = i1 - 1
326  i = i1
327  olnint = nint
328 
329  DO 100 ip = 1, olnint
330  k = 2*i
331  ii = i - offset
332  rgap = wgap( ii )
333  lgap = rgap
334  IF(ii.GT.1) lgap = wgap( ii-1 )
335  gap = min( lgap, rgap )
336  next = iwork( k-1 )
337  left = work( k-1 )
338  right = work( k )
339  mid = half*( left + right )
340 
341 * semiwidth of interval
342  width = right - mid
343  tmp = max( abs( left ), abs( right ) )
344  cvrgd = max(rtol1*gap,rtol2*tmp)
345  IF( ( width.LE.cvrgd ) .OR. ( width.LE.mnwdth ).OR.
346  $ ( iter.EQ.maxitr ) )THEN
347 * reduce number of unconverged intervals
348  nint = nint - 1
349 * Mark interval as converged.
350  iwork( k-1 ) = 0
351  IF( i1.EQ.i ) THEN
352  i1 = next
353  ELSE
354 * Prev holds the last unconverged interval previously examined
355  IF(prev.GE.i1) iwork( 2*prev-1 ) = next
356  END IF
357  i = next
358  GO TO 100
359  END IF
360  prev = i
361 *
362 * Perform one bisection step
363 *
364  negcnt = dlaneg( n, d, lld, mid, pivmin, r )
365  IF( negcnt.LE.i-1 ) THEN
366  work( k-1 ) = mid
367  ELSE
368  work( k ) = mid
369  END IF
370  i = next
371  100 CONTINUE
372  iter = iter + 1
373 * do another loop if there are still unconverged intervals
374 * However, in the last iteration, all intervals are accepted
375 * since this is the best we can do.
376  IF( ( nint.GT.0 ).AND.(iter.LE.maxitr) ) GO TO 80
377 *
378 *
379 * At this point, all the intervals have converged
380  DO 110 i = ifirst, ilast
381  k = 2*i
382  ii = i - offset
383 * All intervals marked by '0' have been refined.
384  IF( iwork( k-1 ).EQ.0 ) THEN
385  w( ii ) = half*( work( k-1 )+work( k ) )
386  werr( ii ) = work( k ) - w( ii )
387  END IF
388  110 CONTINUE
389 *
390  DO 111 i = ifirst+1, ilast
391  k = 2*i
392  ii = i - offset
393  wgap( ii-1 ) = max( zero,
394  $ w(ii) - werr(ii) - w( ii-1 ) - werr( ii-1 ))
395  111 CONTINUE
396 
397  RETURN
398 *
399 * End of DLARRB
400 *
integer function dlaneg(N, D, LLD, SIGMA, PIVMIN, R)
DLANEG computes the Sturm count.
Definition: dlaneg.f:120

Here is the caller graph for this function:

subroutine dlarrc ( character  JOBT,
integer  N,
double precision  VL,
double precision  VU,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision  PIVMIN,
integer  EIGCNT,
integer  LCNT,
integer  RCNT,
integer  INFO 
)

DLARRC computes the number of eigenvalues of the symmetric tridiagonal matrix.

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

Purpose:
 Find the number of eigenvalues of the symmetric tridiagonal matrix T
 that are in the interval (VL,VU] if JOBT = 'T', and of L D L^T
 if JOBT = 'L'.
Parameters
[in]JOBT
          JOBT is CHARACTER*1
          = 'T':  Compute Sturm count for matrix T.
          = 'L':  Compute Sturm count for matrix L D L^T.
[in]N
          N is INTEGER
          The order of the matrix. N > 0.
[in]VL
          VL is DOUBLE PRECISION
[in]VU
          VU is DOUBLE PRECISION
          The lower and upper bounds for the eigenvalues.
[in]D
          D is DOUBLE PRECISION array, dimension (N)
          JOBT = 'T': The N diagonal elements of the tridiagonal matrix T.
          JOBT = 'L': The N diagonal elements of the diagonal matrix D.
[in]E
          E is DOUBLE PRECISION array, dimension (N)
          JOBT = 'T': The N-1 offdiagonal elements of the matrix T.
          JOBT = 'L': The N-1 offdiagonal elements of the matrix L.
[in]PIVMIN
          PIVMIN is DOUBLE PRECISION
          The minimum pivot in the Sturm sequence for T.
[out]EIGCNT
          EIGCNT is INTEGER
          The number of eigenvalues of the symmetric tridiagonal matrix T
          that are in the interval (VL,VU]
[out]LCNT
          LCNT is INTEGER
[out]RCNT
          RCNT is INTEGER
          The left and right negcounts of the interval.
[out]INFO
          INFO is INTEGER
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Contributors:
Beresford Parlett, University of California, Berkeley, USA
Jim Demmel, University of California, Berkeley, USA
Inderjit Dhillon, University of Texas, Austin, USA
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA

Definition at line 138 of file dlarrc.f.

138 *
139 * -- LAPACK auxiliary routine (version 3.4.2) --
140 * -- LAPACK is a software package provided by Univ. of Tennessee, --
141 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
142 * September 2012
143 *
144 * .. Scalar Arguments ..
145  CHARACTER jobt
146  INTEGER eigcnt, info, lcnt, n, rcnt
147  DOUBLE PRECISION pivmin, vl, vu
148 * ..
149 * .. Array Arguments ..
150  DOUBLE PRECISION d( * ), e( * )
151 * ..
152 *
153 * =====================================================================
154 *
155 * .. Parameters ..
156  DOUBLE PRECISION zero
157  parameter( zero = 0.0d0 )
158 * ..
159 * .. Local Scalars ..
160  INTEGER i
161  LOGICAL matt
162  DOUBLE PRECISION lpivot, rpivot, sl, su, tmp, tmp2
163 
164 * ..
165 * .. External Functions ..
166  LOGICAL lsame
167  EXTERNAL lsame
168 * ..
169 * .. Executable Statements ..
170 *
171  info = 0
172  lcnt = 0
173  rcnt = 0
174  eigcnt = 0
175  matt = lsame( jobt, 'T' )
176 
177 
178  IF (matt) THEN
179 * Sturm sequence count on T
180  lpivot = d( 1 ) - vl
181  rpivot = d( 1 ) - vu
182  IF( lpivot.LE.zero ) THEN
183  lcnt = lcnt + 1
184  ENDIF
185  IF( rpivot.LE.zero ) THEN
186  rcnt = rcnt + 1
187  ENDIF
188  DO 10 i = 1, n-1
189  tmp = e(i)**2
190  lpivot = ( d( i+1 )-vl ) - tmp/lpivot
191  rpivot = ( d( i+1 )-vu ) - tmp/rpivot
192  IF( lpivot.LE.zero ) THEN
193  lcnt = lcnt + 1
194  ENDIF
195  IF( rpivot.LE.zero ) THEN
196  rcnt = rcnt + 1
197  ENDIF
198  10 CONTINUE
199  ELSE
200 * Sturm sequence count on L D L^T
201  sl = -vl
202  su = -vu
203  DO 20 i = 1, n - 1
204  lpivot = d( i ) + sl
205  rpivot = d( i ) + su
206  IF( lpivot.LE.zero ) THEN
207  lcnt = lcnt + 1
208  ENDIF
209  IF( rpivot.LE.zero ) THEN
210  rcnt = rcnt + 1
211  ENDIF
212  tmp = e(i) * d(i) * e(i)
213 *
214  tmp2 = tmp / lpivot
215  IF( tmp2.EQ.zero ) THEN
216  sl = tmp - vl
217  ELSE
218  sl = sl*tmp2 - vl
219  END IF
220 *
221  tmp2 = tmp / rpivot
222  IF( tmp2.EQ.zero ) THEN
223  su = tmp - vu
224  ELSE
225  su = su*tmp2 - vu
226  END IF
227  20 CONTINUE
228  lpivot = d( n ) + sl
229  rpivot = d( n ) + su
230  IF( lpivot.LE.zero ) THEN
231  lcnt = lcnt + 1
232  ENDIF
233  IF( rpivot.LE.zero ) THEN
234  rcnt = rcnt + 1
235  ENDIF
236  ENDIF
237  eigcnt = rcnt - lcnt
238 
239  RETURN
240 *
241 * end of DLARRC
242 *
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the caller graph for this function:

subroutine dlarrd ( character  RANGE,
character  ORDER,
integer  N,
double precision  VL,
double precision  VU,
integer  IL,
integer  IU,
double precision, dimension( * )  GERS,
double precision  RELTOL,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( * )  E2,
double precision  PIVMIN,
integer  NSPLIT,
integer, dimension( * )  ISPLIT,
integer  M,
double precision, dimension( * )  W,
double precision, dimension( * )  WERR,
double precision  WL,
double precision  WU,
integer, dimension( * )  IBLOCK,
integer, dimension( * )  INDEXW,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

DLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy.

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

Purpose:
 DLARRD computes the eigenvalues of a symmetric tridiagonal
 matrix T to suitable accuracy. This is an auxiliary code to be
 called from DSTEMR.
 The user may ask for all eigenvalues, all eigenvalues
 in the half-open interval (VL, VU], or the IL-th through IU-th
 eigenvalues.

 To avoid overflow, the matrix must be scaled so that its
 largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
 accuracy, it should not be much smaller than that.

 See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
 Matrix", Report CS41, Computer Science Dept., Stanford
 University, July 21, 1966.
Parameters
[in]RANGE
          RANGE is CHARACTER*1
          = 'A': ("All")   all eigenvalues will be found.
          = 'V': ("Value") all eigenvalues in the half-open interval
                           (VL, VU] will be found.
          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
                           entire matrix) will be found.
[in]ORDER
          ORDER is CHARACTER*1
          = 'B': ("By Block") the eigenvalues will be grouped by
                              split-off block (see IBLOCK, ISPLIT) and
                              ordered from smallest to largest within
                              the block.
          = 'E': ("Entire matrix")
                              the eigenvalues for the entire matrix
                              will be ordered from smallest to
                              largest.
[in]N
          N is INTEGER
          The order of the tridiagonal matrix T.  N >= 0.
[in]VL
          VL is DOUBLE PRECISION
[in]VU
          VU is DOUBLE PRECISION
          If RANGE='V', the lower and upper bounds of the interval to
          be searched for eigenvalues.  Eigenvalues less than or equal
          to VL, or greater than VU, will not be returned.  VL < VU.
          Not referenced if RANGE = 'A' or 'I'.
[in]IL
          IL is INTEGER
[in]IU
          IU is INTEGER
          If RANGE='I', the indices (in ascending order) of the
          smallest and largest eigenvalues 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]GERS
          GERS is DOUBLE PRECISION array, dimension (2*N)
          The N Gerschgorin intervals (the i-th Gerschgorin interval
          is (GERS(2*i-1), GERS(2*i)).
[in]RELTOL
          RELTOL is DOUBLE PRECISION
          The minimum relative width of an interval.  When an interval
          is narrower than RELTOL times the larger (in
          magnitude) endpoint, then it is considered to be
          sufficiently small, i.e., converged.  Note: this should
          always be at least radix*machine epsilon.
[in]D
          D is DOUBLE PRECISION array, dimension (N)
          The n diagonal elements of the tridiagonal matrix T.
[in]E
          E is DOUBLE PRECISION array, dimension (N-1)
          The (n-1) off-diagonal elements of the tridiagonal matrix T.
[in]E2
          E2 is DOUBLE PRECISION array, dimension (N-1)
          The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
[in]PIVMIN
          PIVMIN is DOUBLE PRECISION
          The minimum pivot allowed in the Sturm sequence for T.
[in]NSPLIT
          NSPLIT is INTEGER
          The number of diagonal blocks in the matrix T.
          1 <= NSPLIT <= N.
[in]ISPLIT
          ISPLIT is INTEGER array, dimension (N)
          The splitting points, at which T breaks up into submatrices.
          The first submatrix consists of rows/columns 1 to ISPLIT(1),
          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
          etc., and the NSPLIT-th consists of rows/columns
          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
          (Only the first NSPLIT elements will actually be used, but
          since the user cannot know a priori what value NSPLIT will
          have, N words must be reserved for ISPLIT.)
[out]M
          M is INTEGER
          The actual number of eigenvalues found. 0 <= M <= N.
          (See also the description of INFO=2,3.)
[out]W
          W is DOUBLE PRECISION array, dimension (N)
          On exit, the first M elements of W will contain the
          eigenvalue approximations. DLARRD computes an interval
          I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue
          approximation is given as the interval midpoint
          W(j)= ( a_j + b_j)/2. The corresponding error is bounded by
          WERR(j) = abs( a_j - b_j)/2
[out]WERR
          WERR is DOUBLE PRECISION array, dimension (N)
          The error bound on the corresponding eigenvalue approximation
          in W.
[out]WL
          WL is DOUBLE PRECISION
[out]WU
          WU is DOUBLE PRECISION
          The interval (WL, WU] contains all the wanted eigenvalues.
          If RANGE='V', then WL=VL and WU=VU.
          If RANGE='A', then WL and WU are the global Gerschgorin bounds
                        on the spectrum.
          If RANGE='I', then WL and WU are computed by DLAEBZ from the
                        index range specified.
[out]IBLOCK
          IBLOCK is INTEGER array, dimension (N)
          At each row/column j where E(j) is zero or small, the
          matrix T is considered to split into a block diagonal
          matrix.  On exit, if INFO = 0, IBLOCK(i) specifies to which
          block (from 1 to the number of blocks) the eigenvalue W(i)
          belongs.  (DLARRD may use the remaining N-M elements as
          workspace.)
[out]INDEXW
          INDEXW is INTEGER array, dimension (N)
          The indices of the eigenvalues within each block (submatrix);
          for example, INDEXW(i)= j and IBLOCK(i)=k imply that the
          i-th eigenvalue W(i) is the j-th eigenvalue in block k.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (4*N)
[out]IWORK
          IWORK is INTEGER 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:  some or all of the eigenvalues failed to converge or
                were not computed:
                =1 or 3: Bisection failed to converge for some
                        eigenvalues; these eigenvalues are flagged by a
                        negative block number.  The effect is that the
                        eigenvalues may not be as accurate as the
                        absolute and relative tolerances.  This is
                        generally caused by unexpectedly inaccurate
                        arithmetic.
                =2 or 3: RANGE='I' only: Not all of the eigenvalues
                        IL:IU were found.
                        Effect: M < IU+1-IL
                        Cause:  non-monotonic arithmetic, causing the
                                Sturm sequence to be non-monotonic.
                        Cure:   recalculate, using RANGE='A', and pick
                                out eigenvalues IL:IU.  In some cases,
                                increasing the PARAMETER "FUDGE" may
                                make things work.
                = 4:    RANGE='I', and the Gershgorin interval
                        initially used was too small.  No eigenvalues
                        were computed.
                        Probable cause: your machine has sloppy
                                        floating-point arithmetic.
                        Cure: Increase the PARAMETER "FUDGE",
                              recompile, and try again.
Internal Parameters:
  FUDGE   DOUBLE PRECISION, default = 2
          A "fudge factor" to widen the Gershgorin intervals.  Ideally,
          a value of 1 should work, but on machines with sloppy
          arithmetic, this needs to be larger.  The default for
          publicly released versions should be large enough to handle
          the worst machine around.  Note that this has no effect
          on accuracy of the solution.
Contributors:
W. Kahan, University of California, Berkeley, USA
Beresford Parlett, University of California, Berkeley, USA
Jim Demmel, University of California, Berkeley, USA
Inderjit Dhillon, University of Texas, Austin, USA
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 323 of file dlarrd.f.

323 *
324 * -- LAPACK auxiliary routine (version 3.4.2) --
325 * -- LAPACK is a software package provided by Univ. of Tennessee, --
326 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
327 * September 2012
328 *
329 * .. Scalar Arguments ..
330  CHARACTER order, range
331  INTEGER il, info, iu, m, n, nsplit
332  DOUBLE PRECISION pivmin, reltol, vl, vu, wl, wu
333 * ..
334 * .. Array Arguments ..
335  INTEGER iblock( * ), indexw( * ),
336  $ isplit( * ), iwork( * )
337  DOUBLE PRECISION d( * ), e( * ), e2( * ),
338  $ gers( * ), w( * ), werr( * ), work( * )
339 * ..
340 *
341 * =====================================================================
342 *
343 * .. Parameters ..
344  DOUBLE PRECISION zero, one, two, half, fudge
345  parameter( zero = 0.0d0, one = 1.0d0,
346  $ two = 2.0d0, half = one/two,
347  $ fudge = two )
348  INTEGER allrng, valrng, indrng
349  parameter( allrng = 1, valrng = 2, indrng = 3 )
350 * ..
351 * .. Local Scalars ..
352  LOGICAL ncnvrg, toofew
353  INTEGER i, ib, ibegin, idiscl, idiscu, ie, iend, iinfo,
354  $ im, in, ioff, iout, irange, itmax, itmp1,
355  $ itmp2, iw, iwoff, j, jblk, jdisc, je, jee, nb,
356  $ nwl, nwu
357  DOUBLE PRECISION atoli, eps, gl, gu, rtoli, tmp1, tmp2,
358  $ tnorm, uflow, wkill, wlu, wul
359 
360 * ..
361 * .. Local Arrays ..
362  INTEGER idumma( 1 )
363 * ..
364 * .. External Functions ..
365  LOGICAL lsame
366  INTEGER ilaenv
367  DOUBLE PRECISION dlamch
368  EXTERNAL lsame, ilaenv, dlamch
369 * ..
370 * .. External Subroutines ..
371  EXTERNAL dlaebz
372 * ..
373 * .. Intrinsic Functions ..
374  INTRINSIC abs, int, log, max, min
375 * ..
376 * .. Executable Statements ..
377 *
378  info = 0
379 *
380 * Decode RANGE
381 *
382  IF( lsame( range, 'A' ) ) THEN
383  irange = allrng
384  ELSE IF( lsame( range, 'V' ) ) THEN
385  irange = valrng
386  ELSE IF( lsame( range, 'I' ) ) THEN
387  irange = indrng
388  ELSE
389  irange = 0
390  END IF
391 *
392 * Check for Errors
393 *
394  IF( irange.LE.0 ) THEN
395  info = -1
396  ELSE IF( .NOT.(lsame(order,'B').OR.lsame(order,'E')) ) THEN
397  info = -2
398  ELSE IF( n.LT.0 ) THEN
399  info = -3
400  ELSE IF( irange.EQ.valrng ) THEN
401  IF( vl.GE.vu )
402  $ info = -5
403  ELSE IF( irange.EQ.indrng .AND.
404  $ ( il.LT.1 .OR. il.GT.max( 1, n ) ) ) THEN
405  info = -6
406  ELSE IF( irange.EQ.indrng .AND.
407  $ ( iu.LT.min( n, il ) .OR. iu.GT.n ) ) THEN
408  info = -7
409  END IF
410 *
411  IF( info.NE.0 ) THEN
412  RETURN
413  END IF
414 
415 * Initialize error flags
416  info = 0
417  ncnvrg = .false.
418  toofew = .false.
419 
420 * Quick return if possible
421  m = 0
422  IF( n.EQ.0 ) RETURN
423 
424 * Simplification:
425  IF( irange.EQ.indrng .AND. il.EQ.1 .AND. iu.EQ.n ) irange = 1
426 
427 * Get machine constants
428  eps = dlamch( 'P' )
429  uflow = dlamch( 'U' )
430 
431 
432 * Special Case when N=1
433 * Treat case of 1x1 matrix for quick return
434  IF( n.EQ.1 ) THEN
435  IF( (irange.EQ.allrng).OR.
436  $ ((irange.EQ.valrng).AND.(d(1).GT.vl).AND.(d(1).LE.vu)).OR.
437  $ ((irange.EQ.indrng).AND.(il.EQ.1).AND.(iu.EQ.1)) ) THEN
438  m = 1
439  w(1) = d(1)
440 * The computation error of the eigenvalue is zero
441  werr(1) = zero
442  iblock( 1 ) = 1
443  indexw( 1 ) = 1
444  ENDIF
445  RETURN
446  END IF
447 
448 * NB is the minimum vector length for vector bisection, or 0
449 * if only scalar is to be done.
450  nb = ilaenv( 1, 'DSTEBZ', ' ', n, -1, -1, -1 )
451  IF( nb.LE.1 ) nb = 0
452 
453 * Find global spectral radius
454  gl = d(1)
455  gu = d(1)
456  DO 5 i = 1,n
457  gl = min( gl, gers( 2*i - 1))
458  gu = max( gu, gers(2*i) )
459  5 CONTINUE
460 * Compute global Gerschgorin bounds and spectral diameter
461  tnorm = max( abs( gl ), abs( gu ) )
462  gl = gl - fudge*tnorm*eps*n - fudge*two*pivmin
463  gu = gu + fudge*tnorm*eps*n + fudge*two*pivmin
464 * [JAN/28/2009] remove the line below since SPDIAM variable not use
465 * SPDIAM = GU - GL
466 * Input arguments for DLAEBZ:
467 * The relative tolerance. An interval (a,b] lies within
468 * "relative tolerance" if b-a < RELTOL*max(|a|,|b|),
469  rtoli = reltol
470 * Set the absolute tolerance for interval convergence to zero to force
471 * interval convergence based on relative size of the interval.
472 * This is dangerous because intervals might not converge when RELTOL is
473 * small. But at least a very small number should be selected so that for
474 * strongly graded matrices, the code can get relatively accurate
475 * eigenvalues.
476  atoli = fudge*two*uflow + fudge*two*pivmin
477 
478  IF( irange.EQ.indrng ) THEN
479 
480 * RANGE='I': Compute an interval containing eigenvalues
481 * IL through IU. The initial interval [GL,GU] from the global
482 * Gerschgorin bounds GL and GU is refined by DLAEBZ.
483  itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /
484  $ log( two ) ) + 2
485  work( n+1 ) = gl
486  work( n+2 ) = gl
487  work( n+3 ) = gu
488  work( n+4 ) = gu
489  work( n+5 ) = gl
490  work( n+6 ) = gu
491  iwork( 1 ) = -1
492  iwork( 2 ) = -1
493  iwork( 3 ) = n + 1
494  iwork( 4 ) = n + 1
495  iwork( 5 ) = il - 1
496  iwork( 6 ) = iu
497 *
498  CALL dlaebz( 3, itmax, n, 2, 2, nb, atoli, rtoli, pivmin,
499  $ d, e, e2, iwork( 5 ), work( n+1 ), work( n+5 ), iout,
500  $ iwork, w, iblock, iinfo )
501  IF( iinfo .NE. 0 ) THEN
502  info = iinfo
503  RETURN
504  END IF
505 * On exit, output intervals may not be ordered by ascending negcount
506  IF( iwork( 6 ).EQ.iu ) THEN
507  wl = work( n+1 )
508  wlu = work( n+3 )
509  nwl = iwork( 1 )
510  wu = work( n+4 )
511  wul = work( n+2 )
512  nwu = iwork( 4 )
513  ELSE
514  wl = work( n+2 )
515  wlu = work( n+4 )
516  nwl = iwork( 2 )
517  wu = work( n+3 )
518  wul = work( n+1 )
519  nwu = iwork( 3 )
520  END IF
521 * On exit, the interval [WL, WLU] contains a value with negcount NWL,
522 * and [WUL, WU] contains a value with negcount NWU.
523  IF( nwl.LT.0 .OR. nwl.GE.n .OR. nwu.LT.1 .OR. nwu.GT.n ) THEN
524  info = 4
525  RETURN
526  END IF
527 
528  ELSEIF( irange.EQ.valrng ) THEN
529  wl = vl
530  wu = vu
531 
532  ELSEIF( irange.EQ.allrng ) THEN
533  wl = gl
534  wu = gu
535  ENDIF
536 
537 
538 
539 * Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU.
540 * NWL accumulates the number of eigenvalues .le. WL,
541 * NWU accumulates the number of eigenvalues .le. WU
542  m = 0
543  iend = 0
544  info = 0
545  nwl = 0
546  nwu = 0
547 *
548  DO 70 jblk = 1, nsplit
549  ioff = iend
550  ibegin = ioff + 1
551  iend = isplit( jblk )
552  in = iend - ioff
553 *
554  IF( in.EQ.1 ) THEN
555 * 1x1 block
556  IF( wl.GE.d( ibegin )-pivmin )
557  $ nwl = nwl + 1
558  IF( wu.GE.d( ibegin )-pivmin )
559  $ nwu = nwu + 1
560  IF( irange.EQ.allrng .OR.
561  $ ( wl.LT.d( ibegin )-pivmin
562  $ .AND. wu.GE. d( ibegin )-pivmin ) ) THEN
563  m = m + 1
564  w( m ) = d( ibegin )
565  werr(m) = zero
566 * The gap for a single block doesn't matter for the later
567 * algorithm and is assigned an arbitrary large value
568  iblock( m ) = jblk
569  indexw( m ) = 1
570  END IF
571 
572 * Disabled 2x2 case because of a failure on the following matrix
573 * RANGE = 'I', IL = IU = 4
574 * Original Tridiagonal, d = [
575 * -0.150102010615740E+00
576 * -0.849897989384260E+00
577 * -0.128208148052635E-15
578 * 0.128257718286320E-15
579 * ];
580 * e = [
581 * -0.357171383266986E+00
582 * -0.180411241501588E-15
583 * -0.175152352710251E-15
584 * ];
585 *
586 * ELSE IF( IN.EQ.2 ) THEN
587 ** 2x2 block
588 * DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 )
589 * TMP1 = HALF*(D(IBEGIN)+D(IEND))
590 * L1 = TMP1 - DISC
591 * IF( WL.GE. L1-PIVMIN )
592 * $ NWL = NWL + 1
593 * IF( WU.GE. L1-PIVMIN )
594 * $ NWU = NWU + 1
595 * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE.
596 * $ L1-PIVMIN ) ) THEN
597 * M = M + 1
598 * W( M ) = L1
599 ** The uncertainty of eigenvalues of a 2x2 matrix is very small
600 * WERR( M ) = EPS * ABS( W( M ) ) * TWO
601 * IBLOCK( M ) = JBLK
602 * INDEXW( M ) = 1
603 * ENDIF
604 * L2 = TMP1 + DISC
605 * IF( WL.GE. L2-PIVMIN )
606 * $ NWL = NWL + 1
607 * IF( WU.GE. L2-PIVMIN )
608 * $ NWU = NWU + 1
609 * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE.
610 * $ L2-PIVMIN ) ) THEN
611 * M = M + 1
612 * W( M ) = L2
613 ** The uncertainty of eigenvalues of a 2x2 matrix is very small
614 * WERR( M ) = EPS * ABS( W( M ) ) * TWO
615 * IBLOCK( M ) = JBLK
616 * INDEXW( M ) = 2
617 * ENDIF
618  ELSE
619 * General Case - block of size IN >= 2
620 * Compute local Gerschgorin interval and use it as the initial
621 * interval for DLAEBZ
622  gu = d( ibegin )
623  gl = d( ibegin )
624  tmp1 = zero
625 
626  DO 40 j = ibegin, iend
627  gl = min( gl, gers( 2*j - 1))
628  gu = max( gu, gers(2*j) )
629  40 CONTINUE
630 * [JAN/28/2009]
631 * change SPDIAM by TNORM in lines 2 and 3 thereafter
632 * line 1: remove computation of SPDIAM (not useful anymore)
633 * SPDIAM = GU - GL
634 * GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN
635 * GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN
636  gl = gl - fudge*tnorm*eps*in - fudge*pivmin
637  gu = gu + fudge*tnorm*eps*in + fudge*pivmin
638 *
639  IF( irange.GT.1 ) THEN
640  IF( gu.LT.wl ) THEN
641 * the local block contains none of the wanted eigenvalues
642  nwl = nwl + in
643  nwu = nwu + in
644  GO TO 70
645  END IF
646 * refine search interval if possible, only range (WL,WU] matters
647  gl = max( gl, wl )
648  gu = min( gu, wu )
649  IF( gl.GE.gu )
650  $ GO TO 70
651  END IF
652 
653 * Find negcount of initial interval boundaries GL and GU
654  work( n+1 ) = gl
655  work( n+in+1 ) = gu
656  CALL dlaebz( 1, 0, in, in, 1, nb, atoli, rtoli, pivmin,
657  $ d( ibegin ), e( ibegin ), e2( ibegin ),
658  $ idumma, work( n+1 ), work( n+2*in+1 ), im,
659  $ iwork, w( m+1 ), iblock( m+1 ), iinfo )
660  IF( iinfo .NE. 0 ) THEN
661  info = iinfo
662  RETURN
663  END IF
664 *
665  nwl = nwl + iwork( 1 )
666  nwu = nwu + iwork( in+1 )
667  iwoff = m - iwork( 1 )
668 
669 * Compute Eigenvalues
670  itmax = int( ( log( gu-gl+pivmin )-log( pivmin ) ) /
671  $ log( two ) ) + 2
672  CALL dlaebz( 2, itmax, in, in, 1, nb, atoli, rtoli, pivmin,
673  $ d( ibegin ), e( ibegin ), e2( ibegin ),
674  $ idumma, work( n+1 ), work( n+2*in+1 ), iout,
675  $ iwork, w( m+1 ), iblock( m+1 ), iinfo )
676  IF( iinfo .NE. 0 ) THEN
677  info = iinfo
678  RETURN
679  END IF
680 *
681 * Copy eigenvalues into W and IBLOCK
682 * Use -JBLK for block number for unconverged eigenvalues.
683 * Loop over the number of output intervals from DLAEBZ
684  DO 60 j = 1, iout
685 * eigenvalue approximation is middle point of interval
686  tmp1 = half*( work( j+n )+work( j+in+n ) )
687 * semi length of error interval
688  tmp2 = half*abs( work( j+n )-work( j+in+n ) )
689  IF( j.GT.iout-iinfo ) THEN
690 * Flag non-convergence.
691  ncnvrg = .true.
692  ib = -jblk
693  ELSE
694  ib = jblk
695  END IF
696  DO 50 je = iwork( j ) + 1 + iwoff,
697  $ iwork( j+in ) + iwoff
698  w( je ) = tmp1
699  werr( je ) = tmp2
700  indexw( je ) = je - iwoff
701  iblock( je ) = ib
702  50 CONTINUE
703  60 CONTINUE
704 *
705  m = m + im
706  END IF
707  70 CONTINUE
708 
709 * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
710 * If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
711  IF( irange.EQ.indrng ) THEN
712  idiscl = il - 1 - nwl
713  idiscu = nwu - iu
714 *
715  IF( idiscl.GT.0 ) THEN
716  im = 0
717  DO 80 je = 1, m
718 * Remove some of the smallest eigenvalues from the left so that
719 * at the end IDISCL =0. Move all eigenvalues up to the left.
720  IF( w( je ).LE.wlu .AND. idiscl.GT.0 ) THEN
721  idiscl = idiscl - 1
722  ELSE
723  im = im + 1
724  w( im ) = w( je )
725  werr( im ) = werr( je )
726  indexw( im ) = indexw( je )
727  iblock( im ) = iblock( je )
728  END IF
729  80 CONTINUE
730  m = im
731  END IF
732  IF( idiscu.GT.0 ) THEN
733 * Remove some of the largest eigenvalues from the right so that
734 * at the end IDISCU =0. Move all eigenvalues up to the left.
735  im=m+1
736  DO 81 je = m, 1, -1
737  IF( w( je ).GE.wul .AND. idiscu.GT.0 ) THEN
738  idiscu = idiscu - 1
739  ELSE
740  im = im - 1
741  w( im ) = w( je )
742  werr( im ) = werr( je )
743  indexw( im ) = indexw( je )
744  iblock( im ) = iblock( je )
745  END IF
746  81 CONTINUE
747  jee = 0
748  DO 82 je = im, m
749  jee = jee + 1
750  w( jee ) = w( je )
751  werr( jee ) = werr( je )
752  indexw( jee ) = indexw( je )
753  iblock( jee ) = iblock( je )
754  82 CONTINUE
755  m = m-im+1
756  END IF
757 
758  IF( idiscl.GT.0 .OR. idiscu.GT.0 ) THEN
759 * Code to deal with effects of bad arithmetic. (If N(w) is
760 * monotone non-decreasing, this should never happen.)
761 * Some low eigenvalues to be discarded are not in (WL,WLU],
762 * or high eigenvalues to be discarded are not in (WUL,WU]
763 * so just kill off the smallest IDISCL/largest IDISCU
764 * eigenvalues, by marking the corresponding IBLOCK = 0
765  IF( idiscl.GT.0 ) THEN
766  wkill = wu
767  DO 100 jdisc = 1, idiscl
768  iw = 0
769  DO 90 je = 1, m
770  IF( iblock( je ).NE.0 .AND.
771  $ ( w( je ).LT.wkill .OR. iw.EQ.0 ) ) THEN
772  iw = je
773  wkill = w( je )
774  END IF
775  90 CONTINUE
776  iblock( iw ) = 0
777  100 CONTINUE
778  END IF
779  IF( idiscu.GT.0 ) THEN
780  wkill = wl
781  DO 120 jdisc = 1, idiscu
782  iw = 0
783  DO 110 je = 1, m
784  IF( iblock( je ).NE.0 .AND.
785  $ ( w( je ).GE.wkill .OR. iw.EQ.0 ) ) THEN
786  iw = je
787  wkill = w( je )
788  END IF
789  110 CONTINUE
790  iblock( iw ) = 0
791  120 CONTINUE
792  END IF
793 * Now erase all eigenvalues with IBLOCK set to zero
794  im = 0
795  DO 130 je = 1, m
796  IF( iblock( je ).NE.0 ) THEN
797  im = im + 1
798  w( im ) = w( je )
799  werr( im ) = werr( je )
800  indexw( im ) = indexw( je )
801  iblock( im ) = iblock( je )
802  END IF
803  130 CONTINUE
804  m = im
805  END IF
806  IF( idiscl.LT.0 .OR. idiscu.LT.0 ) THEN
807  toofew = .true.
808  END IF
809  END IF
810 *
811  IF(( irange.EQ.allrng .AND. m.NE.n ).OR.
812  $ ( irange.EQ.indrng .AND. m.NE.iu-il+1 ) ) THEN
813  toofew = .true.
814  END IF
815 
816 * If ORDER='B', do nothing the eigenvalues are already sorted by
817 * block.
818 * If ORDER='E', sort the eigenvalues from smallest to largest
819 
820  IF( lsame(order,'E') .AND. nsplit.GT.1 ) THEN
821  DO 150 je = 1, m - 1
822  ie = 0
823  tmp1 = w( je )
824  DO 140 j = je + 1, m
825  IF( w( j ).LT.tmp1 ) THEN
826  ie = j
827  tmp1 = w( j )
828  END IF
829  140 CONTINUE
830  IF( ie.NE.0 ) THEN
831  tmp2 = werr( ie )
832  itmp1 = iblock( ie )
833  itmp2 = indexw( ie )
834  w( ie ) = w( je )
835  werr( ie ) = werr( je )
836  iblock( ie ) = iblock( je )
837  indexw( ie ) = indexw( je )
838  w( je ) = tmp1
839  werr( je ) = tmp2
840  iblock( je ) = itmp1
841  indexw( je ) = itmp2
842  END IF
843  150 CONTINUE
844  END IF
845 *
846  info = 0
847  IF( ncnvrg )
848  $ info = info + 1
849  IF( toofew )
850  $ info = info + 2
851  RETURN
852 *
853 * End of DLARRD
854 *
subroutine dlaebz(IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, NAB, WORK, IWORK, INFO)
DLAEBZ computes the number of eigenvalues of a real symmetric tridiagonal matrix which are less than ...
Definition: dlaebz.f:321
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dlarre ( character  RANGE,
integer  N,
double precision  VL,
double precision  VU,
integer  IL,
integer  IU,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( * )  E2,
double precision  RTOL1,
double precision  RTOL2,
double precision  SPLTOL,
integer  NSPLIT,
integer, dimension( * )  ISPLIT,
integer  M,
double precision, dimension( * )  W,
double precision, dimension( * )  WERR,
double precision, dimension( * )  WGAP,
integer, dimension( * )  IBLOCK,
integer, dimension( * )  INDEXW,
double precision, dimension( * )  GERS,
double precision  PIVMIN,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

DLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduced block Ti, finds base representations and eigenvalues.

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

Purpose:
 To find the desired eigenvalues of a given real symmetric
 tridiagonal matrix T, DLARRE sets any "small" off-diagonal
 elements to zero, and for each unreduced block T_i, it finds
 (a) a suitable shift at one end of the block's spectrum,
 (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
 (c) eigenvalues of each L_i D_i L_i^T.
 The representations and eigenvalues found are then used by
 DSTEMR to compute the eigenvectors of T.
 The accuracy varies depending on whether bisection is used to
 find a few eigenvalues or the dqds algorithm (subroutine DLASQ2) to
 conpute all and then discard any unwanted one.
 As an added benefit, DLARRE also outputs the n
 Gerschgorin intervals for the matrices L_i D_i L_i^T.
Parameters
[in]RANGE
          RANGE is CHARACTER*1
          = 'A': ("All")   all eigenvalues will be found.
          = 'V': ("Value") all eigenvalues in the half-open interval
                           (VL, VU] will be found.
          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
                           entire matrix) will be found.
[in]N
          N is INTEGER
          The order of the matrix. N > 0.
[in,out]VL
          VL is DOUBLE PRECISION
[in,out]VU
          VU is DOUBLE PRECISION
          If RANGE='V', the lower and upper bounds for the eigenvalues.
          Eigenvalues less than or equal to VL, or greater than VU,
          will not be returned.  VL < VU.
          If RANGE='I' or ='A', DLARRE computes bounds on the desired
          part of the spectrum.
[in]IL
          IL is INTEGER
[in]IU
          IU is INTEGER
          If RANGE='I', the indices (in ascending order) of the
          smallest and largest eigenvalues to be returned.
          1 <= IL <= IU <= N.
[in,out]D
          D is DOUBLE PRECISION array, dimension (N)
          On entry, the N diagonal elements of the tridiagonal
          matrix T.
          On exit, the N diagonal elements of the diagonal
          matrices D_i.
[in,out]E
          E is DOUBLE PRECISION array, dimension (N)
          On entry, the first (N-1) entries contain the subdiagonal
          elements of the tridiagonal matrix T; E(N) need not be set.
          On exit, E contains the subdiagonal elements of the unit
          bidiagonal matrices L_i. The entries E( ISPLIT( I ) ),
          1 <= I <= NSPLIT, contain the base points sigma_i on output.
[in,out]E2
          E2 is DOUBLE PRECISION array, dimension (N)
          On entry, the first (N-1) entries contain the SQUARES of the
          subdiagonal elements of the tridiagonal matrix T;
          E2(N) need not be set.
          On exit, the entries E2( ISPLIT( I ) ),
          1 <= I <= NSPLIT, have been set to zero
[in]RTOL1
          RTOL1 is DOUBLE PRECISION
[in]RTOL2
          RTOL2 is DOUBLE PRECISION
           Parameters for bisection.
           An interval [LEFT,RIGHT] has converged if
           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
[in]SPLTOL
          SPLTOL is DOUBLE PRECISION
          The threshold for splitting.
[out]NSPLIT
          NSPLIT is INTEGER
          The number of blocks T splits into. 1 <= NSPLIT <= N.
[out]ISPLIT
          ISPLIT is INTEGER array, dimension (N)
          The splitting points, at which T breaks up into blocks.
          The first block consists of rows/columns 1 to ISPLIT(1),
          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
          etc., and the NSPLIT-th consists of rows/columns
          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
[out]M
          M is INTEGER
          The total number of eigenvalues (of all L_i D_i L_i^T)
          found.
[out]W
          W is DOUBLE PRECISION array, dimension (N)
          The first M elements contain the eigenvalues. The
          eigenvalues of each of the blocks, L_i D_i L_i^T, are
          sorted in ascending order ( DLARRE may use the
          remaining N-M elements as workspace).
[out]WERR
          WERR is DOUBLE PRECISION array, dimension (N)
          The error bound on the corresponding eigenvalue in W.
[out]WGAP
          WGAP is DOUBLE PRECISION array, dimension (N)
          The separation from the right neighbor eigenvalue in W.
          The gap is only with respect to the eigenvalues of the same block
          as each block has its own representation tree.
          Exception: at the right end of a block we store the left gap
[out]IBLOCK
          IBLOCK is INTEGER array, dimension (N)
          The indices of the blocks (submatrices) associated with the
          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
          W(i) belongs to the first block from the top, =2 if W(i)
          belongs to the second block, etc.
[out]INDEXW
          INDEXW is INTEGER array, dimension (N)
          The indices of the eigenvalues within each block (submatrix);
          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
          i-th eigenvalue W(i) is the 10-th eigenvalue in block 2
[out]GERS
          GERS is DOUBLE PRECISION array, dimension (2*N)
          The N Gerschgorin intervals (the i-th Gerschgorin interval
          is (GERS(2*i-1), GERS(2*i)).
[out]PIVMIN
          PIVMIN is DOUBLE PRECISION
          The minimum pivot in the Sturm sequence for T.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (6*N)
          Workspace.
[out]IWORK
          IWORK is INTEGER array, dimension (5*N)
          Workspace.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          > 0:  A problem occured in DLARRE.
          < 0:  One of the called subroutines signaled an internal problem.
                Needs inspection of the corresponding parameter IINFO
                for further information.

          =-1:  Problem in DLARRD.
          = 2:  No base representation could be found in MAXTRY iterations.
                Increasing MAXTRY and recompilation might be a remedy.
          =-3:  Problem in DLARRB when computing the refined root
                representation for DLASQ2.
          =-4:  Problem in DLARRB when preforming bisection on the
                desired part of the spectrum.
          =-5:  Problem in DLASQ2.
          =-6:  Problem in DLASQ2.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
  The base representations are required to suffer very little
  element growth and consequently define all their eigenvalues to
  high relative accuracy.
Contributors:
Beresford Parlett, University of California, Berkeley, USA
Jim Demmel, University of California, Berkeley, USA
Inderjit Dhillon, University of Texas, Austin, USA
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA

Definition at line 299 of file dlarre.f.

299 *
300 * -- LAPACK auxiliary routine (version 3.4.2) --
301 * -- LAPACK is a software package provided by Univ. of Tennessee, --
302 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
303 * September 2012
304 *
305 * .. Scalar Arguments ..
306  CHARACTER range
307  INTEGER il, info, iu, m, n, nsplit
308  DOUBLE PRECISION pivmin, rtol1, rtol2, spltol, vl, vu
309 * ..
310 * .. Array Arguments ..
311  INTEGER iblock( * ), isplit( * ), iwork( * ),
312  $ indexw( * )
313  DOUBLE PRECISION d( * ), e( * ), e2( * ), gers( * ),
314  $ w( * ),werr( * ), wgap( * ), work( * )
315 * ..
316 *
317 * =====================================================================
318 *
319 * .. Parameters ..
320  DOUBLE PRECISION fac, four, fourth, fudge, half, hndrd,
321  $ maxgrowth, one, pert, two, zero
322  parameter( zero = 0.0d0, one = 1.0d0,
323  $ two = 2.0d0, four=4.0d0,
324  $ hndrd = 100.0d0,
325  $ pert = 8.0d0,
326  $ half = one/two, fourth = one/four, fac= half,
327  $ maxgrowth = 64.0d0, fudge = 2.0d0 )
328  INTEGER maxtry, allrng, indrng, valrng
329  parameter( maxtry = 6, allrng = 1, indrng = 2,
330  $ valrng = 3 )
331 * ..
332 * .. Local Scalars ..
333  LOGICAL forceb, norep, usedqd
334  INTEGER cnt, cnt1, cnt2, i, ibegin, idum, iend, iinfo,
335  $ in, indl, indu, irange, j, jblk, mb, mm,
336  $ wbegin, wend
337  DOUBLE PRECISION avgap, bsrtol, clwdth, dmax, dpivot, eabs,
338  $ emax, eold, eps, gl, gu, isleft, isrght, rtl,
339  $ rtol, s1, s2, safmin, sgndef, sigma, spdiam,
340  $ tau, tmp, tmp1
341 
342 
343 * ..
344 * .. Local Arrays ..
345  INTEGER iseed( 4 )
346 * ..
347 * .. External Functions ..
348  LOGICAL lsame
349  DOUBLE PRECISION dlamch
350  EXTERNAL dlamch, lsame
351 
352 * ..
353 * .. External Subroutines ..
354  EXTERNAL dcopy, dlarnv, dlarra, dlarrb, dlarrc, dlarrd,
355  $ dlasq2
356 * ..
357 * .. Intrinsic Functions ..
358  INTRINSIC abs, max, min
359 
360 * ..
361 * .. Executable Statements ..
362 *
363 
364  info = 0
365 
366 *
367 * Decode RANGE
368 *
369  IF( lsame( range, 'A' ) ) THEN
370  irange = allrng
371  ELSE IF( lsame( range, 'V' ) ) THEN
372  irange = valrng
373  ELSE IF( lsame( range, 'I' ) ) THEN
374  irange = indrng
375  END IF
376 
377  m = 0
378 
379 * Get machine constants
380  safmin = dlamch( 'S' )
381  eps = dlamch( 'P' )
382 
383 * Set parameters
384  rtl = sqrt(eps)
385  bsrtol = sqrt(eps)
386 
387 * Treat case of 1x1 matrix for quick return
388  IF( n.EQ.1 ) THEN
389  IF( (irange.EQ.allrng).OR.
390  $ ((irange.EQ.valrng).AND.(d(1).GT.vl).AND.(d(1).LE.vu)).OR.
391  $ ((irange.EQ.indrng).AND.(il.EQ.1).AND.(iu.EQ.1)) ) THEN
392  m = 1
393  w(1) = d(1)
394 * The computation error of the eigenvalue is zero
395  werr(1) = zero
396  wgap(1) = zero
397  iblock( 1 ) = 1
398  indexw( 1 ) = 1
399  gers(1) = d( 1 )
400  gers(2) = d( 1 )
401  ENDIF
402 * store the shift for the initial RRR, which is zero in this case
403  e(1) = zero
404  RETURN
405  END IF
406 
407 * General case: tridiagonal matrix of order > 1
408 *
409 * Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter.
410 * Compute maximum off-diagonal entry and pivmin.
411  gl = d(1)
412  gu = d(1)
413  eold = zero
414  emax = zero
415  e(n) = zero
416  DO 5 i = 1,n
417  werr(i) = zero
418  wgap(i) = zero
419  eabs = abs( e(i) )
420  IF( eabs .GE. emax ) THEN
421  emax = eabs
422  END IF
423  tmp1 = eabs + eold
424  gers( 2*i-1) = d(i) - tmp1
425  gl = min( gl, gers( 2*i - 1))
426  gers( 2*i ) = d(i) + tmp1
427  gu = max( gu, gers(2*i) )
428  eold = eabs
429  5 CONTINUE
430 * The minimum pivot allowed in the Sturm sequence for T
431  pivmin = safmin * max( one, emax**2 )
432 * Compute spectral diameter. The Gerschgorin bounds give an
433 * estimate that is wrong by at most a factor of SQRT(2)
434  spdiam = gu - gl
435 
436 * Compute splitting points
437  CALL dlarra( n, d, e, e2, spltol, spdiam,
438  $ nsplit, isplit, iinfo )
439 
440 * Can force use of bisection instead of faster DQDS.
441 * Option left in the code for future multisection work.
442  forceb = .false.
443 
444 * Initialize USEDQD, DQDS should be used for ALLRNG unless someone
445 * explicitly wants bisection.
446  usedqd = (( irange.EQ.allrng ) .AND. (.NOT.forceb))
447 
448  IF( (irange.EQ.allrng) .AND. (.NOT. forceb) ) THEN
449 * Set interval [VL,VU] that contains all eigenvalues
450  vl = gl
451  vu = gu
452  ELSE
453 * We call DLARRD to find crude approximations to the eigenvalues
454 * in the desired range. In case IRANGE = INDRNG, we also obtain the
455 * interval (VL,VU] that contains all the wanted eigenvalues.
456 * An interval [LEFT,RIGHT] has converged if
457 * RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT))
458 * DLARRD needs a WORK of size 4*N, IWORK of size 3*N
459  CALL dlarrd( range, 'B', n, vl, vu, il, iu, gers,
460  $ bsrtol, d, e, e2, pivmin, nsplit, isplit,
461  $ mm, w, werr, vl, vu, iblock, indexw,
462  $ work, iwork, iinfo )
463  IF( iinfo.NE.0 ) THEN
464  info = -1
465  RETURN
466  ENDIF
467 * Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0
468  DO 14 i = mm+1,n
469  w( i ) = zero
470  werr( i ) = zero
471  iblock( i ) = 0
472  indexw( i ) = 0
473  14 CONTINUE
474  END IF
475 
476 
477 ***
478 * Loop over unreduced blocks
479  ibegin = 1
480  wbegin = 1
481  DO 170 jblk = 1, nsplit
482  iend = isplit( jblk )
483  in = iend - ibegin + 1
484 
485 * 1 X 1 block
486  IF( in.EQ.1 ) THEN
487  IF( (irange.EQ.allrng).OR.( (irange.EQ.valrng).AND.
488  $ ( d( ibegin ).GT.vl ).AND.( d( ibegin ).LE.vu ) )
489  $ .OR. ( (irange.EQ.indrng).AND.(iblock(wbegin).EQ.jblk))
490  $ ) THEN
491  m = m + 1
492  w( m ) = d( ibegin )
493  werr(m) = zero
494 * The gap for a single block doesn't matter for the later
495 * algorithm and is assigned an arbitrary large value
496  wgap(m) = zero
497  iblock( m ) = jblk
498  indexw( m ) = 1
499  wbegin = wbegin + 1
500  ENDIF
501 * E( IEND ) holds the shift for the initial RRR
502  e( iend ) = zero
503  ibegin = iend + 1
504  GO TO 170
505  END IF
506 *
507 * Blocks of size larger than 1x1
508 *
509 * E( IEND ) will hold the shift for the initial RRR, for now set it =0
510  e( iend ) = zero
511 *
512 * Find local outer bounds GL,GU for the block
513  gl = d(ibegin)
514  gu = d(ibegin)
515  DO 15 i = ibegin , iend
516  gl = min( gers( 2*i-1 ), gl )
517  gu = max( gers( 2*i ), gu )
518  15 CONTINUE
519  spdiam = gu - gl
520 
521  IF(.NOT. ((irange.EQ.allrng).AND.(.NOT.forceb)) ) THEN
522 * Count the number of eigenvalues in the current block.
523  mb = 0
524  DO 20 i = wbegin,mm
525  IF( iblock(i).EQ.jblk ) THEN
526  mb = mb+1
527  ELSE
528  GOTO 21
529  ENDIF
530  20 CONTINUE
531  21 CONTINUE
532 
533  IF( mb.EQ.0) THEN
534 * No eigenvalue in the current block lies in the desired range
535 * E( IEND ) holds the shift for the initial RRR
536  e( iend ) = zero
537  ibegin = iend + 1
538  GO TO 170
539  ELSE
540 
541 * Decide whether dqds or bisection is more efficient
542  usedqd = ( (mb .GT. fac*in) .AND. (.NOT.forceb) )
543  wend = wbegin + mb - 1
544 * Calculate gaps for the current block
545 * In later stages, when representations for individual
546 * eigenvalues are different, we use SIGMA = E( IEND ).
547  sigma = zero
548  DO 30 i = wbegin, wend - 1
549  wgap( i ) = max( zero,
550  $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
551  30 CONTINUE
552  wgap( wend ) = max( zero,
553  $ vu - sigma - (w( wend )+werr( wend )))
554 * Find local index of the first and last desired evalue.
555  indl = indexw(wbegin)
556  indu = indexw( wend )
557  ENDIF
558  ENDIF
559  IF(( (irange.EQ.allrng) .AND. (.NOT. forceb) ).OR.usedqd) THEN
560 * Case of DQDS
561 * Find approximations to the extremal eigenvalues of the block
562  CALL dlarrk( in, 1, gl, gu, d(ibegin),
563  $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
564  IF( iinfo.NE.0 ) THEN
565  info = -1
566  RETURN
567  ENDIF
568  isleft = max(gl, tmp - tmp1
569  $ - hndrd * eps* abs(tmp - tmp1))
570 
571  CALL dlarrk( in, in, gl, gu, d(ibegin),
572  $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
573  IF( iinfo.NE.0 ) THEN
574  info = -1
575  RETURN
576  ENDIF
577  isrght = min(gu, tmp + tmp1
578  $ + hndrd * eps * abs(tmp + tmp1))
579 * Improve the estimate of the spectral diameter
580  spdiam = isrght - isleft
581  ELSE
582 * Case of bisection
583 * Find approximations to the wanted extremal eigenvalues
584  isleft = max(gl, w(wbegin) - werr(wbegin)
585  $ - hndrd * eps*abs(w(wbegin)- werr(wbegin) ))
586  isrght = min(gu,w(wend) + werr(wend)
587  $ + hndrd * eps * abs(w(wend)+ werr(wend)))
588  ENDIF
589 
590 
591 * Decide whether the base representation for the current block
592 * L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I
593 * should be on the left or the right end of the current block.
594 * The strategy is to shift to the end which is "more populated"
595 * Furthermore, decide whether to use DQDS for the computation of
596 * the eigenvalue approximations at the end of DLARRE or bisection.
597 * dqds is chosen if all eigenvalues are desired or the number of
598 * eigenvalues to be computed is large compared to the blocksize.
599  IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
600 * If all the eigenvalues have to be computed, we use dqd
601  usedqd = .true.
602 * INDL is the local index of the first eigenvalue to compute
603  indl = 1
604  indu = in
605 * MB = number of eigenvalues to compute
606  mb = in
607  wend = wbegin + mb - 1
608 * Define 1/4 and 3/4 points of the spectrum
609  s1 = isleft + fourth * spdiam
610  s2 = isrght - fourth * spdiam
611  ELSE
612 * DLARRD has computed IBLOCK and INDEXW for each eigenvalue
613 * approximation.
614 * choose sigma
615  IF( usedqd ) THEN
616  s1 = isleft + fourth * spdiam
617  s2 = isrght - fourth * spdiam
618  ELSE
619  tmp = min(isrght,vu) - max(isleft,vl)
620  s1 = max(isleft,vl) + fourth * tmp
621  s2 = min(isrght,vu) - fourth * tmp
622  ENDIF
623  ENDIF
624 
625 * Compute the negcount at the 1/4 and 3/4 points
626  IF(mb.GT.1) THEN
627  CALL dlarrc( 'T', in, s1, s2, d(ibegin),
628  $ e(ibegin), pivmin, cnt, cnt1, cnt2, iinfo)
629  ENDIF
630 
631  IF(mb.EQ.1) THEN
632  sigma = gl
633  sgndef = one
634  ELSEIF( cnt1 - indl .GE. indu - cnt2 ) THEN
635  IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
636  sigma = max(isleft,gl)
637  ELSEIF( usedqd ) THEN
638 * use Gerschgorin bound as shift to get pos def matrix
639 * for dqds
640  sigma = isleft
641  ELSE
642 * use approximation of the first desired eigenvalue of the
643 * block as shift
644  sigma = max(isleft,vl)
645  ENDIF
646  sgndef = one
647  ELSE
648  IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
649  sigma = min(isrght,gu)
650  ELSEIF( usedqd ) THEN
651 * use Gerschgorin bound as shift to get neg def matrix
652 * for dqds
653  sigma = isrght
654  ELSE
655 * use approximation of the first desired eigenvalue of the
656 * block as shift
657  sigma = min(isrght,vu)
658  ENDIF
659  sgndef = -one
660  ENDIF
661 
662 
663 * An initial SIGMA has been chosen that will be used for computing
664 * T - SIGMA I = L D L^T
665 * Define the increment TAU of the shift in case the initial shift
666 * needs to be refined to obtain a factorization with not too much
667 * element growth.
668  IF( usedqd ) THEN
669 * The initial SIGMA was to the outer end of the spectrum
670 * the matrix is definite and we need not retreat.
671  tau = spdiam*eps*n + two*pivmin
672  tau = max( tau,two*eps*abs(sigma) )
673  ELSE
674  IF(mb.GT.1) THEN
675  clwdth = w(wend) + werr(wend) - w(wbegin) - werr(wbegin)
676  avgap = abs(clwdth / dble(wend-wbegin))
677  IF( sgndef.EQ.one ) THEN
678  tau = half*max(wgap(wbegin),avgap)
679  tau = max(tau,werr(wbegin))
680  ELSE
681  tau = half*max(wgap(wend-1),avgap)
682  tau = max(tau,werr(wend))
683  ENDIF
684  ELSE
685  tau = werr(wbegin)
686  ENDIF
687  ENDIF
688 *
689  DO 80 idum = 1, maxtry
690 * Compute L D L^T factorization of tridiagonal matrix T - sigma I.
691 * Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
692 * pivots in WORK(2*IN+1:3*IN)
693  dpivot = d( ibegin ) - sigma
694  work( 1 ) = dpivot
695  dmax = abs( work(1) )
696  j = ibegin
697  DO 70 i = 1, in - 1
698  work( 2*in+i ) = one / work( i )
699  tmp = e( j )*work( 2*in+i )
700  work( in+i ) = tmp
701  dpivot = ( d( j+1 )-sigma ) - tmp*e( j )
702  work( i+1 ) = dpivot
703  dmax = max( dmax, abs(dpivot) )
704  j = j + 1
705  70 CONTINUE
706 * check for element growth
707  IF( dmax .GT. maxgrowth*spdiam ) THEN
708  norep = .true.
709  ELSE
710  norep = .false.
711  ENDIF
712  IF( usedqd .AND. .NOT.norep ) THEN
713 * Ensure the definiteness of the representation
714 * All entries of D (of L D L^T) must have the same sign
715  DO 71 i = 1, in
716  tmp = sgndef*work( i )
717  IF( tmp.LT.zero ) norep = .true.
718  71 CONTINUE
719  ENDIF
720  IF(norep) THEN
721 * Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin
722 * shift which makes the matrix definite. So we should end up
723 * here really only in the case of IRANGE = VALRNG or INDRNG.
724  IF( idum.EQ.maxtry-1 ) THEN
725  IF( sgndef.EQ.one ) THEN
726 * The fudged Gerschgorin shift should succeed
727  sigma =
728  $ gl - fudge*spdiam*eps*n - fudge*two*pivmin
729  ELSE
730  sigma =
731  $ gu + fudge*spdiam*eps*n + fudge*two*pivmin
732  END IF
733  ELSE
734  sigma = sigma - sgndef * tau
735  tau = two * tau
736  END IF
737  ELSE
738 * an initial RRR is found
739  GO TO 83
740  END IF
741  80 CONTINUE
742 * if the program reaches this point, no base representation could be
743 * found in MAXTRY iterations.
744  info = 2
745  RETURN
746 
747  83 CONTINUE
748 * At this point, we have found an initial base representation
749 * T - SIGMA I = L D L^T with not too much element growth.
750 * Store the shift.
751  e( iend ) = sigma
752 * Store D and L.
753  CALL dcopy( in, work, 1, d( ibegin ), 1 )
754  CALL dcopy( in-1, work( in+1 ), 1, e( ibegin ), 1 )
755 
756 
757  IF(mb.GT.1 ) THEN
758 *
759 * Perturb each entry of the base representation by a small
760 * (but random) relative amount to overcome difficulties with
761 * glued matrices.
762 *
763  DO 122 i = 1, 4
764  iseed( i ) = 1
765  122 CONTINUE
766 
767  CALL dlarnv(2, iseed, 2*in-1, work(1))
768  DO 125 i = 1,in-1
769  d(ibegin+i-1) = d(ibegin+i-1)*(one+eps*pert*work(i))
770  e(ibegin+i-1) = e(ibegin+i-1)*(one+eps*pert*work(in+i))
771  125 CONTINUE
772  d(iend) = d(iend)*(one+eps*four*work(in))
773 *
774  ENDIF
775 *
776 * Don't update the Gerschgorin intervals because keeping track
777 * of the updates would be too much work in DLARRV.
778 * We update W instead and use it to locate the proper Gerschgorin
779 * intervals.
780 
781 * Compute the required eigenvalues of L D L' by bisection or dqds
782  IF ( .NOT.usedqd ) THEN
783 * If DLARRD has been used, shift the eigenvalue approximations
784 * according to their representation. This is necessary for
785 * a uniform DLARRV since dqds computes eigenvalues of the
786 * shifted representation. In DLARRV, W will always hold the
787 * UNshifted eigenvalue approximation.
788  DO 134 j=wbegin,wend
789  w(j) = w(j) - sigma
790  werr(j) = werr(j) + abs(w(j)) * eps
791  134 CONTINUE
792 * call DLARRB to reduce eigenvalue error of the approximations
793 * from DLARRD
794  DO 135 i = ibegin, iend-1
795  work( i ) = d( i ) * e( i )**2
796  135 CONTINUE
797 * use bisection to find EV from INDL to INDU
798  CALL dlarrb(in, d(ibegin), work(ibegin),
799  $ indl, indu, rtol1, rtol2, indl-1,
800  $ w(wbegin), wgap(wbegin), werr(wbegin),
801  $ work( 2*n+1 ), iwork, pivmin, spdiam,
802  $ in, iinfo )
803  IF( iinfo .NE. 0 ) THEN
804  info = -4
805  RETURN
806  END IF
807 * DLARRB computes all gaps correctly except for the last one
808 * Record distance to VU/GU
809  wgap( wend ) = max( zero,
810  $ ( vu-sigma ) - ( w( wend ) + werr( wend ) ) )
811  DO 138 i = indl, indu
812  m = m + 1
813  iblock(m) = jblk
814  indexw(m) = i
815  138 CONTINUE
816  ELSE
817 * Call dqds to get all eigs (and then possibly delete unwanted
818 * eigenvalues).
819 * Note that dqds finds the eigenvalues of the L D L^T representation
820 * of T to high relative accuracy. High relative accuracy
821 * might be lost when the shift of the RRR is subtracted to obtain
822 * the eigenvalues of T. However, T is not guaranteed to define its
823 * eigenvalues to high relative accuracy anyway.
824 * Set RTOL to the order of the tolerance used in DLASQ2
825 * This is an ESTIMATED error, the worst case bound is 4*N*EPS
826 * which is usually too large and requires unnecessary work to be
827 * done by bisection when computing the eigenvectors
828  rtol = log(dble(in)) * four * eps
829  j = ibegin
830  DO 140 i = 1, in - 1
831  work( 2*i-1 ) = abs( d( j ) )
832  work( 2*i ) = e( j )*e( j )*work( 2*i-1 )
833  j = j + 1
834  140 CONTINUE
835  work( 2*in-1 ) = abs( d( iend ) )
836  work( 2*in ) = zero
837  CALL dlasq2( in, work, iinfo )
838  IF( iinfo .NE. 0 ) THEN
839 * If IINFO = -5 then an index is part of a tight cluster
840 * and should be changed. The index is in IWORK(1) and the
841 * gap is in WORK(N+1)
842  info = -5
843  RETURN
844  ELSE
845 * Test that all eigenvalues are positive as expected
846  DO 149 i = 1, in
847  IF( work( i ).LT.zero ) THEN
848  info = -6
849  RETURN
850  ENDIF
851  149 CONTINUE
852  END IF
853  IF( sgndef.GT.zero ) THEN
854  DO 150 i = indl, indu
855  m = m + 1
856  w( m ) = work( in-i+1 )
857  iblock( m ) = jblk
858  indexw( m ) = i
859  150 CONTINUE
860  ELSE
861  DO 160 i = indl, indu
862  m = m + 1
863  w( m ) = -work( i )
864  iblock( m ) = jblk
865  indexw( m ) = i
866  160 CONTINUE
867  END IF
868 
869  DO 165 i = m - mb + 1, m
870 * the value of RTOL below should be the tolerance in DLASQ2
871  werr( i ) = rtol * abs( w(i) )
872  165 CONTINUE
873  DO 166 i = m - mb + 1, m - 1
874 * compute the right gap between the intervals
875  wgap( i ) = max( zero,
876  $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
877  166 CONTINUE
878  wgap( m ) = max( zero,
879  $ ( vu-sigma ) - ( w( m ) + werr( m ) ) )
880  END IF
881 * proceed with next block
882  ibegin = iend + 1
883  wbegin = wend + 1
884  170 CONTINUE
885 *
886 
887  RETURN
888 *
889 * end of DLARRE
890 *
subroutine dlarrc(JOBT, N, VL, VU, D, E, PIVMIN, EIGCNT, LCNT, RCNT, INFO)
DLARRC computes the number of eigenvalues of the symmetric tridiagonal matrix.
Definition: dlarrc.f:138
subroutine dlarrd(RANGE, ORDER, N, VL, VU, IL, IU, GERS, RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, M, W, WERR, WL, WU, IBLOCK, INDEXW, WORK, IWORK, INFO)
DLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy.
Definition: dlarrd.f:323
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dlarrb(N, D, LLD, IFIRST, ILAST, RTOL1, RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK, PIVMIN, SPDIAM, TWIST, INFO)
DLARRB provides limited bisection to locate eigenvalues for more accuracy.
Definition: dlarrb.f:198
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dlarra(N, D, E, E2, SPLTOL, TNRM, NSPLIT, ISPLIT, INFO)
DLARRA computes the splitting points with the specified threshold.
Definition: dlarra.f:138
subroutine dlarrk(N, IW, GL, GU, D, E2, PIVMIN, RELTOL, W, WERR, INFO)
DLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy.
Definition: dlarrk.f:147
subroutine dlasq2(N, Z, INFO)
DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated ...
Definition: dlasq2.f:114
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: dlarnv.f:99
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dlarrf ( integer  N,
double precision, dimension( * )  D,
double precision, dimension( * )  L,
double precision, dimension( * )  LD,
integer  CLSTRT,
integer  CLEND,
double precision, dimension( * )  W,
double precision, dimension( * )  WGAP,
double precision, dimension( * )  WERR,
double precision  SPDIAM,
double precision  CLGAPL,
double precision  CLGAPR,
double precision  PIVMIN,
double precision  SIGMA,
double precision, dimension( * )  DPLUS,
double precision, dimension( * )  LPLUS,
double precision, dimension( * )  WORK,
integer  INFO 
)

DLARRF finds a new relatively robust representation such that at least one of the eigenvalues is relatively isolated.

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

Purpose:
 Given the initial representation L D L^T and its cluster of close
 eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ...
 W( CLEND ), DLARRF finds a new relatively robust representation
 L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the
 eigenvalues of L(+) D(+) L(+)^T is relatively isolated.
Parameters
[in]N
          N is INTEGER
          The order of the matrix (subblock, if the matrix splitted).
[in]D
          D is DOUBLE PRECISION array, dimension (N)
          The N diagonal elements of the diagonal matrix D.
[in]L
          L is DOUBLE PRECISION array, dimension (N-1)
          The (N-1) subdiagonal elements of the unit bidiagonal
          matrix L.
[in]LD
          LD is DOUBLE PRECISION array, dimension (N-1)
          The (N-1) elements L(i)*D(i).
[in]CLSTRT
          CLSTRT is INTEGER
          The index of the first eigenvalue in the cluster.
[in]CLEND
          CLEND is INTEGER
          The index of the last eigenvalue in the cluster.
[in]W
          W is DOUBLE PRECISION array, dimension
          dimension is >=  (CLEND-CLSTRT+1)
          The eigenvalue APPROXIMATIONS of L D L^T in ascending order.
          W( CLSTRT ) through W( CLEND ) form the cluster of relatively
          close eigenalues.
[in,out]WGAP
          WGAP is DOUBLE PRECISION array, dimension
          dimension is >=  (CLEND-CLSTRT+1)
          The separation from the right neighbor eigenvalue in W.
[in]WERR
          WERR is DOUBLE PRECISION array, dimension
          dimension is  >=  (CLEND-CLSTRT+1)
          WERR contain the semiwidth of the uncertainty
          interval of the corresponding eigenvalue APPROXIMATION in W
[in]SPDIAM
          SPDIAM is DOUBLE PRECISION
          estimate of the spectral diameter obtained from the
          Gerschgorin intervals
[in]CLGAPL
          CLGAPL is DOUBLE PRECISION
[in]CLGAPR
          CLGAPR is DOUBLE PRECISION
          absolute gap on each end of the cluster.
          Set by the calling routine to protect against shifts too close
          to eigenvalues outside the cluster.
[in]PIVMIN
          PIVMIN is DOUBLE PRECISION
          The minimum pivot allowed in the Sturm sequence.
[out]SIGMA
          SIGMA is DOUBLE PRECISION
          The shift used to form L(+) D(+) L(+)^T.
[out]DPLUS
          DPLUS is DOUBLE PRECISION array, dimension (N)
          The N diagonal elements of the diagonal matrix D(+).
[out]LPLUS
          LPLUS is DOUBLE PRECISION array, dimension (N-1)
          The first (N-1) elements of LPLUS contain the subdiagonal
          elements of the unit bidiagonal matrix L(+).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (2*N)
          Workspace.
[out]INFO
          INFO is INTEGER
          Signals processing OK (=0) or failure (=1)
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015
Contributors:
Beresford Parlett, University of California, Berkeley, USA
Jim Demmel, University of California, Berkeley, USA
Inderjit Dhillon, University of Texas, Austin, USA
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA

Definition at line 195 of file dlarrf.f.

195 *
196 * -- LAPACK auxiliary routine (version 3.6.0) --
197 * -- LAPACK is a software package provided by Univ. of Tennessee, --
198 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
199 * November 2015
200 *
201 * .. Scalar Arguments ..
202  INTEGER clstrt, clend, info, n
203  DOUBLE PRECISION clgapl, clgapr, pivmin, sigma, spdiam
204 * ..
205 * .. Array Arguments ..
206  DOUBLE PRECISION d( * ), dplus( * ), l( * ), ld( * ),
207  $ lplus( * ), w( * ), wgap( * ), werr( * ), work( * )
208 * ..
209 *
210 * =====================================================================
211 *
212 * .. Parameters ..
213  DOUBLE PRECISION four, maxgrowth1, maxgrowth2, one, quart, two
214  parameter( one = 1.0d0, two = 2.0d0, four = 4.0d0,
215  $ quart = 0.25d0,
216  $ maxgrowth1 = 8.d0,
217  $ maxgrowth2 = 8.d0 )
218 * ..
219 * .. Local Scalars ..
220  LOGICAL dorrr1, forcer, nofail, sawnan1, sawnan2, tryrrr1
221  INTEGER i, indx, ktry, ktrymax, sleft, sright, shift
222  parameter( ktrymax = 1, sleft = 1, sright = 2 )
223  DOUBLE PRECISION avgap, bestshift, clwdth, eps, fact, fail,
224  $ fail2, growthbound, ldelta, ldmax, lsigma,
225  $ max1, max2, mingap, oldp, prod, rdelta, rdmax,
226  $ rrr1, rrr2, rsigma, s, smlgrowth, tmp, znm2
227 * ..
228 * .. External Functions ..
229  LOGICAL disnan
230  DOUBLE PRECISION dlamch
231  EXTERNAL disnan, dlamch
232 * ..
233 * .. External Subroutines ..
234  EXTERNAL dcopy
235 * ..
236 * .. Intrinsic Functions ..
237  INTRINSIC abs
238 * ..
239 * .. Executable Statements ..
240 *
241  info = 0
242  fact = dble(2**ktrymax)
243  eps = dlamch( 'Precision' )
244  shift = 0
245  forcer = .false.
246 
247 
248 * Note that we cannot guarantee that for any of the shifts tried,
249 * the factorization has a small or even moderate element growth.
250 * There could be Ritz values at both ends of the cluster and despite
251 * backing off, there are examples where all factorizations tried
252 * (in IEEE mode, allowing zero pivots & infinities) have INFINITE
253 * element growth.
254 * For this reason, we should use PIVMIN in this subroutine so that at
255 * least the L D L^T factorization exists. It can be checked afterwards
256 * whether the element growth caused bad residuals/orthogonality.
257 
258 * Decide whether the code should accept the best among all
259 * representations despite large element growth or signal INFO=1
260 * Setting NOFAIL to .FALSE. for quick fix for bug 113
261  nofail = .false.
262 *
263 
264 * Compute the average gap length of the cluster
265  clwdth = abs(w(clend)-w(clstrt)) + werr(clend) + werr(clstrt)
266  avgap = clwdth / dble(clend-clstrt)
267  mingap = min(clgapl, clgapr)
268 * Initial values for shifts to both ends of cluster
269  lsigma = min(w( clstrt ),w( clend )) - werr( clstrt )
270  rsigma = max(w( clstrt ),w( clend )) + werr( clend )
271 
272 * Use a small fudge to make sure that we really shift to the outside
273  lsigma = lsigma - abs(lsigma)* four * eps
274  rsigma = rsigma + abs(rsigma)* four * eps
275 
276 * Compute upper bounds for how much to back off the initial shifts
277  ldmax = quart * mingap + two * pivmin
278  rdmax = quart * mingap + two * pivmin
279 
280  ldelta = max(avgap,wgap( clstrt ))/fact
281  rdelta = max(avgap,wgap( clend-1 ))/fact
282 *
283 * Initialize the record of the best representation found
284 *
285  s = dlamch( 'S' )
286  smlgrowth = one / s
287  fail = dble(n-1)*mingap/(spdiam*eps)
288  fail2 = dble(n-1)*mingap/(spdiam*sqrt(eps))
289  bestshift = lsigma
290 *
291 * while (KTRY <= KTRYMAX)
292  ktry = 0
293  growthbound = maxgrowth1*spdiam
294 
295  5 CONTINUE
296  sawnan1 = .false.
297  sawnan2 = .false.
298 * Ensure that we do not back off too much of the initial shifts
299  ldelta = min(ldmax,ldelta)
300  rdelta = min(rdmax,rdelta)
301 
302 * Compute the element growth when shifting to both ends of the cluster
303 * accept the shift if there is no element growth at one of the two ends
304 
305 * Left end
306  s = -lsigma
307  dplus( 1 ) = d( 1 ) + s
308  IF(abs(dplus(1)).LT.pivmin) THEN
309  dplus(1) = -pivmin
310 * Need to set SAWNAN1 because refined RRR test should not be used
311 * in this case
312  sawnan1 = .true.
313  ENDIF
314  max1 = abs( dplus( 1 ) )
315  DO 6 i = 1, n - 1
316  lplus( i ) = ld( i ) / dplus( i )
317  s = s*lplus( i )*l( i ) - lsigma
318  dplus( i+1 ) = d( i+1 ) + s
319  IF(abs(dplus(i+1)).LT.pivmin) THEN
320  dplus(i+1) = -pivmin
321 * Need to set SAWNAN1 because refined RRR test should not be used
322 * in this case
323  sawnan1 = .true.
324  ENDIF
325  max1 = max( max1,abs(dplus(i+1)) )
326  6 CONTINUE
327  sawnan1 = sawnan1 .OR. disnan( max1 )
328 
329  IF( forcer .OR.
330  $ (max1.LE.growthbound .AND. .NOT.sawnan1 ) ) THEN
331  sigma = lsigma
332  shift = sleft
333  GOTO 100
334  ENDIF
335 
336 * Right end
337  s = -rsigma
338  work( 1 ) = d( 1 ) + s
339  IF(abs(work(1)).LT.pivmin) THEN
340  work(1) = -pivmin
341 * Need to set SAWNAN2 because refined RRR test should not be used
342 * in this case
343  sawnan2 = .true.
344  ENDIF
345  max2 = abs( work( 1 ) )
346  DO 7 i = 1, n - 1
347  work( n+i ) = ld( i ) / work( i )
348  s = s*work( n+i )*l( i ) - rsigma
349  work( i+1 ) = d( i+1 ) + s
350  IF(abs(work(i+1)).LT.pivmin) THEN
351  work(i+1) = -pivmin
352 * Need to set SAWNAN2 because refined RRR test should not be used
353 * in this case
354  sawnan2 = .true.
355  ENDIF
356  max2 = max( max2,abs(work(i+1)) )
357  7 CONTINUE
358  sawnan2 = sawnan2 .OR. disnan( max2 )
359 
360  IF( forcer .OR.
361  $ (max2.LE.growthbound .AND. .NOT.sawnan2 ) ) THEN
362  sigma = rsigma
363  shift = sright
364  GOTO 100
365  ENDIF
366 * If we are at this point, both shifts led to too much element growth
367 
368 * Record the better of the two shifts (provided it didn't lead to NaN)
369  IF(sawnan1.AND.sawnan2) THEN
370 * both MAX1 and MAX2 are NaN
371  GOTO 50
372  ELSE
373  IF( .NOT.sawnan1 ) THEN
374  indx = 1
375  IF(max1.LE.smlgrowth) THEN
376  smlgrowth = max1
377  bestshift = lsigma
378  ENDIF
379  ENDIF
380  IF( .NOT.sawnan2 ) THEN
381  IF(sawnan1 .OR. max2.LE.max1) indx = 2
382  IF(max2.LE.smlgrowth) THEN
383  smlgrowth = max2
384  bestshift = rsigma
385  ENDIF
386  ENDIF
387  ENDIF
388 
389 * If we are here, both the left and the right shift led to
390 * element growth. If the element growth is moderate, then
391 * we may still accept the representation, if it passes a
392 * refined test for RRR. This test supposes that no NaN occurred.
393 * Moreover, we use the refined RRR test only for isolated clusters.
394  IF((clwdth.LT.mingap/dble(128)) .AND.
395  $ (min(max1,max2).LT.fail2)
396  $ .AND.(.NOT.sawnan1).AND.(.NOT.sawnan2)) THEN
397  dorrr1 = .true.
398  ELSE
399  dorrr1 = .false.
400  ENDIF
401  tryrrr1 = .true.
402  IF( tryrrr1 .AND. dorrr1 ) THEN
403  IF(indx.EQ.1) THEN
404  tmp = abs( dplus( n ) )
405  znm2 = one
406  prod = one
407  oldp = one
408  DO 15 i = n-1, 1, -1
409  IF( prod .LE. eps ) THEN
410  prod =
411  $ ((dplus(i+1)*work(n+i+1))/(dplus(i)*work(n+i)))*oldp
412  ELSE
413  prod = prod*abs(work(n+i))
414  END IF
415  oldp = prod
416  znm2 = znm2 + prod**2
417  tmp = max( tmp, abs( dplus( i ) * prod ))
418  15 CONTINUE
419  rrr1 = tmp/( spdiam * sqrt( znm2 ) )
420  IF (rrr1.LE.maxgrowth2) THEN
421  sigma = lsigma
422  shift = sleft
423  GOTO 100
424  ENDIF
425  ELSE IF(indx.EQ.2) THEN
426  tmp = abs( work( n ) )
427  znm2 = one
428  prod = one
429  oldp = one
430  DO 16 i = n-1, 1, -1
431  IF( prod .LE. eps ) THEN
432  prod = ((work(i+1)*lplus(i+1))/(work(i)*lplus(i)))*oldp
433  ELSE
434  prod = prod*abs(lplus(i))
435  END IF
436  oldp = prod
437  znm2 = znm2 + prod**2
438  tmp = max( tmp, abs( work( i ) * prod ))
439  16 CONTINUE
440  rrr2 = tmp/( spdiam * sqrt( znm2 ) )
441  IF (rrr2.LE.maxgrowth2) THEN
442  sigma = rsigma
443  shift = sright
444  GOTO 100
445  ENDIF
446  END IF
447  ENDIF
448 
449  50 CONTINUE
450 
451  IF (ktry.LT.ktrymax) THEN
452 * If we are here, both shifts failed also the RRR test.
453 * Back off to the outside
454  lsigma = max( lsigma - ldelta,
455  $ lsigma - ldmax)
456  rsigma = min( rsigma + rdelta,
457  $ rsigma + rdmax )
458  ldelta = two * ldelta
459  rdelta = two * rdelta
460  ktry = ktry + 1
461  GOTO 5
462  ELSE
463 * None of the representations investigated satisfied our
464 * criteria. Take the best one we found.
465  IF((smlgrowth.LT.fail).OR.nofail) THEN
466  lsigma = bestshift
467  rsigma = bestshift
468  forcer = .true.
469  GOTO 5
470  ELSE
471  info = 1
472  RETURN
473  ENDIF
474  END IF
475 
476  100 CONTINUE
477  IF (shift.EQ.sleft) THEN
478  ELSEIF (shift.EQ.sright) THEN
479 * store new L and D back into DPLUS, LPLUS
480  CALL dcopy( n, work, 1, dplus, 1 )
481  CALL dcopy( n-1, work(n+1), 1, lplus, 1 )
482  ENDIF
483 
484  RETURN
485 *
486 * End of DLARRF
487 *
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:61
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dlarrj ( integer  N,
double precision, dimension( * )  D,
double precision, dimension( * )  E2,
integer  IFIRST,
integer  ILAST,
double precision  RTOL,
integer  OFFSET,
double precision, dimension( * )  W,
double precision, dimension( * )  WERR,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
double precision  PIVMIN,
double precision  SPDIAM,
integer  INFO 
)

DLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T.

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

Purpose:
 Given the initial eigenvalue approximations of T, DLARRJ
 does  bisection to refine the eigenvalues of T,
 W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
 guesses for these eigenvalues are input in W, the corresponding estimate
 of the error in these guesses in WERR. During bisection, intervals
 [left, right] are maintained by storing their mid-points and
 semi-widths in the arrays W and WERR respectively.
Parameters
[in]N
          N is INTEGER
          The order of the matrix.
[in]D
          D is DOUBLE PRECISION array, dimension (N)
          The N diagonal elements of T.
[in]E2
          E2 is DOUBLE PRECISION array, dimension (N-1)
          The Squares of the (N-1) subdiagonal elements of T.
[in]IFIRST
          IFIRST is INTEGER
          The index of the first eigenvalue to be computed.
[in]ILAST
          ILAST is INTEGER
          The index of the last eigenvalue to be computed.
[in]RTOL
          RTOL is DOUBLE PRECISION
          Tolerance for the convergence of the bisection intervals.
          An interval [LEFT,RIGHT] has converged if
          RIGHT-LEFT.LT.RTOL*MAX(|LEFT|,|RIGHT|).
[in]OFFSET
          OFFSET is INTEGER
          Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET
          through ILAST-OFFSET elements of these arrays are to be used.
[in,out]W
          W is DOUBLE PRECISION array, dimension (N)
          On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
          estimates of the eigenvalues of L D L^T indexed IFIRST through
          ILAST.
          On output, these estimates are refined.
[in,out]WERR
          WERR is DOUBLE PRECISION array, dimension (N)
          On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
          the errors in the estimates of the corresponding elements in W.
          On output, these errors are refined.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (2*N)
          Workspace.
[out]IWORK
          IWORK is INTEGER array, dimension (2*N)
          Workspace.
[in]PIVMIN
          PIVMIN is DOUBLE PRECISION
          The minimum pivot in the Sturm sequence for T.
[in]SPDIAM
          SPDIAM is DOUBLE PRECISION
          The spectral diameter of T.
[out]INFO
          INFO is INTEGER
          Error flag.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Contributors:
Beresford Parlett, University of California, Berkeley, USA
Jim Demmel, University of California, Berkeley, USA
Inderjit Dhillon, University of Texas, Austin, USA
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA

Definition at line 170 of file dlarrj.f.

170 *
171 * -- LAPACK auxiliary routine (version 3.4.2) --
172 * -- LAPACK is a software package provided by Univ. of Tennessee, --
173 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
174 * September 2012
175 *
176 * .. Scalar Arguments ..
177  INTEGER ifirst, ilast, info, n, offset
178  DOUBLE PRECISION pivmin, rtol, spdiam
179 * ..
180 * .. Array Arguments ..
181  INTEGER iwork( * )
182  DOUBLE PRECISION d( * ), e2( * ), w( * ),
183  $ werr( * ), work( * )
184 * ..
185 *
186 * =====================================================================
187 *
188 * .. Parameters ..
189  DOUBLE PRECISION zero, one, two, half
190  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
191  $ half = 0.5d0 )
192  INTEGER maxitr
193 * ..
194 * .. Local Scalars ..
195  INTEGER cnt, i, i1, i2, ii, iter, j, k, next, nint,
196  $ olnint, p, prev, savi1
197  DOUBLE PRECISION dplus, fac, left, mid, right, s, tmp, width
198 *
199 * ..
200 * .. Intrinsic Functions ..
201  INTRINSIC abs, max
202 * ..
203 * .. Executable Statements ..
204 *
205  info = 0
206 *
207  maxitr = int( ( log( spdiam+pivmin )-log( pivmin ) ) /
208  $ log( two ) ) + 2
209 *
210 * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
211 * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
212 * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
213 * for an unconverged interval is set to the index of the next unconverged
214 * interval, and is -1 or 0 for a converged interval. Thus a linked
215 * list of unconverged intervals is set up.
216 *
217 
218  i1 = ifirst
219  i2 = ilast
220 * The number of unconverged intervals
221  nint = 0
222 * The last unconverged interval found
223  prev = 0
224  DO 75 i = i1, i2
225  k = 2*i
226  ii = i - offset
227  left = w( ii ) - werr( ii )
228  mid = w(ii)
229  right = w( ii ) + werr( ii )
230  width = right - mid
231  tmp = max( abs( left ), abs( right ) )
232 
233 * The following test prevents the test of converged intervals
234  IF( width.LT.rtol*tmp ) THEN
235 * This interval has already converged and does not need refinement.
236 * (Note that the gaps might change through refining the
237 * eigenvalues, however, they can only get bigger.)
238 * Remove it from the list.
239  iwork( k-1 ) = -1
240 * Make sure that I1 always points to the first unconverged interval
241  IF((i.EQ.i1).AND.(i.LT.i2)) i1 = i + 1
242  IF((prev.GE.i1).AND.(i.LE.i2)) iwork( 2*prev-1 ) = i + 1
243  ELSE
244 * unconverged interval found
245  prev = i
246 * Make sure that [LEFT,RIGHT] contains the desired eigenvalue
247 *
248 * Do while( CNT(LEFT).GT.I-1 )
249 *
250  fac = one
251  20 CONTINUE
252  cnt = 0
253  s = left
254  dplus = d( 1 ) - s
255  IF( dplus.LT.zero ) cnt = cnt + 1
256  DO 30 j = 2, n
257  dplus = d( j ) - s - e2( j-1 )/dplus
258  IF( dplus.LT.zero ) cnt = cnt + 1
259  30 CONTINUE
260  IF( cnt.GT.i-1 ) THEN
261  left = left - werr( ii )*fac
262  fac = two*fac
263  GO TO 20
264  END IF
265 *
266 * Do while( CNT(RIGHT).LT.I )
267 *
268  fac = one
269  50 CONTINUE
270  cnt = 0
271  s = right
272  dplus = d( 1 ) - s
273  IF( dplus.LT.zero ) cnt = cnt + 1
274  DO 60 j = 2, n
275  dplus = d( j ) - s - e2( j-1 )/dplus
276  IF( dplus.LT.zero ) cnt = cnt + 1
277  60 CONTINUE
278  IF( cnt.LT.i ) THEN
279  right = right + werr( ii )*fac
280  fac = two*fac
281  GO TO 50
282  END IF
283  nint = nint + 1
284  iwork( k-1 ) = i + 1
285  iwork( k ) = cnt
286  END IF
287  work( k-1 ) = left
288  work( k ) = right
289  75 CONTINUE
290 
291 
292  savi1 = i1
293 *
294 * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
295 * and while (ITER.LT.MAXITR)
296 *
297  iter = 0
298  80 CONTINUE
299  prev = i1 - 1
300  i = i1
301  olnint = nint
302 
303  DO 100 p = 1, olnint
304  k = 2*i
305  ii = i - offset
306  next = iwork( k-1 )
307  left = work( k-1 )
308  right = work( k )
309  mid = half*( left + right )
310 
311 * semiwidth of interval
312  width = right - mid
313  tmp = max( abs( left ), abs( right ) )
314 
315  IF( ( width.LT.rtol*tmp ) .OR.
316  $ (iter.EQ.maxitr) )THEN
317 * reduce number of unconverged intervals
318  nint = nint - 1
319 * Mark interval as converged.
320  iwork( k-1 ) = 0
321  IF( i1.EQ.i ) THEN
322  i1 = next
323  ELSE
324 * Prev holds the last unconverged interval previously examined
325  IF(prev.GE.i1) iwork( 2*prev-1 ) = next
326  END IF
327  i = next
328  GO TO 100
329  END IF
330  prev = i
331 *
332 * Perform one bisection step
333 *
334  cnt = 0
335  s = mid
336  dplus = d( 1 ) - s
337  IF( dplus.LT.zero ) cnt = cnt + 1
338  DO 90 j = 2, n
339  dplus = d( j ) - s - e2( j-1 )/dplus
340  IF( dplus.LT.zero ) cnt = cnt + 1
341  90 CONTINUE
342  IF( cnt.LE.i-1 ) THEN
343  work( k-1 ) = mid
344  ELSE
345  work( k ) = mid
346  END IF
347  i = next
348 
349  100 CONTINUE
350  iter = iter + 1
351 * do another loop if there are still unconverged intervals
352 * However, in the last iteration, all intervals are accepted
353 * since this is the best we can do.
354  IF( ( nint.GT.0 ).AND.(iter.LE.maxitr) ) GO TO 80
355 *
356 *
357 * At this point, all the intervals have converged
358  DO 110 i = savi1, ilast
359  k = 2*i
360  ii = i - offset
361 * All intervals marked by '0' have been refined.
362  IF( iwork( k-1 ).EQ.0 ) THEN
363  w( ii ) = half*( work( k-1 )+work( k ) )
364  werr( ii ) = work( k ) - w( ii )
365  END IF
366  110 CONTINUE
367 *
368 
369  RETURN
370 *
371 * End of DLARRJ
372 *

Here is the caller graph for this function:

subroutine dlarrk ( integer  N,
integer  IW,
double precision  GL,
double precision  GU,
double precision, dimension( * )  D,
double precision, dimension( * )  E2,
double precision  PIVMIN,
double precision  RELTOL,
double precision  W,
double precision  WERR,
integer  INFO 
)

DLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy.

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

Purpose:
 DLARRK computes one eigenvalue of a symmetric tridiagonal
 matrix T to suitable accuracy. This is an auxiliary code to be
 called from DSTEMR.

 To avoid overflow, the matrix must be scaled so that its
 largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
 accuracy, it should not be much smaller than that.

 See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
 Matrix", Report CS41, Computer Science Dept., Stanford
 University, July 21, 1966.
Parameters
[in]N
          N is INTEGER
          The order of the tridiagonal matrix T.  N >= 0.
[in]IW
          IW is INTEGER
          The index of the eigenvalues to be returned.
[in]GL
          GL is DOUBLE PRECISION
[in]GU
          GU is DOUBLE PRECISION
          An upper and a lower bound on the eigenvalue.
[in]D
          D is DOUBLE PRECISION array, dimension (N)
          The n diagonal elements of the tridiagonal matrix T.
[in]E2
          E2 is DOUBLE PRECISION array, dimension (N-1)
          The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
[in]PIVMIN
          PIVMIN is DOUBLE PRECISION
          The minimum pivot allowed in the Sturm sequence for T.
[in]RELTOL
          RELTOL is DOUBLE PRECISION
          The minimum relative width of an interval.  When an interval
          is narrower than RELTOL times the larger (in
          magnitude) endpoint, then it is considered to be
          sufficiently small, i.e., converged.  Note: this should
          always be at least radix*machine epsilon.
[out]W
          W is DOUBLE PRECISION
[out]WERR
          WERR is DOUBLE PRECISION
          The error bound on the corresponding eigenvalue approximation
          in W.
[out]INFO
          INFO is INTEGER
          = 0:       Eigenvalue converged
          = -1:      Eigenvalue did NOT converge
Internal Parameters:
  FUDGE   DOUBLE PRECISION, default = 2
          A "fudge factor" to widen the Gershgorin intervals.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 147 of file dlarrk.f.

147 *
148 * -- LAPACK auxiliary routine (version 3.4.2) --
149 * -- LAPACK is a software package provided by Univ. of Tennessee, --
150 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
151 * September 2012
152 *
153 * .. Scalar Arguments ..
154  INTEGER info, iw, n
155  DOUBLE PRECISION pivmin, reltol, gl, gu, w, werr
156 * ..
157 * .. Array Arguments ..
158  DOUBLE PRECISION d( * ), e2( * )
159 * ..
160 *
161 * =====================================================================
162 *
163 * .. Parameters ..
164  DOUBLE PRECISION fudge, half, two, zero
165  parameter( half = 0.5d0, two = 2.0d0,
166  $ fudge = two, zero = 0.0d0 )
167 * ..
168 * .. Local Scalars ..
169  INTEGER i, it, itmax, negcnt
170  DOUBLE PRECISION atoli, eps, left, mid, right, rtoli, tmp1,
171  $ tmp2, tnorm
172 * ..
173 * .. External Functions ..
174  DOUBLE PRECISION dlamch
175  EXTERNAL dlamch
176 * ..
177 * .. Intrinsic Functions ..
178  INTRINSIC abs, int, log, max
179 * ..
180 * .. Executable Statements ..
181 *
182 * Get machine constants
183  eps = dlamch( 'P' )
184 
185  tnorm = max( abs( gl ), abs( gu ) )
186  rtoli = reltol
187  atoli = fudge*two*pivmin
188 
189  itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /
190  $ log( two ) ) + 2
191 
192  info = -1
193 
194  left = gl - fudge*tnorm*eps*n - fudge*two*pivmin
195  right = gu + fudge*tnorm*eps*n + fudge*two*pivmin
196  it = 0
197 
198  10 CONTINUE
199 *
200 * Check if interval converged or maximum number of iterations reached
201 *
202  tmp1 = abs( right - left )
203  tmp2 = max( abs(right), abs(left) )
204  IF( tmp1.LT.max( atoli, pivmin, rtoli*tmp2 ) ) THEN
205  info = 0
206  GOTO 30
207  ENDIF
208  IF(it.GT.itmax)
209  $ GOTO 30
210 
211 *
212 * Count number of negative pivots for mid-point
213 *
214  it = it + 1
215  mid = half * (left + right)
216  negcnt = 0
217  tmp1 = d( 1 ) - mid
218  IF( abs( tmp1 ).LT.pivmin )
219  $ tmp1 = -pivmin
220  IF( tmp1.LE.zero )
221  $ negcnt = negcnt + 1
222 *
223  DO 20 i = 2, n
224  tmp1 = d( i ) - e2( i-1 ) / tmp1 - mid
225  IF( abs( tmp1 ).LT.pivmin )
226  $ tmp1 = -pivmin
227  IF( tmp1.LE.zero )
228  $ negcnt = negcnt + 1
229  20 CONTINUE
230 
231  IF(negcnt.GE.iw) THEN
232  right = mid
233  ELSE
234  left = mid
235  ENDIF
236  GOTO 10
237 
238  30 CONTINUE
239 *
240 * Converged or maximum number of iterations reached
241 *
242  w = half * (left + right)
243  werr = half * abs( right - left )
244 
245  RETURN
246 *
247 * End of DLARRK
248 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65

Here is the caller graph for this function:

subroutine dlarrr ( integer  N,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
integer  INFO 
)

DLARRR performs tests to decide whether the symmetric tridiagonal matrix T warrants expensive computations which guarantee high relative accuracy in the eigenvalues.

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

Purpose:
 Perform tests to decide whether the symmetric tridiagonal matrix T
 warrants expensive computations which guarantee high relative accuracy
 in the eigenvalues.
Parameters
[in]N
          N is INTEGER
          The order of the matrix. N > 0.
[in]D
          D is DOUBLE PRECISION array, dimension (N)
          The N diagonal elements of the tridiagonal matrix T.
[in,out]E
          E is DOUBLE PRECISION array, dimension (N)
          On entry, the first (N-1) entries contain the subdiagonal
          elements of the tridiagonal matrix T; E(N) is set to ZERO.
[out]INFO
          INFO is INTEGER
          INFO = 0(default) : the matrix warrants computations preserving
                              relative accuracy.
          INFO = 1          : the matrix warrants computations guaranteeing
                              only absolute accuracy.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Contributors:
Beresford Parlett, University of California, Berkeley, USA
Jim Demmel, University of California, Berkeley, USA
Inderjit Dhillon, University of Texas, Austin, USA
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA

Definition at line 96 of file dlarrr.f.

96 *
97 * -- LAPACK auxiliary routine (version 3.4.2) --
98 * -- LAPACK is a software package provided by Univ. of Tennessee, --
99 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
100 * September 2012
101 *
102 * .. Scalar Arguments ..
103  INTEGER n, info
104 * ..
105 * .. Array Arguments ..
106  DOUBLE PRECISION d( * ), e( * )
107 * ..
108 *
109 *
110 * =====================================================================
111 *
112 * .. Parameters ..
113  DOUBLE PRECISION zero, relcond
114  parameter( zero = 0.0d0,
115  $ relcond = 0.999d0 )
116 * ..
117 * .. Local Scalars ..
118  INTEGER i
119  LOGICAL yesrel
120  DOUBLE PRECISION eps, safmin, smlnum, rmin, tmp, tmp2,
121  $ offdig, offdig2
122 
123 * ..
124 * .. External Functions ..
125  DOUBLE PRECISION dlamch
126  EXTERNAL dlamch
127 * ..
128 * .. Intrinsic Functions ..
129  INTRINSIC abs
130 * ..
131 * .. Executable Statements ..
132 *
133 * As a default, do NOT go for relative-accuracy preserving computations.
134  info = 1
135 
136  safmin = dlamch( 'Safe minimum' )
137  eps = dlamch( 'Precision' )
138  smlnum = safmin / eps
139  rmin = sqrt( smlnum )
140 
141 * Tests for relative accuracy
142 *
143 * Test for scaled diagonal dominance
144 * Scale the diagonal entries to one and check whether the sum of the
145 * off-diagonals is less than one
146 *
147 * The sdd relative error bounds have a 1/(1- 2*x) factor in them,
148 * x = max(OFFDIG + OFFDIG2), so when x is close to 1/2, no relative
149 * accuracy is promised. In the notation of the code fragment below,
150 * 1/(1 - (OFFDIG + OFFDIG2)) is the condition number.
151 * We don't think it is worth going into "sdd mode" unless the relative
152 * condition number is reasonable, not 1/macheps.
153 * The threshold should be compatible with other thresholds used in the
154 * code. We set OFFDIG + OFFDIG2 <= .999 =: RELCOND, it corresponds
155 * to losing at most 3 decimal digits: 1 / (1 - (OFFDIG + OFFDIG2)) <= 1000
156 * instead of the current OFFDIG + OFFDIG2 < 1
157 *
158  yesrel = .true.
159  offdig = zero
160  tmp = sqrt(abs(d(1)))
161  IF (tmp.LT.rmin) yesrel = .false.
162  IF(.NOT.yesrel) GOTO 11
163  DO 10 i = 2, n
164  tmp2 = sqrt(abs(d(i)))
165  IF (tmp2.LT.rmin) yesrel = .false.
166  IF(.NOT.yesrel) GOTO 11
167  offdig2 = abs(e(i-1))/(tmp*tmp2)
168  IF(offdig+offdig2.GE.relcond) yesrel = .false.
169  IF(.NOT.yesrel) GOTO 11
170  tmp = tmp2
171  offdig = offdig2
172  10 CONTINUE
173  11 CONTINUE
174 
175  IF( yesrel ) THEN
176  info = 0
177  RETURN
178  ELSE
179  ENDIF
180 *
181 
182 *
183 * *** MORE TO BE IMPLEMENTED ***
184 *
185 
186 *
187 * Test if the lower bidiagonal matrix L from T = L D L^T
188 * (zero shift facto) is well conditioned
189 *
190 
191 *
192 * Test if the upper bidiagonal matrix U from T = U D U^T
193 * (zero shift facto) is well conditioned.
194 * In this case, the matrix needs to be flipped and, at the end
195 * of the eigenvector computation, the flip needs to be applied
196 * to the computed eigenvectors (and the support)
197 *
198 
199 *
200  RETURN
201 *
202 * END OF DLARRR
203 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65

Here is the caller graph for this function:

subroutine dlartg ( double precision  F,
double precision  G,
double precision  CS,
double precision  SN,
double precision  R 
)

DLARTG generates a plane rotation with real cosine and real sine.

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

Purpose:
 DLARTG generate a plane rotation so that

    [  CS  SN  ]  .  [ F ]  =  [ R ]   where CS**2 + SN**2 = 1.
    [ -SN  CS  ]     [ G ]     [ 0 ]

 This is a slower, more accurate version of the BLAS1 routine DROTG,
 with the following other differences:
    F and G are unchanged on return.
    If G=0, then CS=1 and SN=0.
    If F=0 and (G .ne. 0), then CS=0 and SN=1 without doing any
       floating point operations (saves work in DBDSQR when
       there are zeros on the diagonal).

 If F exceeds G in magnitude, CS will be positive.
Parameters
[in]F
          F is DOUBLE PRECISION
          The first component of vector to be rotated.
[in]G
          G is DOUBLE PRECISION
          The second component of vector to be rotated.
[out]CS
          CS is DOUBLE PRECISION
          The cosine of the rotation.
[out]SN
          SN is DOUBLE PRECISION
          The sine of the rotation.
[out]R
          R is DOUBLE PRECISION
          The nonzero component of the rotated vector.

  This version has a few statements commented out for thread safety
  (machine parameters are computed on each entry). 10 feb 03, SJH.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 99 of file dlartg.f.

99 *
100 * -- LAPACK auxiliary routine (version 3.4.2) --
101 * -- LAPACK is a software package provided by Univ. of Tennessee, --
102 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
103 * September 2012
104 *
105 * .. Scalar Arguments ..
106  DOUBLE PRECISION cs, f, g, r, sn
107 * ..
108 *
109 * =====================================================================
110 *
111 * .. Parameters ..
112  DOUBLE PRECISION zero
113  parameter( zero = 0.0d0 )
114  DOUBLE PRECISION one
115  parameter( one = 1.0d0 )
116  DOUBLE PRECISION two
117  parameter( two = 2.0d0 )
118 * ..
119 * .. Local Scalars ..
120 * LOGICAL FIRST
121  INTEGER count, i
122  DOUBLE PRECISION eps, f1, g1, safmin, safmn2, safmx2, scale
123 * ..
124 * .. External Functions ..
125  DOUBLE PRECISION dlamch
126  EXTERNAL dlamch
127 * ..
128 * .. Intrinsic Functions ..
129  INTRINSIC abs, int, log, max, sqrt
130 * ..
131 * .. Save statement ..
132 * SAVE FIRST, SAFMX2, SAFMIN, SAFMN2
133 * ..
134 * .. Data statements ..
135 * DATA FIRST / .TRUE. /
136 * ..
137 * .. Executable Statements ..
138 *
139 * IF( FIRST ) THEN
140  safmin = dlamch( 'S' )
141  eps = dlamch( 'E' )
142  safmn2 = dlamch( 'B' )**int( log( safmin / eps ) /
143  $ log( dlamch( 'B' ) ) / two )
144  safmx2 = one / safmn2
145 * FIRST = .FALSE.
146 * END IF
147  IF( g.EQ.zero ) THEN
148  cs = one
149  sn = zero
150  r = f
151  ELSE IF( f.EQ.zero ) THEN
152  cs = zero
153  sn = one
154  r = g
155  ELSE
156  f1 = f
157  g1 = g
158  scale = max( abs( f1 ), abs( g1 ) )
159  IF( scale.GE.safmx2 ) THEN
160  count = 0
161  10 CONTINUE
162  count = count + 1
163  f1 = f1*safmn2
164  g1 = g1*safmn2
165  scale = max( abs( f1 ), abs( g1 ) )
166  IF( scale.GE.safmx2 )
167  $ GO TO 10
168  r = sqrt( f1**2+g1**2 )
169  cs = f1 / r
170  sn = g1 / r
171  DO 20 i = 1, count
172  r = r*safmx2
173  20 CONTINUE
174  ELSE IF( scale.LE.safmn2 ) THEN
175  count = 0
176  30 CONTINUE
177  count = count + 1
178  f1 = f1*safmx2
179  g1 = g1*safmx2
180  scale = max( abs( f1 ), abs( g1 ) )
181  IF( scale.LE.safmn2 )
182  $ GO TO 30
183  r = sqrt( f1**2+g1**2 )
184  cs = f1 / r
185  sn = g1 / r
186  DO 40 i = 1, count
187  r = r*safmn2
188  40 CONTINUE
189  ELSE
190  r = sqrt( f1**2+g1**2 )
191  cs = f1 / r
192  sn = g1 / r
193  END IF
194  IF( abs( f ).GT.abs( g ) .AND. cs.LT.zero ) THEN
195  cs = -cs
196  sn = -sn
197  r = -r
198  END IF
199  END IF
200  RETURN
201 *
202 * End of DLARTG
203 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65

Here is the caller graph for this function:

subroutine dlartgp ( double precision  F,
double precision  G,
double precision  CS,
double precision  SN,
double precision  R 
)

DLARTGP generates a plane rotation so that the diagonal is nonnegative.

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

Purpose:
 DLARTGP generates a plane rotation so that

    [  CS  SN  ]  .  [ F ]  =  [ R ]   where CS**2 + SN**2 = 1.
    [ -SN  CS  ]     [ G ]     [ 0 ]

 This is a slower, more accurate version of the Level 1 BLAS routine DROTG,
 with the following other differences:
    F and G are unchanged on return.
    If G=0, then CS=(+/-)1 and SN=0.
    If F=0 and (G .ne. 0), then CS=0 and SN=(+/-)1.

 The sign is chosen so that R >= 0.
Parameters
[in]F
          F is DOUBLE PRECISION
          The first component of vector to be rotated.
[in]G
          G is DOUBLE PRECISION
          The second component of vector to be rotated.
[out]CS
          CS is DOUBLE PRECISION
          The cosine of the rotation.
[out]SN
          SN is DOUBLE PRECISION
          The sine of the rotation.
[out]R
          R is DOUBLE PRECISION
          The nonzero component of the rotated vector.

  This version has a few statements commented out for thread safety
  (machine parameters are computed on each entry). 10 feb 03, SJH.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 97 of file dlartgp.f.

97 *
98 * -- LAPACK auxiliary routine (version 3.4.2) --
99 * -- LAPACK is a software package provided by Univ. of Tennessee, --
100 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
101 * September 2012
102 *
103 * .. Scalar Arguments ..
104  DOUBLE PRECISION cs, f, g, r, sn
105 * ..
106 *
107 * =====================================================================
108 *
109 * .. Parameters ..
110  DOUBLE PRECISION zero
111  parameter( zero = 0.0d0 )
112  DOUBLE PRECISION one
113  parameter( one = 1.0d0 )
114  DOUBLE PRECISION two
115  parameter( two = 2.0d0 )
116 * ..
117 * .. Local Scalars ..
118 * LOGICAL FIRST
119  INTEGER count, i
120  DOUBLE PRECISION eps, f1, g1, safmin, safmn2, safmx2, scale
121 * ..
122 * .. External Functions ..
123  DOUBLE PRECISION dlamch
124  EXTERNAL dlamch
125 * ..
126 * .. Intrinsic Functions ..
127  INTRINSIC abs, int, log, max, sign, sqrt
128 * ..
129 * .. Save statement ..
130 * SAVE FIRST, SAFMX2, SAFMIN, SAFMN2
131 * ..
132 * .. Data statements ..
133 * DATA FIRST / .TRUE. /
134 * ..
135 * .. Executable Statements ..
136 *
137 * IF( FIRST ) THEN
138  safmin = dlamch( 'S' )
139  eps = dlamch( 'E' )
140  safmn2 = dlamch( 'B' )**int( log( safmin / eps ) /
141  $ log( dlamch( 'B' ) ) / two )
142  safmx2 = one / safmn2
143 * FIRST = .FALSE.
144 * END IF
145  IF( g.EQ.zero ) THEN
146  cs = sign( one, f )
147  sn = zero
148  r = abs( f )
149  ELSE IF( f.EQ.zero ) THEN
150  cs = zero
151  sn = sign( one, g )
152  r = abs( g )
153  ELSE
154  f1 = f
155  g1 = g
156  scale = max( abs( f1 ), abs( g1 ) )
157  IF( scale.GE.safmx2 ) THEN
158  count = 0
159  10 CONTINUE
160  count = count + 1
161  f1 = f1*safmn2
162  g1 = g1*safmn2
163  scale = max( abs( f1 ), abs( g1 ) )
164  IF( scale.GE.safmx2 )
165  $ GO TO 10
166  r = sqrt( f1**2+g1**2 )
167  cs = f1 / r
168  sn = g1 / r
169  DO 20 i = 1, count
170  r = r*safmx2
171  20 CONTINUE
172  ELSE IF( scale.LE.safmn2 ) THEN
173  count = 0
174  30 CONTINUE
175  count = count + 1
176  f1 = f1*safmx2
177  g1 = g1*safmx2
178  scale = max( abs( f1 ), abs( g1 ) )
179  IF( scale.LE.safmn2 )
180  $ GO TO 30
181  r = sqrt( f1**2+g1**2 )
182  cs = f1 / r
183  sn = g1 / r
184  DO 40 i = 1, count
185  r = r*safmn2
186  40 CONTINUE
187  ELSE
188  r = sqrt( f1**2+g1**2 )
189  cs = f1 / r
190  sn = g1 / r
191  END IF
192  IF( r.LT.zero ) THEN
193  cs = -cs
194  sn = -sn
195  r = -r
196  END IF
197  END IF
198  RETURN
199 *
200 * End of DLARTGP
201 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65

Here is the caller graph for this function:

subroutine dlaruv ( integer, dimension( 4 )  ISEED,
integer  N,
double precision, dimension( n )  X 
)

DLARUV returns a vector of n random real numbers from a uniform distribution.

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

Purpose:
 DLARUV returns a vector of n random real numbers from a uniform (0,1)
 distribution (n <= 128).

 This is an auxiliary routine called by DLARNV and ZLARNV.
Parameters
[in,out]ISEED
          ISEED is INTEGER array, dimension (4)
          On entry, the seed of the random number generator; the array
          elements must be between 0 and 4095, and ISEED(4) must be
          odd.
          On exit, the seed is updated.
[in]N
          N is INTEGER
          The number of random numbers to be generated. N <= 128.
[out]X
          X is DOUBLE PRECISION array, dimension (N)
          The generated random numbers.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
  This routine uses a multiplicative congruential method with modulus
  2**48 and multiplier 33952834046453 (see G.S.Fishman,
  'Multiplicative congruential random number generators with modulus
  2**b: an exhaustive analysis for b = 32 and a partial analysis for
  b = 48', Math. Comp. 189, pp 331-344, 1990).

  48-bit integers are stored in 4 integer array elements with 12 bits
  per element. Hence the routine is portable across machines with
  integers of 32 bits or more.

Definition at line 97 of file dlaruv.f.

97 *
98 * -- LAPACK auxiliary routine (version 3.4.2) --
99 * -- LAPACK is a software package provided by Univ. of Tennessee, --
100 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
101 * September 2012
102 *
103 * .. Scalar Arguments ..
104  INTEGER n
105 * ..
106 * .. Array Arguments ..
107  INTEGER iseed( 4 )
108  DOUBLE PRECISION x( n )
109 * ..
110 *
111 * =====================================================================
112 *
113 * .. Parameters ..
114  DOUBLE PRECISION one
115  parameter( one = 1.0d0 )
116  INTEGER lv, ipw2
117  DOUBLE PRECISION r
118  parameter( lv = 128, ipw2 = 4096, r = one / ipw2 )
119 * ..
120 * .. Local Scalars ..
121  INTEGER i, i1, i2, i3, i4, it1, it2, it3, it4, j
122 * ..
123 * .. Local Arrays ..
124  INTEGER mm( lv, 4 )
125 * ..
126 * .. Intrinsic Functions ..
127  INTRINSIC dble, min, mod
128 * ..
129 * .. Data statements ..
130  DATA ( mm( 1, j ), j = 1, 4 ) / 494, 322, 2508,
131  $ 2549 /
132  DATA ( mm( 2, j ), j = 1, 4 ) / 2637, 789, 3754,
133  $ 1145 /
134  DATA ( mm( 3, j ), j = 1, 4 ) / 255, 1440, 1766,
135  $ 2253 /
136  DATA ( mm( 4, j ), j = 1, 4 ) / 2008, 752, 3572,
137  $ 305 /
138  DATA ( mm( 5, j ), j = 1, 4 ) / 1253, 2859, 2893,
139  $ 3301 /
140  DATA ( mm( 6, j ), j = 1, 4 ) / 3344, 123, 307,
141  $ 1065 /
142  DATA ( mm( 7, j ), j = 1, 4 ) / 4084, 1848, 1297,
143  $ 3133 /
144  DATA ( mm( 8, j ), j = 1, 4 ) / 1739, 643, 3966,
145  $ 2913 /
146  DATA ( mm( 9, j ), j = 1, 4 ) / 3143, 2405, 758,
147  $ 3285 /
148  DATA ( mm( 10, j ), j = 1, 4 ) / 3468, 2638, 2598,
149  $ 1241 /
150  DATA ( mm( 11, j ), j = 1, 4 ) / 688, 2344, 3406,
151  $ 1197 /
152  DATA ( mm( 12, j ), j = 1, 4 ) / 1657, 46, 2922,
153  $ 3729 /
154  DATA ( mm( 13, j ), j = 1, 4 ) / 1238, 3814, 1038,
155  $ 2501 /
156  DATA ( mm( 14, j ), j = 1, 4 ) / 3166, 913, 2934,
157  $ 1673 /
158  DATA ( mm( 15, j ), j = 1, 4 ) / 1292, 3649, 2091,
159  $ 541 /
160  DATA ( mm( 16, j ), j = 1, 4 ) / 3422, 339, 2451,
161  $ 2753 /
162  DATA ( mm( 17, j ), j = 1, 4 ) / 1270, 3808, 1580,
163  $ 949 /
164  DATA ( mm( 18, j ), j = 1, 4 ) / 2016, 822, 1958,
165  $ 2361 /
166  DATA ( mm( 19, j ), j = 1, 4 ) / 154, 2832, 2055,
167  $ 1165 /
168  DATA ( mm( 20, j ), j = 1, 4 ) / 2862, 3078, 1507,
169  $ 4081 /
170  DATA ( mm( 21, j ), j = 1, 4 ) / 697, 3633, 1078,
171  $ 2725 /
172  DATA ( mm( 22, j ), j = 1, 4 ) / 1706, 2970, 3273,
173  $ 3305 /
174  DATA ( mm( 23, j ), j = 1, 4 ) / 491, 637, 17,
175  $ 3069 /
176  DATA ( mm( 24, j ), j = 1, 4 ) / 931, 2249, 854,
177  $ 3617 /
178  DATA ( mm( 25, j ), j = 1, 4 ) / 1444, 2081, 2916,
179  $ 3733 /
180  DATA ( mm( 26, j ), j = 1, 4 ) / 444, 4019, 3971,
181  $ 409 /
182  DATA ( mm( 27, j ), j = 1, 4 ) / 3577, 1478, 2889,
183  $ 2157 /
184  DATA ( mm( 28, j ), j = 1, 4 ) / 3944, 242, 3831,
185  $ 1361 /
186  DATA ( mm( 29, j ), j = 1, 4 ) / 2184, 481, 2621,
187  $ 3973 /
188  DATA ( mm( 30, j ), j = 1, 4 ) / 1661, 2075, 1541,
189  $ 1865 /
190  DATA ( mm( 31, j ), j = 1, 4 ) / 3482, 4058, 893,
191  $ 2525 /
192  DATA ( mm( 32, j ), j = 1, 4 ) / 657, 622, 736,
193  $ 1409 /
194  DATA ( mm( 33, j ), j = 1, 4 ) / 3023, 3376, 3992,
195  $ 3445 /
196  DATA ( mm( 34, j ), j = 1, 4 ) / 3618, 812, 787,
197  $ 3577 /
198  DATA ( mm( 35, j ), j = 1, 4 ) / 1267, 234, 2125,
199  $ 77 /
200  DATA ( mm( 36, j ), j = 1, 4 ) / 1828, 641, 2364,
201  $ 3761 /
202  DATA ( mm( 37, j ), j = 1, 4 ) / 164, 4005, 2460,
203  $ 2149 /
204  DATA ( mm( 38, j ), j = 1, 4 ) / 3798, 1122, 257,
205  $ 1449 /
206  DATA ( mm( 39, j ), j = 1, 4 ) / 3087, 3135, 1574,
207  $ 3005 /
208  DATA ( mm( 40, j ), j = 1, 4 ) / 2400, 2640, 3912,
209  $ 225 /
210  DATA ( mm( 41, j ), j = 1, 4 ) / 2870, 2302, 1216,
211  $ 85 /
212  DATA ( mm( 42, j ), j = 1, 4 ) / 3876, 40, 3248,
213  $ 3673 /
214  DATA ( mm( 43, j ), j = 1, 4 ) / 1905, 1832, 3401,
215  $ 3117 /
216  DATA ( mm( 44, j ), j = 1, 4 ) / 1593, 2247, 2124,
217  $ 3089 /
218  DATA ( mm( 45, j ), j = 1, 4 ) / 1797, 2034, 2762,
219  $ 1349 /
220  DATA ( mm( 46, j ), j = 1, 4 ) / 1234, 2637, 149,
221  $ 2057 /
222  DATA ( mm( 47, j ), j = 1, 4 ) / 3460, 1287, 2245,
223  $ 413 /
224  DATA ( mm( 48, j ), j = 1, 4 ) / 328, 1691, 166,
225  $ 65 /
226  DATA ( mm( 49, j ), j = 1, 4 ) / 2861, 496, 466,
227  $ 1845 /
228  DATA ( mm( 50, j ), j = 1, 4 ) / 1950, 1597, 4018,
229  $ 697 /
230  DATA ( mm( 51, j ), j = 1, 4 ) / 617, 2394, 1399,
231  $ 3085 /
232  DATA ( mm( 52, j ), j = 1, 4 ) / 2070, 2584, 190,
233  $ 3441 /
234  DATA ( mm( 53, j ), j = 1, 4 ) / 3331, 1843, 2879,
235  $ 1573 /
236  DATA ( mm( 54, j ), j = 1, 4 ) / 769, 336, 153,
237  $ 3689 /
238  DATA ( mm( 55, j ), j = 1, 4 ) / 1558, 1472, 2320,
239  $ 2941 /
240  DATA ( mm( 56, j ), j = 1, 4 ) / 2412, 2407, 18,
241  $ 929 /
242  DATA ( mm( 57, j ), j = 1, 4 ) / 2800, 433, 712,
243  $ 533 /
244  DATA ( mm( 58, j ), j = 1, 4 ) / 189, 2096, 2159,
245  $ 2841 /
246  DATA ( mm( 59, j ), j = 1, 4 ) / 287, 1761, 2318,
247  $ 4077 /
248  DATA ( mm( 60, j ), j = 1, 4 ) / 2045, 2810, 2091,
249  $ 721 /
250  DATA ( mm( 61, j ), j = 1, 4 ) / 1227, 566, 3443,
251  $ 2821 /
252  DATA ( mm( 62, j ), j = 1, 4 ) / 2838, 442, 1510,
253  $ 2249 /
254  DATA ( mm( 63, j ), j = 1, 4 ) / 209, 41, 449,
255  $ 2397 /
256  DATA ( mm( 64, j ), j = 1, 4 ) / 2770, 1238, 1956,
257  $ 2817 /
258  DATA ( mm( 65, j ), j = 1, 4 ) / 3654, 1086, 2201,
259  $ 245 /
260  DATA ( mm( 66, j ), j = 1, 4 ) / 3993, 603, 3137,
261  $ 1913 /
262  DATA ( mm( 67, j ), j = 1, 4 ) / 192, 840, 3399,
263  $ 1997 /
264  DATA ( mm( 68, j ), j = 1, 4 ) / 2253, 3168, 1321,
265  $ 3121 /
266  DATA ( mm( 69, j ), j = 1, 4 ) / 3491, 1499, 2271,
267  $ 997 /
268  DATA ( mm( 70, j ), j = 1, 4 ) / 2889, 1084, 3667,
269  $ 1833 /
270  DATA ( mm( 71, j ), j = 1, 4 ) / 2857, 3438, 2703,
271  $ 2877 /
272  DATA ( mm( 72, j ), j = 1, 4 ) / 2094, 2408, 629,
273  $ 1633 /
274  DATA ( mm( 73, j ), j = 1, 4 ) / 1818, 1589, 2365,
275  $ 981 /
276  DATA ( mm( 74, j ), j = 1, 4 ) / 688, 2391, 2431,
277  $ 2009 /
278  DATA ( mm( 75, j ), j = 1, 4 ) / 1407, 288, 1113,
279  $ 941 /
280  DATA ( mm( 76, j ), j = 1, 4 ) / 634, 26, 3922,
281  $ 2449 /
282  DATA ( mm( 77, j ), j = 1, 4 ) / 3231, 512, 2554,
283  $ 197 /
284  DATA ( mm( 78, j ), j = 1, 4 ) / 815, 1456, 184,
285  $ 2441 /
286  DATA ( mm( 79, j ), j = 1, 4 ) / 3524, 171, 2099,
287  $ 285 /
288  DATA ( mm( 80, j ), j = 1, 4 ) / 1914, 1677, 3228,
289  $ 1473 /
290  DATA ( mm( 81, j ), j = 1, 4 ) / 516, 2657, 4012,
291  $ 2741 /
292  DATA ( mm( 82, j ), j = 1, 4 ) / 164, 2270, 1921,
293  $ 3129 /
294  DATA ( mm( 83, j ), j = 1, 4 ) / 303, 2587, 3452,
295  $ 909 /
296  DATA ( mm( 84, j ), j = 1, 4 ) / 2144, 2961, 3901,
297  $ 2801 /
298  DATA ( mm( 85, j ), j = 1, 4 ) / 3480, 1970, 572,
299  $ 421 /
300  DATA ( mm( 86, j ), j = 1, 4 ) / 119, 1817, 3309,
301  $ 4073 /
302  DATA ( mm( 87, j ), j = 1, 4 ) / 3357, 676, 3171,
303  $ 2813 /
304  DATA ( mm( 88, j ), j = 1, 4 ) / 837, 1410, 817,
305  $ 2337 /
306  DATA ( mm( 89, j ), j = 1, 4 ) / 2826, 3723, 3039,
307  $ 1429 /
308  DATA ( mm( 90, j ), j = 1, 4 ) / 2332, 2803, 1696,
309  $ 1177 /
310  DATA ( mm( 91, j ), j = 1, 4 ) / 2089, 3185, 1256,
311  $ 1901 /
312  DATA ( mm( 92, j ), j = 1, 4 ) / 3780, 184, 3715,
313  $ 81 /
314  DATA ( mm( 93, j ), j = 1, 4 ) / 1700, 663, 2077,
315  $ 1669 /
316  DATA ( mm( 94, j ), j = 1, 4 ) / 3712, 499, 3019,
317  $ 2633 /
318  DATA ( mm( 95, j ), j = 1, 4 ) / 150, 3784, 1497,
319  $ 2269 /
320  DATA ( mm( 96, j ), j = 1, 4 ) / 2000, 1631, 1101,
321  $ 129 /
322  DATA ( mm( 97, j ), j = 1, 4 ) / 3375, 1925, 717,
323  $ 1141 /
324  DATA ( mm( 98, j ), j = 1, 4 ) / 1621, 3912, 51,
325  $ 249 /
326  DATA ( mm( 99, j ), j = 1, 4 ) / 3090, 1398, 981,
327  $ 3917 /
328  DATA ( mm( 100, j ), j = 1, 4 ) / 3765, 1349, 1978,
329  $ 2481 /
330  DATA ( mm( 101, j ), j = 1, 4 ) / 1149, 1441, 1813,
331  $ 3941 /
332  DATA ( mm( 102, j ), j = 1, 4 ) / 3146, 2224, 3881,
333  $ 2217 /
334  DATA ( mm( 103, j ), j = 1, 4 ) / 33, 2411, 76,
335  $ 2749 /
336  DATA ( mm( 104, j ), j = 1, 4 ) / 3082, 1907, 3846,
337  $ 3041 /
338  DATA ( mm( 105, j ), j = 1, 4 ) / 2741, 3192, 3694,
339  $ 1877 /
340  DATA ( mm( 106, j ), j = 1, 4 ) / 359, 2786, 1682,
341  $ 345 /
342  DATA ( mm( 107, j ), j = 1, 4 ) / 3316, 382, 124,
343  $ 2861 /
344  DATA ( mm( 108, j ), j = 1, 4 ) / 1749, 37, 1660,
345  $ 1809 /
346  DATA ( mm( 109, j ), j = 1, 4 ) / 185, 759, 3997,
347  $ 3141 /
348  DATA ( mm( 110, j ), j = 1, 4 ) / 2784, 2948, 479,
349  $ 2825 /
350  DATA ( mm( 111, j ), j = 1, 4 ) / 2202, 1862, 1141,
351  $ 157 /
352  DATA ( mm( 112, j ), j = 1, 4 ) / 2199, 3802, 886,
353  $ 2881 /
354  DATA ( mm( 113, j ), j = 1, 4 ) / 1364, 2423, 3514,
355  $ 3637 /
356  DATA ( mm( 114, j ), j = 1, 4 ) / 1244, 2051, 1301,
357  $ 1465 /
358  DATA ( mm( 115, j ), j = 1, 4 ) / 2020, 2295, 3604,
359  $ 2829 /
360  DATA ( mm( 116, j ), j = 1, 4 ) / 3160, 1332, 1888,
361  $ 2161 /
362  DATA ( mm( 117, j ), j = 1, 4 ) / 2785, 1832, 1836,
363  $ 3365 /
364  DATA ( mm( 118, j ), j = 1, 4 ) / 2772, 2405, 1990,
365  $ 361 /
366  DATA ( mm( 119, j ), j = 1, 4 ) / 1217, 3638, 2058,
367  $ 2685 /
368  DATA ( mm( 120, j ), j = 1, 4 ) / 1822, 3661, 692,
369  $ 3745 /
370  DATA ( mm( 121, j ), j = 1, 4 ) / 1245, 327, 1194,
371  $ 2325 /
372  DATA ( mm( 122, j ), j = 1, 4 ) / 2252, 3660, 20,
373  $ 3609 /
374  DATA ( mm( 123, j ), j = 1, 4 ) / 3904, 716, 3285,
375  $ 3821 /
376  DATA ( mm( 124, j ), j = 1, 4 ) / 2774, 1842, 2046,
377  $ 3537 /
378  DATA ( mm( 125, j ), j = 1, 4 ) / 997, 3987, 2107,
379  $ 517 /
380  DATA ( mm( 126, j ), j = 1, 4 ) / 2573, 1368, 3508,
381  $ 3017 /
382  DATA ( mm( 127, j ), j = 1, 4 ) / 1148, 1848, 3525,
383  $ 2141 /
384  DATA ( mm( 128, j ), j = 1, 4 ) / 545, 2366, 3801,
385  $ 1537 /
386 * ..
387 * .. Executable Statements ..
388 *
389  i1 = iseed( 1 )
390  i2 = iseed( 2 )
391  i3 = iseed( 3 )
392  i4 = iseed( 4 )
393 *
394  DO 10 i = 1, min( n, lv )
395 *
396  20 CONTINUE
397 *
398 * Multiply the seed by i-th power of the multiplier modulo 2**48
399 *
400  it4 = i4*mm( i, 4 )
401  it3 = it4 / ipw2
402  it4 = it4 - ipw2*it3
403  it3 = it3 + i3*mm( i, 4 ) + i4*mm( i, 3 )
404  it2 = it3 / ipw2
405  it3 = it3 - ipw2*it2
406  it2 = it2 + i2*mm( i, 4 ) + i3*mm( i, 3 ) + i4*mm( i, 2 )
407  it1 = it2 / ipw2
408  it2 = it2 - ipw2*it1
409  it1 = it1 + i1*mm( i, 4 ) + i2*mm( i, 3 ) + i3*mm( i, 2 ) +
410  $ i4*mm( i, 1 )
411  it1 = mod( it1, ipw2 )
412 *
413 * Convert 48-bit integer to a real number in the interval (0,1)
414 *
415  x( i ) = r*( dble( it1 )+r*( dble( it2 )+r*( dble( it3 )+r*
416  $ dble( it4 ) ) ) )
417 *
418  IF (x( i ).EQ.1.0d0) THEN
419 * If a real number has n bits of precision, and the first
420 * n bits of the 48-bit integer above happen to be all 1 (which
421 * will occur about once every 2**n calls), then X( I ) will
422 * be rounded to exactly 1.0.
423 * Since X( I ) is not supposed to return exactly 0.0 or 1.0,
424 * the statistically correct thing to do in this situation is
425 * simply to iterate again.
426 * N.B. the case X( I ) = 0.0 should not be possible.
427  i1 = i1 + 2
428  i2 = i2 + 2
429  i3 = i3 + 2
430  i4 = i4 + 2
431  GOTO 20
432  END IF
433 *
434  10 CONTINUE
435 *
436 * Return final value of seed
437 *
438  iseed( 1 ) = it1
439  iseed( 2 ) = it2
440  iseed( 3 ) = it3
441  iseed( 4 ) = it4
442  RETURN
443 *
444 * End of DLARUV
445 *

Here is the caller graph for this function:

subroutine dlas2 ( double precision  F,
double precision  G,
double precision  H,
double precision  SSMIN,
double precision  SSMAX 
)

DLAS2 computes singular values of a 2-by-2 triangular matrix.

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

Purpose:
 DLAS2  computes the singular values of the 2-by-2 matrix
    [  F   G  ]
    [  0   H  ].
 On return, SSMIN is the smaller singular value and SSMAX is the
 larger singular value.
Parameters
[in]F
          F is DOUBLE PRECISION
          The (1,1) element of the 2-by-2 matrix.
[in]G
          G is DOUBLE PRECISION
          The (1,2) element of the 2-by-2 matrix.
[in]H
          H is DOUBLE PRECISION
          The (2,2) element of the 2-by-2 matrix.
[out]SSMIN
          SSMIN is DOUBLE PRECISION
          The smaller singular value.
[out]SSMAX
          SSMAX is DOUBLE PRECISION
          The larger singular value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
  Barring over/underflow, all output quantities are correct to within
  a few units in the last place (ulps), even in the absence of a guard
  digit in addition/subtraction.

  In IEEE arithmetic, the code works correctly if one matrix element is
  infinite.

  Overflow will not occur unless the largest singular value itself
  overflows, or is within a few ulps of overflow. (On machines with
  partial overflow, like the Cray, overflow may occur if the largest
  singular value is within a factor of 2 of overflow.)

  Underflow is harmless if underflow is gradual. Otherwise, results
  may correspond to a matrix modified by perturbations of size near
  the underflow threshold.

Definition at line 109 of file dlas2.f.

109 *
110 * -- LAPACK auxiliary routine (version 3.4.2) --
111 * -- LAPACK is a software package provided by Univ. of Tennessee, --
112 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
113 * September 2012
114 *
115 * .. Scalar Arguments ..
116  DOUBLE PRECISION f, g, h, ssmax, ssmin
117 * ..
118 *
119 * ====================================================================
120 *
121 * .. Parameters ..
122  DOUBLE PRECISION zero
123  parameter( zero = 0.0d0 )
124  DOUBLE PRECISION one
125  parameter( one = 1.0d0 )
126  DOUBLE PRECISION two
127  parameter( two = 2.0d0 )
128 * ..
129 * .. Local Scalars ..
130  DOUBLE PRECISION as, at, au, c, fa, fhmn, fhmx, ga, ha
131 * ..
132 * .. Intrinsic Functions ..
133  INTRINSIC abs, max, min, sqrt
134 * ..
135 * .. Executable Statements ..
136 *
137  fa = abs( f )
138  ga = abs( g )
139  ha = abs( h )
140  fhmn = min( fa, ha )
141  fhmx = max( fa, ha )
142  IF( fhmn.EQ.zero ) THEN
143  ssmin = zero
144  IF( fhmx.EQ.zero ) THEN
145  ssmax = ga
146  ELSE
147  ssmax = max( fhmx, ga )*sqrt( one+
148  $ ( min( fhmx, ga ) / max( fhmx, ga ) )**2 )
149  END IF
150  ELSE
151  IF( ga.LT.fhmx ) THEN
152  as = one + fhmn / fhmx
153  at = ( fhmx-fhmn ) / fhmx
154  au = ( ga / fhmx )**2
155  c = two / ( sqrt( as*as+au )+sqrt( at*at+au ) )
156  ssmin = fhmn*c
157  ssmax = fhmx / c
158  ELSE
159  au = fhmx / ga
160  IF( au.EQ.zero ) THEN
161 *
162 * Avoid possible harmful underflow if exponent range
163 * asymmetric (true SSMIN may not underflow even if
164 * AU underflows)
165 *
166  ssmin = ( fhmn*fhmx ) / ga
167  ssmax = ga
168  ELSE
169  as = one + fhmn / fhmx
170  at = ( fhmx-fhmn ) / fhmx
171  c = one / ( sqrt( one+( as*au )**2 )+
172  $ sqrt( one+( at*au )**2 ) )
173  ssmin = ( fhmn*c )*au
174  ssmin = ssmin + ssmin
175  ssmax = ga / ( c+c )
176  END IF
177  END IF
178  END IF
179  RETURN
180 *
181 * End of DLAS2
182 *

Here is the caller graph for this function:

subroutine dlascl ( character  TYPE,
integer  KL,
integer  KU,
double precision  CFROM,
double precision  CTO,
integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer  INFO 
)

DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.

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

Purpose:
 DLASCL multiplies the M by N real matrix A by the real scalar
 CTO/CFROM.  This is done without over/underflow as long as the final
 result CTO*A(I,J)/CFROM does not over/underflow. TYPE specifies that
 A may be full, upper triangular, lower triangular, upper Hessenberg,
 or banded.
Parameters
[in]TYPE
          TYPE is CHARACTER*1
          TYPE indices the storage type of the input matrix.
          = 'G':  A is a full matrix.
          = 'L':  A is a lower triangular matrix.
          = 'U':  A is an upper triangular matrix.
          = 'H':  A is an upper Hessenberg matrix.
          = 'B':  A is a symmetric band matrix with lower bandwidth KL
                  and upper bandwidth KU and with the only the lower
                  half stored.
          = 'Q':  A is a symmetric band matrix with lower bandwidth KL
                  and upper bandwidth KU and with the only the upper
                  half stored.
          = 'Z':  A is a band matrix with lower bandwidth KL and upper
                  bandwidth KU. See DGBTRF for storage details.
[in]KL
          KL is INTEGER
          The lower bandwidth of A.  Referenced only if TYPE = 'B',
          'Q' or 'Z'.
[in]KU
          KU is INTEGER
          The upper bandwidth of A.  Referenced only if TYPE = 'B',
          'Q' or 'Z'.
[in]CFROM
          CFROM is DOUBLE PRECISION
[in]CTO
          CTO is DOUBLE PRECISION

          The matrix A is multiplied by CTO/CFROM. A(I,J) is computed
          without over/underflow if the final result CTO*A(I,J)/CFROM
          can be represented without over/underflow.  CFROM must be
          nonzero.
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is DOUBLE PRECISION array, dimension (LDA,N)
          The matrix to be multiplied by CTO/CFROM.  See TYPE for the
          storage type.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]INFO
          INFO is INTEGER
          0  - successful exit
          <0 - if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 141 of file dlascl.f.

141 *
142 * -- LAPACK auxiliary routine (version 3.4.2) --
143 * -- LAPACK is a software package provided by Univ. of Tennessee, --
144 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
145 * September 2012
146 *
147 * .. Scalar Arguments ..
148  CHARACTER type
149  INTEGER info, kl, ku, lda, m, n
150  DOUBLE PRECISION cfrom, cto
151 * ..
152 * .. Array Arguments ..
153  DOUBLE PRECISION a( lda, * )
154 * ..
155 *
156 * =====================================================================
157 *
158 * .. Parameters ..
159  DOUBLE PRECISION zero, one
160  parameter( zero = 0.0d0, one = 1.0d0 )
161 * ..
162 * .. Local Scalars ..
163  LOGICAL done
164  INTEGER i, itype, j, k1, k2, k3, k4
165  DOUBLE PRECISION bignum, cfrom1, cfromc, cto1, ctoc, mul, smlnum
166 * ..
167 * .. External Functions ..
168  LOGICAL lsame, disnan
169  DOUBLE PRECISION dlamch
170  EXTERNAL lsame, dlamch, disnan
171 * ..
172 * .. Intrinsic Functions ..
173  INTRINSIC abs, max, min
174 * ..
175 * .. External Subroutines ..
176  EXTERNAL xerbla
177 * ..
178 * .. Executable Statements ..
179 *
180 * Test the input arguments
181 *
182  info = 0
183 *
184  IF( lsame( TYPE, 'G' ) ) then
185  itype = 0
186  ELSE IF( lsame( TYPE, 'L' ) ) then
187  itype = 1
188  ELSE IF( lsame( TYPE, 'U' ) ) then
189  itype = 2
190  ELSE IF( lsame( TYPE, 'H' ) ) then
191  itype = 3
192  ELSE IF( lsame( TYPE, 'B' ) ) then
193  itype = 4
194  ELSE IF( lsame( TYPE, 'Q' ) ) then
195  itype = 5
196  ELSE IF( lsame( TYPE, 'Z' ) ) then
197  itype = 6
198  ELSE
199  itype = -1
200  END IF
201 *
202  IF( itype.EQ.-1 ) THEN
203  info = -1
204  ELSE IF( cfrom.EQ.zero .OR. disnan(cfrom) ) THEN
205  info = -4
206  ELSE IF( disnan(cto) ) THEN
207  info = -5
208  ELSE IF( m.LT.0 ) THEN
209  info = -6
210  ELSE IF( n.LT.0 .OR. ( itype.EQ.4 .AND. n.NE.m ) .OR.
211  $ ( itype.EQ.5 .AND. n.NE.m ) ) THEN
212  info = -7
213  ELSE IF( itype.LE.3 .AND. lda.LT.max( 1, m ) ) THEN
214  info = -9
215  ELSE IF( itype.GE.4 ) THEN
216  IF( kl.LT.0 .OR. kl.GT.max( m-1, 0 ) ) THEN
217  info = -2
218  ELSE IF( ku.LT.0 .OR. ku.GT.max( n-1, 0 ) .OR.
219  $ ( ( itype.EQ.4 .OR. itype.EQ.5 ) .AND. kl.NE.ku ) )
220  $ THEN
221  info = -3
222  ELSE IF( ( itype.EQ.4 .AND. lda.LT.kl+1 ) .OR.
223  $ ( itype.EQ.5 .AND. lda.LT.ku+1 ) .OR.
224  $ ( itype.EQ.6 .AND. lda.LT.2*kl+ku+1 ) ) THEN
225  info = -9
226  END IF
227  END IF
228 *
229  IF( info.NE.0 ) THEN
230  CALL xerbla( 'DLASCL', -info )
231  RETURN
232  END IF
233 *
234 * Quick return if possible
235 *
236  IF( n.EQ.0 .OR. m.EQ.0 )
237  $ RETURN
238 *
239 * Get machine parameters
240 *
241  smlnum = dlamch( 'S' )
242  bignum = one / smlnum
243 *
244  cfromc = cfrom
245  ctoc = cto
246 *
247  10 CONTINUE
248  cfrom1 = cfromc*smlnum
249  IF( cfrom1.EQ.cfromc ) THEN
250 ! CFROMC is an inf. Multiply by a correctly signed zero for
251 ! finite CTOC, or a NaN if CTOC is infinite.
252  mul = ctoc / cfromc
253  done = .true.
254  cto1 = ctoc
255  ELSE
256  cto1 = ctoc / bignum
257  IF( cto1.EQ.ctoc ) THEN
258 ! CTOC is either 0 or an inf. In both cases, CTOC itself
259 ! serves as the correct multiplication factor.
260  mul = ctoc
261  done = .true.
262  cfromc = one
263  ELSE IF( abs( cfrom1 ).GT.abs( ctoc ) .AND. ctoc.NE.zero ) THEN
264  mul = smlnum
265  done = .false.
266  cfromc = cfrom1
267  ELSE IF( abs( cto1 ).GT.abs( cfromc ) ) THEN
268  mul = bignum
269  done = .false.
270  ctoc = cto1
271  ELSE
272  mul = ctoc / cfromc
273  done = .true.
274  END IF
275  END IF
276 *
277  IF( itype.EQ.0 ) THEN
278 *
279 * Full matrix
280 *
281  DO 30 j = 1, n
282  DO 20 i = 1, m
283  a( i, j ) = a( i, j )*mul
284  20 CONTINUE
285  30 CONTINUE
286 *
287  ELSE IF( itype.EQ.1 ) THEN
288 *
289 * Lower triangular matrix
290 *
291  DO 50 j = 1, n
292  DO 40 i = j, m
293  a( i, j ) = a( i, j )*mul
294  40 CONTINUE
295  50 CONTINUE
296 *
297  ELSE IF( itype.EQ.2 ) THEN
298 *
299 * Upper triangular matrix
300 *
301  DO 70 j = 1, n
302  DO 60 i = 1, min( j, m )
303  a( i, j ) = a( i, j )*mul
304  60 CONTINUE
305  70 CONTINUE
306 *
307  ELSE IF( itype.EQ.3 ) THEN
308 *
309 * Upper Hessenberg matrix
310 *
311  DO 90 j = 1, n
312  DO 80 i = 1, min( j+1, m )
313  a( i, j ) = a( i, j )*mul
314  80 CONTINUE
315  90 CONTINUE
316 *
317  ELSE IF( itype.EQ.4 ) THEN
318 *
319 * Lower half of a symmetric band matrix
320 *
321  k3 = kl + 1
322  k4 = n + 1
323  DO 110 j = 1, n
324  DO 100 i = 1, min( k3, k4-j )
325  a( i, j ) = a( i, j )*mul
326  100 CONTINUE
327  110 CONTINUE
328 *
329  ELSE IF( itype.EQ.5 ) THEN
330 *
331 * Upper half of a symmetric band matrix
332 *
333  k1 = ku + 2
334  k3 = ku + 1
335  DO 130 j = 1, n
336  DO 120 i = max( k1-j, 1 ), k3
337  a( i, j ) = a( i, j )*mul
338  120 CONTINUE
339  130 CONTINUE
340 *
341  ELSE IF( itype.EQ.6 ) THEN
342 *
343 * Band matrix
344 *
345  k1 = kl + ku + 2
346  k2 = kl + 1
347  k3 = 2*kl + ku + 1
348  k4 = kl + ku + 1 + m
349  DO 150 j = 1, n
350  DO 140 i = max( k1-j, k2 ), min( k3, k4-j )
351  a( i, j ) = a( i, j )*mul
352  140 CONTINUE
353  150 CONTINUE
354 *
355  END IF
356 *
357  IF( .NOT.done )
358  $ GO TO 10
359 *
360  RETURN
361 *
362 * End of DLASCL
363 *
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:61
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

subroutine dlasd0 ( integer  N,
integer  SQRE,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( ldu, * )  U,
integer  LDU,
double precision, dimension( ldvt, * )  VT,
integer  LDVT,
integer  SMLSIZ,
integer, dimension( * )  IWORK,
double precision, dimension( * )  WORK,
integer  INFO 
)

DLASD0 computes the singular values of a real upper bidiagonal n-by-m matrix B with diagonal d and off-diagonal e. Used by sbdsdc.

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

Purpose:
 Using a divide and conquer approach, DLASD0 computes the singular
 value decomposition (SVD) of a real upper bidiagonal N-by-M
 matrix B with diagonal D and offdiagonal E, where M = N + SQRE.
 The algorithm computes orthogonal matrices U and VT such that
 B = U * S * VT. The singular values S are overwritten on D.

 A related subroutine, DLASDA, computes only the singular values,
 and optionally, the singular vectors in compact form.
Parameters
[in]N
          N is INTEGER
         On entry, the row dimension of the upper bidiagonal matrix.
         This is also the dimension of the main diagonal array D.
[in]SQRE
          SQRE is INTEGER
         Specifies the column dimension of the bidiagonal matrix.
         = 0: The bidiagonal matrix has column dimension M = N;
         = 1: The bidiagonal matrix has column dimension M = N+1;
[in,out]D
          D is DOUBLE PRECISION array, dimension (N)
         On entry D contains the main diagonal of the bidiagonal
         matrix.
         On exit D, if INFO = 0, contains its singular values.
[in]E
          E is DOUBLE PRECISION array, dimension (M-1)
         Contains the subdiagonal entries of the bidiagonal matrix.
         On exit, E has been destroyed.
[out]U
          U is DOUBLE PRECISION array, dimension at least (LDQ, N)
         On exit, U contains the left singular vectors.
[in]LDU
          LDU is INTEGER
         On entry, leading dimension of U.
[out]VT
          VT is DOUBLE PRECISION array, dimension at least (LDVT, M)
         On exit, VT**T contains the right singular vectors.
[in]LDVT
          LDVT is INTEGER
         On entry, leading dimension of VT.
[in]SMLSIZ
          SMLSIZ is INTEGER
         On entry, maximum size of the subproblems at the
         bottom of the computation tree.
[out]IWORK
          IWORK is INTEGER work array.
         Dimension must be at least (8 * N)
[out]WORK
          WORK is DOUBLE PRECISION work array.
         Dimension must be at least (3 * M**2 + 2 * M)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  if INFO = 1, a singular value did not converge
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015
Contributors:
Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 154 of file dlasd0.f.

154 *
155 * -- LAPACK auxiliary routine (version 3.6.0) --
156 * -- LAPACK is a software package provided by Univ. of Tennessee, --
157 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
158 * November 2015
159 *
160 * .. Scalar Arguments ..
161  INTEGER info, ldu, ldvt, n, smlsiz, sqre
162 * ..
163 * .. Array Arguments ..
164  INTEGER iwork( * )
165  DOUBLE PRECISION d( * ), e( * ), u( ldu, * ), vt( ldvt, * ),
166  $ work( * )
167 * ..
168 *
169 * =====================================================================
170 *
171 * .. Local Scalars ..
172  INTEGER i, i1, ic, idxq, idxqc, im1, inode, itemp, iwk,
173  $ j, lf, ll, lvl, m, ncc, nd, ndb1, ndiml, ndimr,
174  $ nl, nlf, nlp1, nlvl, nr, nrf, nrp1, sqrei
175  DOUBLE PRECISION alpha, beta
176 * ..
177 * .. External Subroutines ..
178  EXTERNAL dlasd1, dlasdq, dlasdt, xerbla
179 * ..
180 * .. Executable Statements ..
181 *
182 * Test the input parameters.
183 *
184  info = 0
185 *
186  IF( n.LT.0 ) THEN
187  info = -1
188  ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
189  info = -2
190  END IF
191 *
192  m = n + sqre
193 *
194  IF( ldu.LT.n ) THEN
195  info = -6
196  ELSE IF( ldvt.LT.m ) THEN
197  info = -8
198  ELSE IF( smlsiz.LT.3 ) THEN
199  info = -9
200  END IF
201  IF( info.NE.0 ) THEN
202  CALL xerbla( 'DLASD0', -info )
203  RETURN
204  END IF
205 *
206 * If the input matrix is too small, call DLASDQ to find the SVD.
207 *
208  IF( n.LE.smlsiz ) THEN
209  CALL dlasdq( 'U', sqre, n, m, n, 0, d, e, vt, ldvt, u, ldu, u,
210  $ ldu, work, info )
211  RETURN
212  END IF
213 *
214 * Set up the computation tree.
215 *
216  inode = 1
217  ndiml = inode + n
218  ndimr = ndiml + n
219  idxq = ndimr + n
220  iwk = idxq + n
221  CALL dlasdt( n, nlvl, nd, iwork( inode ), iwork( ndiml ),
222  $ iwork( ndimr ), smlsiz )
223 *
224 * For the nodes on bottom level of the tree, solve
225 * their subproblems by DLASDQ.
226 *
227  ndb1 = ( nd+1 ) / 2
228  ncc = 0
229  DO 30 i = ndb1, nd
230 *
231 * IC : center row of each node
232 * NL : number of rows of left subproblem
233 * NR : number of rows of right subproblem
234 * NLF: starting row of the left subproblem
235 * NRF: starting row of the right subproblem
236 *
237  i1 = i - 1
238  ic = iwork( inode+i1 )
239  nl = iwork( ndiml+i1 )
240  nlp1 = nl + 1
241  nr = iwork( ndimr+i1 )
242  nrp1 = nr + 1
243  nlf = ic - nl
244  nrf = ic + 1
245  sqrei = 1
246  CALL dlasdq( 'U', sqrei, nl, nlp1, nl, ncc, d( nlf ), e( nlf ),
247  $ vt( nlf, nlf ), ldvt, u( nlf, nlf ), ldu,
248  $ u( nlf, nlf ), ldu, work, info )
249  IF( info.NE.0 ) THEN
250  RETURN
251  END IF
252  itemp = idxq + nlf - 2
253  DO 10 j = 1, nl
254  iwork( itemp+j ) = j
255  10 CONTINUE
256  IF( i.EQ.nd ) THEN
257  sqrei = sqre
258  ELSE
259  sqrei = 1
260  END IF
261  nrp1 = nr + sqrei
262  CALL dlasdq( 'U', sqrei, nr, nrp1, nr, ncc, d( nrf ), e( nrf ),
263  $ vt( nrf, nrf ), ldvt, u( nrf, nrf ), ldu,
264  $ u( nrf, nrf ), ldu, work, info )
265  IF( info.NE.0 ) THEN
266  RETURN
267  END IF
268  itemp = idxq + ic
269  DO 20 j = 1, nr
270  iwork( itemp+j-1 ) = j
271  20 CONTINUE
272  30 CONTINUE
273 *
274 * Now conquer each subproblem bottom-up.
275 *
276  DO 50 lvl = nlvl, 1, -1
277 *
278 * Find the first node LF and last node LL on the
279 * current level LVL.
280 *
281  IF( lvl.EQ.1 ) THEN
282  lf = 1
283  ll = 1
284  ELSE
285  lf = 2**( lvl-1 )
286  ll = 2*lf - 1
287  END IF
288  DO 40 i = lf, ll
289  im1 = i - 1
290  ic = iwork( inode+im1 )
291  nl = iwork( ndiml+im1 )
292  nr = iwork( ndimr+im1 )
293  nlf = ic - nl
294  IF( ( sqre.EQ.0 ) .AND. ( i.EQ.ll ) ) THEN
295  sqrei = sqre
296  ELSE
297  sqrei = 1
298  END IF
299  idxqc = idxq + nlf - 1
300  alpha = d( ic )
301  beta = e( ic )
302  CALL dlasd1( nl, nr, sqrei, d( nlf ), alpha, beta,
303  $ u( nlf, nlf ), ldu, vt( nlf, nlf ), ldvt,
304  $ iwork( idxqc ), iwork( iwk ), work, info )
305 *
306 * Report the possible convergence failure.
307 *
308  IF( info.NE.0 ) THEN
309  RETURN
310  END IF
311  40 CONTINUE
312  50 CONTINUE
313 *
314  RETURN
315 *
316 * End of DLASD0
317 *
subroutine dlasdq(UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e...
Definition: dlasdq.f:213
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlasdt(N, LVL, ND, INODE, NDIML, NDIMR, MSUB)
DLASDT creates a tree of subproblems for bidiagonal divide and conquer. Used by sbdsdc.
Definition: dlasdt.f:107
subroutine dlasd1(NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT, IDXQ, IWORK, WORK, INFO)
DLASD1 computes the SVD of an upper bidiagonal matrix B of the specified size. Used by sbdsdc...
Definition: dlasd1.f:206

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dlasd1 ( integer  NL,
integer  NR,
integer  SQRE,
double precision, dimension( * )  D,
double precision  ALPHA,