## table of contents

auxOTHERcomputational(3) | LAPACK | auxOTHERcomputational(3) |

# NAME¶

auxOTHERcomputational - auxiliary Computational routines

# SYNOPSIS¶

## Functions¶

character *1 function **chla_transtype** (TRANS)

**CHLA_TRANSTYPE** subroutine **dbdsdc** (UPLO, COMPQ, N, D, E, U, LDU,
VT, LDVT, Q, IQ, WORK, IWORK, INFO)

**DBDSDC** subroutine **dbdsqr** (UPLO, N, NCVT, NRU, NCC, D, E, VT,
LDVT, U, LDU, C, LDC, WORK, INFO)

**DBDSQR** subroutine **ddisna** (JOB, M, N, D, SEP, INFO)

**DDISNA** subroutine **dlaed0** (ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE,
LDQS, WORK, IWORK, INFO)

**DLAED0** used by DSTEDC. Computes all eigenvalues and corresponding
eigenvectors of an unreduced symmetric tridiagonal matrix using the divide
and conquer method. subroutine **dlaed1** (N, D, Q, LDQ, INDXQ, RHO,
CUTPNT, WORK, IWORK, INFO)

**DLAED1** used by DSTEDC. Computes the updated eigensystem of a diagonal
matrix after modification by a rank-one symmetric matrix. Used when the
original matrix is tridiagonal. subroutine **dlaed2** (K, N, N1, D, Q,
LDQ, INDXQ, RHO, Z, DLAMDA, W, Q2, INDX, INDXC, INDXP, COLTYP, INFO)

**DLAED2** used by DSTEDC. Merges eigenvalues and deflates secular
equation. Used when the original matrix is tridiagonal. subroutine
**dlaed3** (K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX, CTOT, W, S, INFO)

**DLAED3** used by DSTEDC. Finds the roots of the secular equation and
updates the eigenvectors. Used when the original matrix is tridiagonal.
subroutine **dlaed4** (N, I, D, Z, DELTA, RHO, DLAM, INFO)

**DLAED4** used by DSTEDC. Finds a single root of the secular equation.
subroutine **dlaed5** (I, D, Z, DELTA, RHO, DLAM)

**DLAED5** used by DSTEDC. Solves the 2-by-2 secular equation. subroutine
**dlaed6** (KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO)

**DLAED6** used by DSTEDC. Computes one Newton step in solution of the
secular equation. subroutine **dlaed7** (ICOMPQ, N, QSIZ, TLVLS, CURLVL,
CURPBM, D, Q, LDQ, INDXQ, RHO, CUTPNT, QSTORE, QPTR, PRMPTR, PERM, GIVPTR,
GIVCOL, GIVNUM, WORK, IWORK, INFO)

**DLAED7** used by DSTEDC. Computes the updated eigensystem of a diagonal
matrix after modification by a rank-one symmetric matrix. Used when the
original matrix is dense. subroutine **dlaed8** (ICOMPQ, K, N, QSIZ, D,
Q, LDQ, INDXQ, RHO, CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR, GIVCOL,
GIVNUM, INDXP, INDX, INFO)

**DLAED8** used by DSTEDC. Merges eigenvalues and deflates secular
equation. Used when the original matrix is dense. subroutine **dlaed9**
(K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W, S, LDS, INFO)

**DLAED9** used by DSTEDC. Finds the roots of the secular equation and
updates the eigenvectors. Used when the original matrix is dense. subroutine
**dlaeda** (N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR, GIVCOL,
GIVNUM, Q, QPTR, Z, ZTEMP, INFO)

**DLAEDA** used by DSTEDC. Computes the Z vector determining the rank-one
modification of the diagonal matrix. Used when the original matrix is dense.
subroutine **dlagtf** (N, A, LAMBDA, B, C, TOL, D, IN, INFO)

**DLAGTF** computes an LU factorization of a matrix T-λI, where T is
a general tridiagonal matrix, and λ a scalar, using partial pivoting
with row interchanges. subroutine **dlamrg** (N1, N2, A, DTRD1, DTRD2,
INDEX)

**DLAMRG** creates a permutation list to merge the entries of two
independently sorted sets into a single set sorted in ascending order.
subroutine **dlartgs** (X, Y, SIGMA, CS, SN)

**DLARTGS** generates a plane rotation designed to introduce a bulge in
implicit QR iteration for the bidiagonal SVD problem. subroutine
**dlasq1** (N, D, E, WORK, INFO)

**DLASQ1** computes the singular values of a real square bidiagonal matrix.
Used by sbdsqr. subroutine **dlasq2** (N, Z, INFO)

**DLASQ2** computes all the eigenvalues of the symmetric positive definite
tridiagonal matrix associated with the qd Array Z to high relative accuracy.
Used by sbdsqr and sstegr. subroutine **dlasq3** (I0, N0, Z, PP, DMIN,
SIGMA, DESIG, QMAX, NFAIL, ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
DN2, G, TAU)

**DLASQ3** checks for deflation, computes a shift and calls dqds. Used by
sbdsqr. subroutine **dlasq4** (I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2,
DN, DN1, DN2, TAU, TTYPE, G)

**DLASQ4** computes an approximation to the smallest eigenvalue using
values of d from the previous transform. Used by sbdsqr. subroutine
**dlasq5** (I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN, DNM1,
DNM2, IEEE, EPS)

**DLASQ5** computes one dqds transform in ping-pong form. Used by sbdsqr
and sstegr. subroutine **dlasq6** (I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN,
DNM1, DNM2)

**DLASQ6** computes one dqd transform in ping-pong form. Used by sbdsqr and
sstegr. subroutine **dlasrt** (ID, N, D, INFO)

**DLASRT** sorts numbers in increasing or decreasing order. subroutine
**dstebz** (RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W,
IBLOCK, ISPLIT, WORK, IWORK, INFO)

**DSTEBZ** subroutine **dstedc** (COMPZ, N, D, E, Z, LDZ, WORK, LWORK,
IWORK, LIWORK, INFO)

**DSTEDC** subroutine **dsteqr** (COMPZ, N, D, E, Z, LDZ, WORK, INFO)

**DSTEQR** subroutine **dsterf** (N, D, E, INFO)

**DSTERF** integer function **iladiag** (DIAG)

**ILADIAG** integer function **ilaprec** (PREC)

**ILAPREC** integer function **ilatrans** (TRANS)

**ILATRANS** integer function **ilauplo** (UPLO)

**ILAUPLO** subroutine **sbdsdc** (UPLO, COMPQ, N, D, E, U, LDU, VT,
LDVT, Q, IQ, WORK, IWORK, INFO)

**SBDSDC** subroutine **sbdsqr** (UPLO, N, NCVT, NRU, NCC, D, E, VT,
LDVT, U, LDU, C, LDC, WORK, INFO)

**SBDSQR** subroutine **sdisna** (JOB, M, N, D, SEP, INFO)

**SDISNA** subroutine **slaed0** (ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE,
LDQS, WORK, IWORK, INFO)

**SLAED0** used by SSTEDC. Computes all eigenvalues and corresponding
eigenvectors of an unreduced symmetric tridiagonal matrix using the divide
and conquer method. subroutine **slaed1** (N, D, Q, LDQ, INDXQ, RHO,
CUTPNT, WORK, IWORK, INFO)

**SLAED1** used by SSTEDC. Computes the updated eigensystem of a diagonal
matrix after modification by a rank-one symmetric matrix. Used when the
original matrix is tridiagonal. subroutine **slaed2** (K, N, N1, D, Q,
LDQ, INDXQ, RHO, Z, DLAMDA, W, Q2, INDX, INDXC, INDXP, COLTYP, INFO)

**SLAED2** used by SSTEDC. Merges eigenvalues and deflates secular
equation. Used when the original matrix is tridiagonal. subroutine
**slaed3** (K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX, CTOT, W, S, INFO)

**SLAED3** used by SSTEDC. Finds the roots of the secular equation and
updates the eigenvectors. Used when the original matrix is tridiagonal.
subroutine **slaed4** (N, I, D, Z, DELTA, RHO, DLAM, INFO)

**SLAED4** used by SSTEDC. Finds a single root of the secular equation.
subroutine **slaed5** (I, D, Z, DELTA, RHO, DLAM)

**SLAED5** used by SSTEDC. Solves the 2-by-2 secular equation. subroutine
**slaed6** (KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO)

**SLAED6** used by SSTEDC. Computes one Newton step in solution of the
secular equation. subroutine **slaed7** (ICOMPQ, N, QSIZ, TLVLS, CURLVL,
CURPBM, D, Q, LDQ, INDXQ, RHO, CUTPNT, QSTORE, QPTR, PRMPTR, PERM, GIVPTR,
GIVCOL, GIVNUM, WORK, IWORK, INFO)

**SLAED7** used by SSTEDC. Computes the updated eigensystem of a diagonal
matrix after modification by a rank-one symmetric matrix. Used when the
original matrix is dense. subroutine **slaed8** (ICOMPQ, K, N, QSIZ, D,
Q, LDQ, INDXQ, RHO, CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR, GIVCOL,
GIVNUM, INDXP, INDX, INFO)

**SLAED8** used by SSTEDC. Merges eigenvalues and deflates secular
equation. Used when the original matrix is dense. subroutine **slaed9**
(K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W, S, LDS, INFO)

**SLAED9** used by SSTEDC. Finds the roots of the secular equation and
updates the eigenvectors. Used when the original matrix is dense. subroutine
**slaeda** (N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR, GIVCOL,
GIVNUM, Q, QPTR, Z, ZTEMP, INFO)

**SLAEDA** used by SSTEDC. Computes the Z vector determining the rank-one
modification of the diagonal matrix. Used when the original matrix is dense.
subroutine **slagtf** (N, A, LAMBDA, B, C, TOL, D, IN, INFO)

**SLAGTF** computes an LU factorization of a matrix T-λI, where T is
a general tridiagonal matrix, and λ a scalar, using partial pivoting
with row interchanges. subroutine **slamrg** (N1, N2, A, STRD1, STRD2,
INDEX)

**SLAMRG** creates a permutation list to merge the entries of two
independently sorted sets into a single set sorted in ascending order.
subroutine **slartgs** (X, Y, SIGMA, CS, SN)

**SLARTGS** generates a plane rotation designed to introduce a bulge in
implicit QR iteration for the bidiagonal SVD problem. subroutine
**slasq1** (N, D, E, WORK, INFO)

**SLASQ1** computes the singular values of a real square bidiagonal matrix.
Used by sbdsqr. subroutine **slasq2** (N, Z, INFO)

**SLASQ2** computes all the eigenvalues of the symmetric positive definite
tridiagonal matrix associated with the qd Array Z to high relative accuracy.
Used by sbdsqr and sstegr. subroutine **slasq3** (I0, N0, Z, PP, DMIN,
SIGMA, DESIG, QMAX, NFAIL, ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
DN2, G, TAU)

**SLASQ3** checks for deflation, computes a shift and calls dqds. Used by
sbdsqr. subroutine **slasq4** (I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2,
DN, DN1, DN2, TAU, TTYPE, G)

**SLASQ4** computes an approximation to the smallest eigenvalue using
values of d from the previous transform. Used by sbdsqr. subroutine
**slasq5** (I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN, DNM1,
DNM2, IEEE, EPS)

** SLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and
sstegr. ** subroutine **slasq6** (I0, N0, Z, PP, DMIN, DMIN1, DMIN2,
DN, DNM1, DNM2)

**SLASQ6** computes one dqd transform in ping-pong form. Used by sbdsqr and
sstegr. subroutine **slasrt** (ID, N, D, INFO)

**SLASRT** sorts numbers in increasing or decreasing order. subroutine
**spttrf** (N, D, E, INFO)

**SPTTRF** subroutine **sstebz** (RANGE, ORDER, N, VL, VU, IL, IU,
ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)

**SSTEBZ** subroutine **sstedc** (COMPZ, N, D, E, Z, LDZ, WORK, LWORK,
IWORK, LIWORK, INFO)

**SSTEDC** subroutine **ssteqr** (COMPZ, N, D, E, Z, LDZ, WORK, INFO)

**SSTEQR** subroutine **ssterf** (N, D, E, INFO)

**SSTERF**

# Detailed Description¶

This is the group of auxiliary Computational routines

# Function Documentation¶

## character*1 function chla_transtype (integer TRANS)¶

**CHLA_TRANSTYPE**

**Purpose:**

This subroutine translates from a BLAST-specified integer constant to

the character string specifying a transposition operation.

CHLA_TRANSTYPE returns an CHARACTER*1. If CHLA_TRANSTYPE is 'X',

then input is not an integer indicating a transposition operator.

Otherwise CHLA_TRANSTYPE returns the constant value corresponding to

TRANS.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine dbdsdc (character UPLO, character COMPQ, integer N, double precision, dimension( * ) D, double precision, dimension( * ) E, double precision, dimension( ldu, * ) U, integer LDU, double precision, dimension( ldvt, * ) VT, integer LDVT, double precision, dimension( * ) Q, integer, dimension( * ) IQ, double precision, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)¶

**DBDSDC**

**Purpose:**

DBDSDC computes the singular value decomposition (SVD) of a real

N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT,

using a divide and conquer method, where S is a diagonal matrix

with non-negative diagonal elements (the singular values of B), and

U and VT are orthogonal matrices of left and right singular vectors,

respectively. DBDSDC can be used to compute all singular values,

and optionally, singular vectors or singular vectors in compact form.

This code makes very mild assumptions about floating point

arithmetic. It will work on machines with a guard digit in

add/subtract, or on those binary machines without guard digits

which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.

It could conceivably fail on hexadecimal or decimal machines

without guard digits, but we know of none. See DLASD3 for details.

The code currently calls DLASDQ if singular values only are desired.

However, it can be slightly modified to compute singular values

using the divide and conquer method.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

= 'U': B is upper bidiagonal.

= 'L': B is lower bidiagonal.

*COMPQ*

COMPQ is CHARACTER*1

Specifies whether singular vectors are to be computed

as follows:

= 'N': Compute singular values only;

= 'P': Compute singular values and compute singular

vectors in compact form;

= 'I': Compute singular values and singular vectors.

*N*

N is INTEGER

The order of the matrix B. N >= 0.

*D*

D is DOUBLE PRECISION array, dimension (N)

On entry, the n diagonal elements of the bidiagonal matrix B.

On exit, if INFO=0, the singular values of B.

*E*

E is DOUBLE PRECISION array, dimension (N-1)

On entry, the elements of E contain the offdiagonal

elements of the bidiagonal matrix whose SVD is desired.

On exit, E has been destroyed.

*U*

U is DOUBLE PRECISION array, dimension (LDU,N)

If COMPQ = 'I', then:

On exit, if INFO = 0, U contains the left singular vectors

of the bidiagonal matrix.

For other values of COMPQ, U is not referenced.

*LDU*

LDU is INTEGER

The leading dimension of the array U. LDU >= 1.

If singular vectors are desired, then LDU >= max( 1, N ).

*VT*

VT is DOUBLE PRECISION array, dimension (LDVT,N)

If COMPQ = 'I', then:

On exit, if INFO = 0, VT**T contains the right singular

vectors of the bidiagonal matrix.

For other values of COMPQ, VT is not referenced.

*LDVT*

LDVT is INTEGER

The leading dimension of the array VT. LDVT >= 1.

If singular vectors are desired, then LDVT >= max( 1, N ).

*Q*

Q is DOUBLE PRECISION array, dimension (LDQ)

If COMPQ = 'P', then:

On exit, if INFO = 0, Q and IQ contain the left

and right singular vectors in a compact form,

requiring O(N log N) space instead of 2*N**2.

In particular, Q contains all the DOUBLE PRECISION data in

LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))

words of memory, where SMLSIZ is returned by ILAENV and

is equal to the maximum size of the subproblems at the

bottom of the computation tree (usually about 25).

For other values of COMPQ, Q is not referenced.

*IQ*

IQ is INTEGER array, dimension (LDIQ)

If COMPQ = 'P', then:

On exit, if INFO = 0, Q and IQ contain the left

and right singular vectors in a compact form,

requiring O(N log N) space instead of 2*N**2.

In particular, IQ contains all INTEGER data in

LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))

words of memory, where SMLSIZ is returned by ILAENV and

is equal to the maximum size of the subproblems at the

bottom of the computation tree (usually about 25).

For other values of COMPQ, IQ is not referenced.

*WORK*

WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))

If COMPQ = 'N' then LWORK >= (4 * N).

If COMPQ = 'P' then LWORK >= (6 * N).

If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).

*IWORK*

IWORK is INTEGER array, dimension (8*N)

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

> 0: The algorithm failed to compute a singular value.

The update process of divide and conquer failed.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine dbdsqr (character UPLO, integer N, integer NCVT, integer NRU, integer NCC, double precision, dimension( * ) D, double precision, dimension( * ) E, double precision, dimension( ldvt, * ) VT, integer LDVT, double precision, dimension( ldu, * ) U, integer LDU, double precision, dimension( ldc, * ) C, integer LDC, double precision, dimension( * ) WORK, integer INFO)¶

**DBDSQR**

**Purpose:**

DBDSQR computes the singular values and, optionally, the right and/or

left singular vectors from the singular value decomposition (SVD) of

a real N-by-N (upper or lower) bidiagonal matrix B using the implicit

zero-shift QR algorithm. The SVD of B has the form

B = Q * S * P**T

where S is the diagonal matrix of singular values, Q is an orthogonal

matrix of left singular vectors, and P is an orthogonal matrix of

right singular vectors. If left singular vectors are requested, this

subroutine actually returns U*Q instead of Q, and, if right singular

vectors are requested, this subroutine returns P**T*VT instead of

P**T, for given real input matrices U and VT. When U and VT are the

orthogonal matrices that reduce a general matrix A to bidiagonal

form: A = U*B*VT, as computed by DGEBRD, then

A = (U*Q) * S * (P**T*VT)

is the SVD of A. Optionally, the subroutine may also compute Q**T*C

for a given real input matrix C.

See "Computing Small Singular Values of Bidiagonal Matrices With

Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,

LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,

no. 5, pp. 873-912, Sept 1990) and

"Accurate singular values and differential qd algorithms," by

B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics

Department, University of California at Berkeley, July 1992

for a detailed description of the algorithm.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

= 'U': B is upper bidiagonal;

= 'L': B is lower bidiagonal.

*N*

N is INTEGER

The order of the matrix B. N >= 0.

*NCVT*

NCVT is INTEGER

The number of columns of the matrix VT. NCVT >= 0.

*NRU*

NRU is INTEGER

The number of rows of the matrix U. NRU >= 0.

*NCC*

NCC is INTEGER

The number of columns of the matrix C. NCC >= 0.

*D*

D is DOUBLE PRECISION array, dimension (N)

On entry, the n diagonal elements of the bidiagonal matrix B.

On exit, if INFO=0, the singular values of B in decreasing

order.

*E*

E is DOUBLE PRECISION array, dimension (N-1)

On entry, the N-1 offdiagonal elements of the bidiagonal

matrix B.

On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E

will contain the diagonal and superdiagonal elements of a

bidiagonal matrix orthogonally equivalent to the one given

as input.

*VT*

VT is DOUBLE PRECISION array, dimension (LDVT, NCVT)

On entry, an N-by-NCVT matrix VT.

On exit, VT is overwritten by P**T * VT.

Not referenced if NCVT = 0.

*LDVT*

LDVT is INTEGER

The leading dimension of the array VT.

LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.

*U*

U is DOUBLE PRECISION array, dimension (LDU, N)

On entry, an NRU-by-N matrix U.

On exit, U is overwritten by U * Q.

Not referenced if NRU = 0.

*LDU*

LDU is INTEGER

The leading dimension of the array U. LDU >= max(1,NRU).

*C*

C is DOUBLE PRECISION array, dimension (LDC, NCC)

On entry, an N-by-NCC matrix C.

On exit, C is overwritten by Q**T * C.

Not referenced if NCC = 0.

*LDC*

LDC is INTEGER

The leading dimension of the array C.

LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.

*WORK*

WORK is DOUBLE PRECISION array, dimension (4*(N-1))

*INFO*

INFO is INTEGER

= 0: successful exit

< 0: If INFO = -i, the i-th argument had an illegal value

> 0:

if NCVT = NRU = NCC = 0,

= 1, a split was marked by a positive value in E

= 2, current block of Z not diagonalized after 30*N

iterations (in inner while loop)

= 3, termination criterion of outer while loop not met

(program created more than N unreduced blocks)

else NCVT = NRU = NCC = 0,

the algorithm did not converge; D and E contain the

elements of a bidiagonal matrix which is orthogonally

similar to the input matrix B; if INFO = i, i

elements of E have not converged to zero.

**Internal Parameters:**

TOLMUL DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))

TOLMUL controls the convergence criterion of the QR loop.

If it is positive, TOLMUL*EPS is the desired relative

precision in the computed singular values.

If it is negative, abs(TOLMUL*EPS*sigma_max) is the

desired absolute accuracy in the computed singular

values (corresponds to relative accuracy

abs(TOLMUL*EPS) in the largest singular value.

abs(TOLMUL) should be between 1 and 1/EPS, and preferably

between 10 (for fast convergence) and .1/EPS

(for there to be some accuracy in the results).

Default is to lose at either one eighth or 2 of the

available decimal digits in each computed singular value

(whichever is smaller).

MAXITR INTEGER, default = 6

MAXITR controls the maximum number of passes of the

algorithm through its inner loop. The algorithms stops

(and so fails to converge) if the number of passes

through the inner loop exceeds MAXITR*N**2.

**Note:**

Bug report from Cezary Dendek.

On March 23rd 2017, the INTEGER variable MAXIT = MAXITR*N**2 is

removed since it can overflow pretty easily (for N larger or equal

than 18,919). We instead use MAXITDIVN = MAXITR*N.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine ddisna (character JOB, integer M, integer N, double precision, dimension( * ) D, double precision, dimension( * ) SEP, integer INFO)¶

**DDISNA**

**Purpose:**

DDISNA computes the reciprocal condition numbers for the eigenvectors

of a real symmetric or complex Hermitian matrix or for the left or

right singular vectors of a general m-by-n matrix. The reciprocal

condition number is the 'gap' between the corresponding eigenvalue or

singular value and the nearest other one.

The bound on the error, measured by angle in radians, in the I-th

computed vector is given by

DLAMCH( 'E' ) * ( ANORM / SEP( I ) )

where ANORM = 2-norm(A) = max( abs( D(j) ) ). SEP(I) is not allowed

to be smaller than DLAMCH( 'E' )*ANORM in order to limit the size of

the error bound.

DDISNA may also be used to compute error bounds for eigenvectors of

the generalized symmetric definite eigenproblem.

**Parameters**

*JOB*

JOB is CHARACTER*1

Specifies for which problem the reciprocal condition numbers

should be computed:

= 'E': the eigenvectors of a symmetric/Hermitian matrix;

= 'L': the left singular vectors of a general matrix;

= 'R': the right singular vectors of a general matrix.

*M*

M is INTEGER

The number of rows of the matrix. M >= 0.

*N*

N is INTEGER

If JOB = 'L' or 'R', the number of columns of the matrix,

in which case N >= 0. Ignored if JOB = 'E'.

*D*

D is DOUBLE PRECISION array, dimension (M) if JOB = 'E'

dimension (min(M,N)) if JOB = 'L' or 'R'

The eigenvalues (if JOB = 'E') or singular values (if JOB =

'L' or 'R') of the matrix, in either increasing or decreasing

order. If singular values, they must be non-negative.

*SEP*

SEP is DOUBLE PRECISION array, dimension (M) if JOB = 'E'

dimension (min(M,N)) if JOB = 'L' or 'R'

The reciprocal condition numbers of the vectors.

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine dlaed0 (integer ICOMPQ, integer QSIZ, integer N, double precision, dimension( * ) D, double precision, dimension( * ) E, double precision, dimension( ldq, * ) Q, integer LDQ, double precision, dimension( ldqs, * ) QSTORE, integer LDQS, double precision, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)¶

**DLAED0** used by DSTEDC. Computes all eigenvalues and
corresponding eigenvectors of an unreduced symmetric tridiagonal matrix
using the divide and conquer method.

**Purpose:**

DLAED0 computes all eigenvalues and corresponding eigenvectors of a

symmetric tridiagonal matrix using the divide and conquer method.

**Parameters**

*ICOMPQ*

ICOMPQ is INTEGER

= 0: Compute eigenvalues only.

= 1: Compute eigenvectors of original dense symmetric matrix

also. On entry, Q contains the orthogonal matrix used

to reduce the original matrix to tridiagonal form.

= 2: Compute eigenvalues and eigenvectors of tridiagonal

matrix.

*QSIZ*

QSIZ is INTEGER

The dimension of the orthogonal matrix used to reduce

the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.

*N*

N is INTEGER

The dimension of the symmetric tridiagonal matrix. N >= 0.

*D*

D is DOUBLE PRECISION array, dimension (N)

On entry, the main diagonal of the tridiagonal matrix.

On exit, its eigenvalues.

*E*

E is DOUBLE PRECISION array, dimension (N-1)

The off-diagonal elements of the tridiagonal matrix.

On exit, E has been destroyed.

*Q*

Q is DOUBLE PRECISION array, dimension (LDQ, N)

On entry, Q must contain an N-by-N orthogonal matrix.

If ICOMPQ = 0 Q is not referenced.

If ICOMPQ = 1 On entry, Q is a subset of the columns of the

orthogonal matrix used to reduce the full

matrix to tridiagonal form corresponding to

the subset of the full matrix which is being

decomposed at this time.

If ICOMPQ = 2 On entry, Q will be the identity matrix.

On exit, Q contains the eigenvectors of the

tridiagonal matrix.

*LDQ*

LDQ is INTEGER

The leading dimension of the array Q. If eigenvectors are

desired, then LDQ >= max(1,N). In any case, LDQ >= 1.

*QSTORE*

QSTORE is DOUBLE PRECISION array, dimension (LDQS, N)

Referenced only when ICOMPQ = 1. Used to store parts of

the eigenvector matrix when the updating matrix multiplies

take place.

*LDQS*

LDQS is INTEGER

The leading dimension of the array QSTORE. If ICOMPQ = 1,

then LDQS >= max(1,N). In any case, LDQS >= 1.

*WORK*

WORK is DOUBLE PRECISION array,

If ICOMPQ = 0 or 1, the dimension of WORK must be at least

1 + 3*N + 2*N*lg N + 3*N**2

( lg( N ) = smallest integer k

such that 2^k >= N )

If ICOMPQ = 2, the dimension of WORK must be at least

4*N + N**2.

*IWORK*

IWORK is INTEGER array,

If ICOMPQ = 0 or 1, the dimension of IWORK must be at least

6 + 6*N + 5*N*lg N.

( lg( N ) = smallest integer k

such that 2^k >= N )

If ICOMPQ = 2, the dimension of IWORK must be at least

3 + 5*N.

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

> 0: The algorithm failed to compute an eigenvalue while

working on the submatrix lying in rows and columns

INFO/(N+1) through mod(INFO,N+1).

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine dlaed1 (integer N, double precision, dimension( * ) D, double precision, dimension( ldq, * ) Q, integer LDQ, integer, dimension( * ) INDXQ, double precision RHO, integer CUTPNT, double precision, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)¶

**DLAED1** used by DSTEDC. Computes the updated eigensystem of
a diagonal matrix after modification by a rank-one symmetric matrix. Used
when the original matrix is tridiagonal.

**Purpose:**

DLAED1 computes the updated eigensystem of a diagonal

matrix after modification by a rank-one symmetric matrix. This

routine is used only for the eigenproblem which requires all

eigenvalues and eigenvectors of a tridiagonal matrix. DLAED7 handles

the case in which eigenvalues only or eigenvalues and eigenvectors

of a full symmetric matrix (which was reduced to tridiagonal form)

are desired.

T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)

where Z = Q**T*u, u is a vector of length N with ones in the

CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.

The eigenvectors of the original matrix are stored in Q, and the

eigenvalues are in D. The algorithm consists of three stages:

The first stage consists of deflating the size of the problem

when there are multiple eigenvalues or if there is a zero in

the Z vector. For each such occurrence the dimension of the

secular equation problem is reduced by one. This stage is

performed by the routine DLAED2.

The second stage consists of calculating the updated

eigenvalues. This is done by finding the roots of the secular

equation via the routine DLAED4 (as called by DLAED3).

This routine also calculates the eigenvectors of the current

problem.

The final stage consists of computing the updated eigenvectors

directly using the updated eigenvalues. The eigenvectors for

the current problem are multiplied with the eigenvectors from

the overall problem.

**Parameters**

*N*

N is INTEGER

The dimension of the symmetric tridiagonal matrix. N >= 0.

*D*

D is DOUBLE PRECISION array, dimension (N)

On entry, the eigenvalues of the rank-1-perturbed matrix.

On exit, the eigenvalues of the repaired matrix.

*Q*

Q is DOUBLE PRECISION array, dimension (LDQ,N)

On entry, the eigenvectors of the rank-1-perturbed matrix.

On exit, the eigenvectors of the repaired tridiagonal matrix.

*LDQ*

LDQ is INTEGER

The leading dimension of the array Q. LDQ >= max(1,N).

*INDXQ*

INDXQ is INTEGER array, dimension (N)

On entry, the permutation which separately sorts the two

subproblems in D into ascending order.

On exit, the permutation which will reintegrate the

subproblems back into sorted order,

i.e. D( INDXQ( I = 1, N ) ) will be in ascending order.

*RHO*

RHO is DOUBLE PRECISION

The subdiagonal entry used to create the rank-1 modification.

*CUTPNT*

CUTPNT is INTEGER

The location of the last eigenvalue in the leading sub-matrix.

min(1,N) <= CUTPNT <= N/2.

*WORK*

WORK is DOUBLE PRECISION array, dimension (4*N + N**2)

*IWORK*

IWORK is INTEGER array, dimension (4*N)

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

> 0: if INFO = 1, an eigenvalue did not converge

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

Modified by Francoise Tisseur, University of Tennessee

## subroutine dlaed2 (integer K, integer N, integer N1, double precision, dimension( * ) D, double precision, dimension( ldq, * ) Q, integer LDQ, integer, dimension( * ) INDXQ, double precision RHO, double precision, dimension( * ) Z, double precision, dimension( * ) DLAMDA, double precision, dimension( * ) W, double precision, dimension( * ) Q2, integer, dimension( * ) INDX, integer, dimension( * ) INDXC, integer, dimension( * ) INDXP, integer, dimension( * ) COLTYP, integer INFO)¶

**DLAED2** used by DSTEDC. Merges eigenvalues and deflates
secular equation. Used when the original matrix is tridiagonal.

**Purpose:**

DLAED2 merges the two sets of eigenvalues together into a single

sorted set. Then it tries to deflate the size of the problem.

There are two ways in which deflation can occur: when two or more

eigenvalues are close together or if there is a tiny entry in the

Z vector. For each such occurrence the order of the related secular

equation problem is reduced by one.

**Parameters**

*K*

K is INTEGER

The number of non-deflated eigenvalues, and the order of the

related secular equation. 0 <= K <=N.

*N*

N is INTEGER

The dimension of the symmetric tridiagonal matrix. N >= 0.

*N1*

N1 is INTEGER

The location of the last eigenvalue in the leading sub-matrix.

min(1,N) <= N1 <= N/2.

*D*

D is DOUBLE PRECISION array, dimension (N)

On entry, D contains the eigenvalues of the two submatrices to

be combined.

On exit, D contains the trailing (N-K) updated eigenvalues

(those which were deflated) sorted into increasing order.

*Q*

Q is DOUBLE PRECISION array, dimension (LDQ, N)

On entry, Q contains the eigenvectors of two submatrices in

the two square blocks with corners at (1,1), (N1,N1)

and (N1+1, N1+1), (N,N).

On exit, Q contains the trailing (N-K) updated eigenvectors

(those which were deflated) in its last N-K columns.

*LDQ*

LDQ is INTEGER

The leading dimension of the array Q. LDQ >= max(1,N).

*INDXQ*

INDXQ is INTEGER array, dimension (N)

The permutation which separately sorts the two sub-problems

in D into ascending order. Note that elements in the second

half of this permutation must first have N1 added to their

values. Destroyed on exit.

*RHO*

RHO is DOUBLE PRECISION

On entry, the off-diagonal element associated with the rank-1

cut which originally split the two submatrices which are now

being recombined.

On exit, RHO has been modified to the value required by

DLAED3.

*Z*

Z is DOUBLE PRECISION array, dimension (N)

On entry, Z contains the updating vector (the last

row of the first sub-eigenvector matrix and the first row of

the second sub-eigenvector matrix).

On exit, the contents of Z have been destroyed by the updating

process.

*DLAMDA*

DLAMDA is DOUBLE PRECISION array, dimension (N)

A copy of the first K eigenvalues which will be used by

DLAED3 to form the secular equation.

*W*

W is DOUBLE PRECISION array, dimension (N)

The first k values of the final deflation-altered z-vector

which will be passed to DLAED3.

*Q2*

Q2 is DOUBLE PRECISION array, dimension (N1**2+(N-N1)**2)

A copy of the first K eigenvectors which will be used by

DLAED3 in a matrix multiply (DGEMM) to solve for the new

eigenvectors.

*INDX*

INDX is INTEGER array, dimension (N)

The permutation used to sort the contents of DLAMDA into

ascending order.

*INDXC*

INDXC is INTEGER array, dimension (N)

The permutation used to arrange the columns of the deflated

Q matrix into three groups: the first group contains non-zero

elements only at and above N1, the second contains

non-zero elements only below N1, and the third is dense.

*INDXP*

INDXP is INTEGER array, dimension (N)

The permutation used to place deflated values of D at the end

of the array. INDXP(1:K) points to the nondeflated D-values

and INDXP(K+1:N) points to the deflated eigenvalues.

*COLTYP*

COLTYP is INTEGER array, dimension (N)

During execution, a label which will indicate which of the

following types a column in the Q2 matrix is:

1 : non-zero in the upper half only;

2 : dense;

3 : non-zero in the lower half only;

4 : deflated.

On exit, COLTYP(i) is the number of columns of type i,

for i=1 to 4 only.

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

Modified by Francoise Tisseur, University of Tennessee

## subroutine dlaed3 (integer K, integer N, integer N1, double precision, dimension( * ) D, double precision, dimension( ldq, * ) Q, integer LDQ, double precision RHO, double precision, dimension( * ) DLAMDA, double precision, dimension( * ) Q2, integer, dimension( * ) INDX, integer, dimension( * ) CTOT, double precision, dimension( * ) W, double precision, dimension( * ) S, integer INFO)¶

**DLAED3** used by DSTEDC. Finds the roots of the secular
equation and updates the eigenvectors. Used when the original matrix is
tridiagonal.

**Purpose:**

DLAED3 finds the roots of the secular equation, as defined by the

values in D, W, and RHO, between 1 and K. It makes the

appropriate calls to DLAED4 and then updates the eigenvectors by

multiplying the matrix of eigenvectors of the pair of eigensystems

being combined by the matrix of eigenvectors of the K-by-K system

which is solved here.

This code makes very mild assumptions about floating point

arithmetic. It will work on machines with a guard digit in

add/subtract, or on those binary machines without guard digits

which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.

It could conceivably fail on hexadecimal or decimal machines

without guard digits, but we know of none.

**Parameters**

*K*

K is INTEGER

The number of terms in the rational function to be solved by

DLAED4. K >= 0.

*N*

N is INTEGER

The number of rows and columns in the Q matrix.

N >= K (deflation may result in N>K).

*N1*

N1 is INTEGER

The location of the last eigenvalue in the leading submatrix.

min(1,N) <= N1 <= N/2.

*D*

D is DOUBLE PRECISION array, dimension (N)

D(I) contains the updated eigenvalues for

1 <= I <= K.

*Q*

Q is DOUBLE PRECISION array, dimension (LDQ,N)

Initially the first K columns are used as workspace.

On output the columns 1 to K contain

the updated eigenvectors.

*LDQ*

LDQ is INTEGER

The leading dimension of the array Q. LDQ >= max(1,N).

*RHO*

RHO is DOUBLE PRECISION

The value of the parameter in the rank one update equation.

RHO >= 0 required.

*DLAMDA*

DLAMDA is DOUBLE PRECISION array, dimension (K)

The first K elements of this array contain the old roots

of the deflated updating problem. These are the poles

of the secular equation. May be changed on output by

having lowest order bit set to zero on Cray X-MP, Cray Y-MP,

Cray-2, or Cray C-90, as described above.

*Q2*

Q2 is DOUBLE PRECISION array, dimension (LDQ2*N)

The first K columns of this matrix contain the non-deflated

eigenvectors for the split problem.

*INDX*

INDX is INTEGER array, dimension (N)

The permutation used to arrange the columns of the deflated

Q matrix into three groups (see DLAED2).

The rows of the eigenvectors found by DLAED4 must be likewise

permuted before the matrix multiply can take place.

*CTOT*

CTOT is INTEGER array, dimension (4)

A count of the total number of the various types of columns

in Q, as described in INDX. The fourth column type is any

column which has been deflated.

*W*

W is DOUBLE PRECISION array, dimension (K)

The first K elements of this array contain the components

of the deflation-adjusted updating vector. Destroyed on

output.

*S*

S is DOUBLE PRECISION array, dimension (N1 + 1)*K

Will contain the eigenvectors of the repaired matrix which

will be multiplied by the previously accumulated eigenvectors

to update the system.

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

> 0: if INFO = 1, an eigenvalue did not converge

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

Modified by Francoise Tisseur, University of Tennessee

## subroutine dlaed4 (integer N, integer I, double precision, dimension( * ) D, double precision, dimension( * ) Z, double precision, dimension( * ) DELTA, double precision RHO, double precision DLAM, integer INFO)¶

**DLAED4** used by DSTEDC. Finds a single root of the secular
equation.

**Purpose:**

This subroutine computes the I-th updated eigenvalue of a symmetric

rank-one modification to a diagonal matrix whose elements are

given in the array d, and that

D(i) < D(j) for i < j

and that RHO > 0. This is arranged by the calling routine, and is

no loss in generality. The rank-one modified system is thus

diag( D ) + RHO * Z * Z_transpose.

where we assume the Euclidean norm of Z is 1.

The method consists of approximating the rational functions in the

secular equation by simpler interpolating rational functions.

**Parameters**

*N*

N is INTEGER

The length of all arrays.

*I*

I is INTEGER

The index of the eigenvalue to be computed. 1 <= I <= N.

*D*

D is DOUBLE PRECISION array, dimension (N)

The original eigenvalues. It is assumed that they are in

order, D(I) < D(J) for I < J.

*Z*

Z is DOUBLE PRECISION array, dimension (N)

The components of the updating vector.

*DELTA*

DELTA is DOUBLE PRECISION array, dimension (N)

If N > 2, DELTA contains (D(j) - lambda_I) in its j-th

component. If N = 1, then DELTA(1) = 1. If N = 2, see DLAED5

for detail. The vector DELTA contains the information necessary

to construct the eigenvectors by DLAED3 and DLAED9.

*RHO*

RHO is DOUBLE PRECISION

The scalar in the symmetric updating formula.

*DLAM*

DLAM is DOUBLE PRECISION

The computed lambda_I, the I-th updated eigenvalue.

*INFO*

INFO is INTEGER

= 0: successful exit

> 0: if INFO = 1, the updating process failed.

**Internal Parameters:**

Logical variable ORGATI (origin-at-i?) is used for distinguishing

whether D(i) or D(i+1) is treated as the origin.

ORGATI = .true. origin at i

ORGATI = .false. origin at i+1

Logical variable SWTCH3 (switch-for-3-poles?) is for noting

if we are working with THREE poles!

MAXIT is the maximum number of iterations allowed for each

eigenvalue.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine dlaed5 (integer I, double precision, dimension( 2 ) D, double precision, dimension( 2 ) Z, double precision, dimension( 2 ) DELTA, double precision RHO, double precision DLAM)¶

**DLAED5** used by DSTEDC. Solves the 2-by-2 secular
equation.

**Purpose:**

This subroutine computes the I-th eigenvalue of a symmetric rank-one

modification of a 2-by-2 diagonal matrix

diag( D ) + RHO * Z * transpose(Z) .

The diagonal elements in the array D are assumed to satisfy

D(i) < D(j) for i < j .

We also assume RHO > 0 and that the Euclidean norm of the vector

Z is one.

**Parameters**

*I*

I is INTEGER

The index of the eigenvalue to be computed. I = 1 or I = 2.

*D*

D is DOUBLE PRECISION array, dimension (2)

The original eigenvalues. We assume D(1) < D(2).

*Z*

Z is DOUBLE PRECISION array, dimension (2)

The components of the updating vector.

*DELTA*

DELTA is DOUBLE PRECISION array, dimension (2)

The vector DELTA contains the information necessary

to construct the eigenvectors.

*RHO*

RHO is DOUBLE PRECISION

The scalar in the symmetric updating formula.

*DLAM*

DLAM is DOUBLE PRECISION

The computed lambda_I, the I-th updated eigenvalue.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine dlaed6 (integer KNITER, logical ORGATI, double precision RHO, double precision, dimension( 3 ) D, double precision, dimension( 3 ) Z, double precision FINIT, double precision TAU, integer INFO)¶

**DLAED6** used by DSTEDC. Computes one Newton step in solution
of the secular equation.

**Purpose:**

DLAED6 computes the positive or negative root (closest to the origin)

of

z(1) z(2) z(3)

f(x) = rho + --------- + ---------- + ---------

d(1)-x d(2)-x d(3)-x

It is assumed that

if ORGATI = .true. the root is between d(2) and d(3);

otherwise it is between d(1) and d(2)

This routine will be called by DLAED4 when necessary. In most cases,

the root sought is the smallest in magnitude, though it might not be

in some extremely rare situations.

**Parameters**

*KNITER*

KNITER is INTEGER

Refer to DLAED4 for its significance.

*ORGATI*

ORGATI is LOGICAL

If ORGATI is true, the needed root is between d(2) and

d(3); otherwise it is between d(1) and d(2). See

DLAED4 for further details.

*RHO*

RHO is DOUBLE PRECISION

Refer to the equation f(x) above.

*D*

D is DOUBLE PRECISION array, dimension (3)

D satisfies d(1) < d(2) < d(3).

*Z*

Z is DOUBLE PRECISION array, dimension (3)

Each of the elements in z must be positive.

*FINIT*

FINIT is DOUBLE PRECISION

The value of f at 0. It is more accurate than the one

evaluated inside this routine (if someone wants to do

so).

*TAU*

TAU is DOUBLE PRECISION

The root of the equation f(x).

*INFO*

INFO is INTEGER

= 0: successful exit

> 0: if INFO = 1, failure to converge

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

10/02/03: This version has a few statements commented out for thread

safety (machine parameters are computed on each entry). SJH.

05/10/06: Modified from a new version of Ren-Cang Li, use

Gragg-Thornton-Warner cubic convergent scheme for better stability.

**Contributors:**

## subroutine dlaed7 (integer ICOMPQ, integer N, integer QSIZ, integer TLVLS, integer CURLVL, integer CURPBM, double precision, dimension( * ) D, double precision, dimension( ldq, * ) Q, integer LDQ, integer, dimension( * ) INDXQ, double precision RHO, integer CUTPNT, double precision, dimension( * ) QSTORE, integer, dimension( * ) QPTR, integer, dimension( * ) PRMPTR, integer, dimension( * ) PERM, integer, dimension( * ) GIVPTR, integer, dimension( 2, * ) GIVCOL, double precision, dimension( 2, * ) GIVNUM, double precision, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)¶

**DLAED7** used by DSTEDC. Computes the updated eigensystem of
a diagonal matrix after modification by a rank-one symmetric matrix. Used
when the original matrix is dense.

**Purpose:**

DLAED7 computes the updated eigensystem of a diagonal

matrix after modification by a rank-one symmetric matrix. This

routine is used only for the eigenproblem which requires all

eigenvalues and optionally eigenvectors of a dense symmetric matrix

that has been reduced to tridiagonal form. DLAED1 handles

the case in which all eigenvalues and eigenvectors of a symmetric

tridiagonal matrix are desired.

T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)

where Z = Q**Tu, u is a vector of length N with ones in the

CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.

The eigenvectors of the original matrix are stored in Q, and the

eigenvalues are in D. The algorithm consists of three stages:

The first stage consists of deflating the size of the problem

when there are multiple eigenvalues or if there is a zero in

the Z vector. For each such occurrence the dimension of the

secular equation problem is reduced by one. This stage is

performed by the routine DLAED8.

The second stage consists of calculating the updated

eigenvalues. This is done by finding the roots of the secular

equation via the routine DLAED4 (as called by DLAED9).

This routine also calculates the eigenvectors of the current

problem.

The final stage consists of computing the updated eigenvectors

directly using the updated eigenvalues. The eigenvectors for

the current problem are multiplied with the eigenvectors from

the overall problem.

**Parameters**

*ICOMPQ*

ICOMPQ is INTEGER

= 0: Compute eigenvalues only.

= 1: Compute eigenvectors of original dense symmetric matrix

also. On entry, Q contains the orthogonal matrix used

to reduce the original matrix to tridiagonal form.

*N*

N is INTEGER

The dimension of the symmetric tridiagonal matrix. N >= 0.

*QSIZ*

QSIZ is INTEGER

The dimension of the orthogonal matrix used to reduce

the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.

*TLVLS*

TLVLS is INTEGER

The total number of merging levels in the overall divide and

conquer tree.

*CURLVL*

CURLVL is INTEGER

The current level in the overall merge routine,

0 <= CURLVL <= TLVLS.

*CURPBM*

CURPBM is INTEGER

The current problem in the current level in the overall

merge routine (counting from upper left to lower right).

*D*

D is DOUBLE PRECISION array, dimension (N)

On entry, the eigenvalues of the rank-1-perturbed matrix.

On exit, the eigenvalues of the repaired matrix.

*Q*

Q is DOUBLE PRECISION array, dimension (LDQ, N)

On entry, the eigenvectors of the rank-1-perturbed matrix.

On exit, the eigenvectors of the repaired tridiagonal matrix.

*LDQ*

LDQ is INTEGER

The leading dimension of the array Q. LDQ >= max(1,N).

*INDXQ*

INDXQ is INTEGER array, dimension (N)

The permutation which will reintegrate the subproblem just

solved back into sorted order, i.e., D( INDXQ( I = 1, N ) )

will be in ascending order.

*RHO*

RHO is DOUBLE PRECISION

The subdiagonal element used to create the rank-1

modification.

*CUTPNT*

CUTPNT is INTEGER

Contains the location of the last eigenvalue in the leading

sub-matrix. min(1,N) <= CUTPNT <= N.

*QSTORE*

QSTORE is DOUBLE PRECISION array, dimension (N**2+1)

Stores eigenvectors of submatrices encountered during

divide and conquer, packed together. QPTR points to

beginning of the submatrices.

*QPTR*

QPTR is INTEGER array, dimension (N+2)

List of indices pointing to beginning of submatrices stored

in QSTORE. The submatrices are numbered starting at the

bottom left of the divide and conquer tree, from left to

right and bottom to top.

*PRMPTR*

PRMPTR is INTEGER array, dimension (N lg N)

Contains a list of pointers which indicate where in PERM a

level's permutation is stored. PRMPTR(i+1) - PRMPTR(i)

indicates the size of the permutation and also the size of

the full, non-deflated problem.

*PERM*

PERM is INTEGER array, dimension (N lg N)

Contains the permutations (from deflation and sorting) to be

applied to each eigenblock.

*GIVPTR*

GIVPTR is INTEGER array, dimension (N lg N)

Contains a list of pointers which indicate where in GIVCOL a

level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i)

indicates the number of Givens rotations.

*GIVCOL*

GIVCOL is INTEGER array, dimension (2, N lg N)

Each pair of numbers indicates a pair of columns to take place

in a Givens rotation.

*GIVNUM*

GIVNUM is DOUBLE PRECISION array, dimension (2, N lg N)

Each number indicates the S value to be used in the

corresponding Givens rotation.

*WORK*

WORK is DOUBLE PRECISION array, dimension (3*N+2*QSIZ*N)

*IWORK*

IWORK is INTEGER array, dimension (4*N)

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

> 0: if INFO = 1, an eigenvalue did not converge

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine dlaed8 (integer ICOMPQ, integer K, integer N, integer QSIZ, double precision, dimension( * ) D, double precision, dimension( ldq, * ) Q, integer LDQ, integer, dimension( * ) INDXQ, double precision RHO, integer CUTPNT, double precision, dimension( * ) Z, double precision, dimension( * ) DLAMDA, double precision, dimension( ldq2, * ) Q2, integer LDQ2, double precision, dimension( * ) W, integer, dimension( * ) PERM, integer GIVPTR, integer, dimension( 2, * ) GIVCOL, double precision, dimension( 2, * ) GIVNUM, integer, dimension( * ) INDXP, integer, dimension( * ) INDX, integer INFO)¶

**DLAED8** used by DSTEDC. Merges eigenvalues and deflates
secular equation. Used when the original matrix is dense.

**Purpose:**

DLAED8 merges the two sets of eigenvalues together into a single

sorted set. Then it tries to deflate the size of the problem.

There are two ways in which deflation can occur: when two or more

eigenvalues are close together or if there is a tiny element in the

Z vector. For each such occurrence the order of the related secular

equation problem is reduced by one.

**Parameters**

*ICOMPQ*

ICOMPQ is INTEGER

= 0: Compute eigenvalues only.

= 1: Compute eigenvectors of original dense symmetric matrix

also. On entry, Q contains the orthogonal matrix used

to reduce the original matrix to tridiagonal form.

*K*

K is INTEGER

The number of non-deflated eigenvalues, and the order of the

related secular equation.

*N*

N is INTEGER

The dimension of the symmetric tridiagonal matrix. N >= 0.

*QSIZ*

QSIZ is INTEGER

The dimension of the orthogonal matrix used to reduce

the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.

*D*

D is DOUBLE PRECISION array, dimension (N)

On entry, the eigenvalues of the two submatrices to be

combined. On exit, the trailing (N-K) updated eigenvalues

(those which were deflated) sorted into increasing order.

*Q*

Q is DOUBLE PRECISION array, dimension (LDQ,N)

If ICOMPQ = 0, Q is not referenced. Otherwise,

on entry, Q contains the eigenvectors of the partially solved

system which has been previously updated in matrix

multiplies with other partially solved eigensystems.

On exit, Q contains the trailing (N-K) updated eigenvectors

(those which were deflated) in its last N-K columns.

*LDQ*

LDQ is INTEGER

The leading dimension of the array Q. LDQ >= max(1,N).

*INDXQ*

INDXQ is INTEGER array, dimension (N)

The permutation which separately sorts the two sub-problems

in D into ascending order. Note that elements in the second

half of this permutation must first have CUTPNT added to

their values in order to be accurate.

*RHO*

RHO is DOUBLE PRECISION

On entry, the off-diagonal element associated with the rank-1

cut which originally split the two submatrices which are now

being recombined.

On exit, RHO has been modified to the value required by

DLAED3.

*CUTPNT*

CUTPNT is INTEGER

The location of the last eigenvalue in the leading

sub-matrix. min(1,N) <= CUTPNT <= N.

*Z*

Z is DOUBLE PRECISION array, dimension (N)

On entry, Z contains the updating vector (the last row of

the first sub-eigenvector matrix and the first row of the

second sub-eigenvector matrix).

On exit, the contents of Z are destroyed by the updating

process.

*DLAMDA*

DLAMDA is DOUBLE PRECISION array, dimension (N)

A copy of the first K eigenvalues which will be used by

DLAED3 to form the secular equation.

*Q2*

Q2 is DOUBLE PRECISION array, dimension (LDQ2,N)

If ICOMPQ = 0, Q2 is not referenced. Otherwise,

a copy of the first K eigenvectors which will be used by

DLAED7 in a matrix multiply (DGEMM) to update the new

eigenvectors.

*LDQ2*

LDQ2 is INTEGER

The leading dimension of the array Q2. LDQ2 >= max(1,N).

*W*

W is DOUBLE PRECISION array, dimension (N)

The first k values of the final deflation-altered z-vector and

will be passed to DLAED3.

*PERM*

PERM is INTEGER array, dimension (N)

The permutations (from deflation and sorting) to be applied

to each eigenblock.

*GIVPTR*

GIVPTR is INTEGER

The number of Givens rotations which took place in this

subproblem.

*GIVCOL*

GIVCOL is INTEGER array, dimension (2, N)

Each pair of numbers indicates a pair of columns to take place

in a Givens rotation.

*GIVNUM*

GIVNUM is DOUBLE PRECISION array, dimension (2, N)

Each number indicates the S value to be used in the

corresponding Givens rotation.

*INDXP*

INDXP is INTEGER array, dimension (N)

The permutation used to place deflated values of D at the end

of the array. INDXP(1:K) points to the nondeflated D-values

and INDXP(K+1:N) points to the deflated eigenvalues.

*INDX*

INDX is INTEGER array, dimension (N)

The permutation used to sort the contents of D into ascending

order.

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine dlaed9 (integer K, integer KSTART, integer KSTOP, integer N, double precision, dimension( * ) D, double precision, dimension( ldq, * ) Q, integer LDQ, double precision RHO, double precision, dimension( * ) DLAMDA, double precision, dimension( * ) W, double precision, dimension( lds, * ) S, integer LDS, integer INFO)¶

**DLAED9** used by DSTEDC. Finds the roots of the secular
equation and updates the eigenvectors. Used when the original matrix is
dense.

**Purpose:**

DLAED9 finds the roots of the secular equation, as defined by the

values in D, Z, and RHO, between KSTART and KSTOP. It makes the

appropriate calls to DLAED4 and then stores the new matrix of

eigenvectors for use in calculating the next level of Z vectors.

**Parameters**

*K*

K is INTEGER

The number of terms in the rational function to be solved by

DLAED4. K >= 0.

*KSTART*

KSTART is INTEGER

*KSTOP*

KSTOP is INTEGER

The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP

are to be computed. 1 <= KSTART <= KSTOP <= K.

*N*

N is INTEGER

The number of rows and columns in the Q matrix.

N >= K (delation may result in N > K).

*D*

D is DOUBLE PRECISION array, dimension (N)

D(I) contains the updated eigenvalues

for KSTART <= I <= KSTOP.

*Q*

Q is DOUBLE PRECISION array, dimension (LDQ,N)

*LDQ*

LDQ is INTEGER

The leading dimension of the array Q. LDQ >= max( 1, N ).

*RHO*

RHO is DOUBLE PRECISION

The value of the parameter in the rank one update equation.

RHO >= 0 required.

*DLAMDA*

DLAMDA is DOUBLE PRECISION array, dimension (K)

The first K elements of this array contain the old roots

of the deflated updating problem. These are the poles

of the secular equation.

*W*

W is DOUBLE PRECISION array, dimension (K)

The first K elements of this array contain the components

of the deflation-adjusted updating vector.

*S*

S is DOUBLE PRECISION array, dimension (LDS, K)

Will contain the eigenvectors of the repaired matrix which

will be stored for subsequent Z vector calculation and

multiplied by the previously accumulated eigenvectors

to update the system.

*LDS*

LDS is INTEGER

The leading dimension of S. LDS >= max( 1, K ).

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

> 0: if INFO = 1, an eigenvalue did not converge

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine dlaeda (integer N, integer TLVLS, integer CURLVL, integer CURPBM, integer, dimension( * ) PRMPTR, integer, dimension( * ) PERM, integer, dimension( * ) GIVPTR, integer, dimension( 2, * ) GIVCOL, double precision, dimension( 2, * ) GIVNUM, double precision, dimension( * ) Q, integer, dimension( * ) QPTR, double precision, dimension( * ) Z, double precision, dimension( * ) ZTEMP, integer INFO)¶

**DLAEDA** used by DSTEDC. Computes the Z vector determining
the rank-one modification of the diagonal matrix. Used when the original
matrix is dense.

**Purpose:**

DLAEDA computes the Z vector corresponding to the merge step in the

CURLVLth step of the merge process with TLVLS steps for the CURPBMth

problem.

**Parameters**

*N*

N is INTEGER

The dimension of the symmetric tridiagonal matrix. N >= 0.

*TLVLS*

TLVLS is INTEGER

The total number of merging levels in the overall divide and

conquer tree.

*CURLVL*

CURLVL is INTEGER

The current level in the overall merge routine,

0 <= curlvl <= tlvls.

*CURPBM*

CURPBM is INTEGER

The current problem in the current level in the overall

merge routine (counting from upper left to lower right).

*PRMPTR*

PRMPTR is INTEGER array, dimension (N lg N)

Contains a list of pointers which indicate where in PERM a

level's permutation is stored. PRMPTR(i+1) - PRMPTR(i)

indicates the size of the permutation and incidentally the

size of the full, non-deflated problem.

*PERM*

PERM is INTEGER array, dimension (N lg N)

Contains the permutations (from deflation and sorting) to be

applied to each eigenblock.

*GIVPTR*

GIVPTR is INTEGER array, dimension (N lg N)

Contains a list of pointers which indicate where in GIVCOL a

level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i)

indicates the number of Givens rotations.

*GIVCOL*

GIVCOL is INTEGER array, dimension (2, N lg N)

Each pair of numbers indicates a pair of columns to take place

in a Givens rotation.

*GIVNUM*

GIVNUM is DOUBLE PRECISION array, dimension (2, N lg N)

Each number indicates the S value to be used in the

corresponding Givens rotation.

*Q*

Q is DOUBLE PRECISION array, dimension (N**2)

Contains the square eigenblocks from previous levels, the

starting positions for blocks are given by QPTR.

*QPTR*

QPTR is INTEGER array, dimension (N+2)

Contains a list of pointers which indicate where in Q an

eigenblock is stored. SQRT( QPTR(i+1) - QPTR(i) ) indicates

the size of the block.

*Z*

Z is DOUBLE PRECISION array, dimension (N)

On output this vector contains the updating vector (the last

row of the first sub-eigenvector matrix and the first row of

the second sub-eigenvector matrix).

*ZTEMP*

ZTEMP is DOUBLE PRECISION array, dimension (N)

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine dlagtf (integer N, double precision, dimension( * ) A, double precision LAMBDA, double precision, dimension( * ) B, double precision, dimension( * ) C, double precision TOL, double precision, dimension( * ) D, integer, dimension( * ) IN, integer INFO)¶

**DLAGTF** computes an LU factorization of a matrix
T-λI, where T is a general tridiagonal matrix, and λ a scalar,
using partial pivoting with row interchanges.

**Purpose:**

DLAGTF factorizes the matrix (T - lambda*I), where T is an n by n

tridiagonal matrix and lambda is a scalar, as

T - lambda*I = PLU,

where P is a permutation matrix, L is a unit lower tridiagonal matrix

with at most one non-zero sub-diagonal elements per column and U is

an upper triangular matrix with at most two non-zero super-diagonal

elements per column.

The factorization is obtained by Gaussian elimination with partial

pivoting and implicit row scaling.

The parameter LAMBDA is included in the routine so that DLAGTF may

be used, in conjunction with DLAGTS, to obtain eigenvectors of T by

inverse iteration.

**Parameters**

*N*

N is INTEGER

The order of the matrix T.

*A*

A is DOUBLE PRECISION array, dimension (N)

On entry, A must contain the diagonal elements of T.

On exit, A is overwritten by the n diagonal elements of the

upper triangular matrix U of the factorization of T.

*LAMBDA*

LAMBDA is DOUBLE PRECISION

On entry, the scalar lambda.

*B*

B is DOUBLE PRECISION array, dimension (N-1)

On entry, B must contain the (n-1) super-diagonal elements of

T.

On exit, B is overwritten by the (n-1) super-diagonal

elements of the matrix U of the factorization of T.

*C*

C is DOUBLE PRECISION array, dimension (N-1)

On entry, C must contain the (n-1) sub-diagonal elements of

T.

On exit, C is overwritten by the (n-1) sub-diagonal elements

of the matrix L of the factorization of T.

*TOL*

TOL is DOUBLE PRECISION

On entry, a relative tolerance used to indicate whether or

not the matrix (T - lambda*I) is nearly singular. TOL should

normally be chose as approximately the largest relative error

in the elements of T. For example, if the elements of T are

correct to about 4 significant figures, then TOL should be

set to about 5*10**(-4). If TOL is supplied as less than eps,

where eps is the relative machine precision, then the value

eps is used in place of TOL.

*D*

D is DOUBLE PRECISION array, dimension (N-2)

On exit, D is overwritten by the (n-2) second super-diagonal

elements of the matrix U of the factorization of T.

*IN*

IN is INTEGER array, dimension (N)

On exit, IN contains details of the permutation matrix P. If

an interchange occurred at the kth step of the elimination,

then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)

returns the smallest positive integer j such that

abs( u(j,j) ) <= norm( (T - lambda*I)(j) )*TOL,

where norm( A(j) ) denotes the sum of the absolute values of

the jth row of the matrix A. If no such j exists then IN(n)

is returned as zero. If IN(n) is returned as positive, then a

diagonal element of U is small, indicating that

(T - lambda*I) is singular or nearly singular,

*INFO*

INFO is INTEGER

= 0: successful exit

< 0: if INFO = -k, the kth argument had an illegal value

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine dlamrg (integer N1, integer N2, double precision, dimension( * ) A, integer DTRD1, integer DTRD2, integer, dimension( * ) INDEX)¶

**DLAMRG** creates a permutation list to merge the entries of
two independently sorted sets into a single set sorted in ascending
order.

**Purpose:**

DLAMRG will create a permutation list which will merge the elements

of A (which is composed of two independently sorted sets) into a

single set which is sorted in ascending order.

**Parameters**

*N1*

N1 is INTEGER

*N2*

N2 is INTEGER

These arguments contain the respective lengths of the two

sorted lists to be merged.

*A*

A is DOUBLE PRECISION array, dimension (N1+N2)

The first N1 elements of A contain a list of numbers which

are sorted in either ascending or descending order. Likewise

for the final N2 elements.

*DTRD1*

DTRD1 is INTEGER

*DTRD2*

DTRD2 is INTEGER

These are the strides to be taken through the array A.

Allowable strides are 1 and -1. They indicate whether a

subset of A is sorted in ascending (DTRDx = 1) or descending

(DTRDx = -1) order.

*INDEX*

INDEX is INTEGER array, dimension (N1+N2)

On exit this array will contain a permutation such that

if B( I ) = A( INDEX( I ) ) for I=1,N1+N2, then B will be

sorted in ascending order.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine dlartgs (double precision X, double precision Y, double precision SIGMA, double precision CS, double precision SN)¶

**DLARTGS** generates a plane rotation designed to introduce a
bulge in implicit QR iteration for the bidiagonal SVD problem.

**Purpose:**

DLARTGS generates a plane rotation designed to introduce a bulge in

Golub-Reinsch-style implicit QR iteration for the bidiagonal SVD

problem. X and Y are the top-row entries, and SIGMA is the shift.

The computed CS and SN define a plane rotation satisfying

[ CS SN ] . [ X^2 - SIGMA ] = [ R ],

[ -SN CS ] [ X * Y ] [ 0 ]

with R nonnegative. If X^2 - SIGMA and X * Y are 0, then the

rotation is by PI/2.

**Parameters**

*X*

X is DOUBLE PRECISION

The (1,1) entry of an upper bidiagonal matrix.

*Y*

Y is DOUBLE PRECISION

The (1,2) entry of an upper bidiagonal matrix.

*SIGMA*

SIGMA is DOUBLE PRECISION

The shift.

*CS*

CS is DOUBLE PRECISION

The cosine of the rotation.

*SN*

SN is DOUBLE PRECISION

The sine of the rotation.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine dlasq1 (integer N, double precision, dimension( * ) D, double precision, dimension( * ) E, double precision, dimension( * ) WORK, integer INFO)¶

**DLASQ1** computes the singular values of a real square
bidiagonal matrix. Used by sbdsqr.

**Purpose:**

DLASQ1 computes the singular values of a real N-by-N bidiagonal

matrix with diagonal D and off-diagonal E. The singular values

are computed to high relative accuracy, in the absence of

denormalization, underflow and overflow. The algorithm was first

presented in

"Accurate singular values and differential qd algorithms" by K. V.

Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,

1994,

and the present implementation is described in "An implementation of

the dqds Algorithm (Positive Case)", LAPACK Working Note.

**Parameters**

*N*

N is INTEGER

The number of rows and columns in the matrix. N >= 0.

*D*

D is DOUBLE PRECISION array, dimension (N)

On entry, D contains the diagonal elements of the

bidiagonal matrix whose SVD is desired. On normal exit,

D contains the singular values in decreasing order.

*E*

E is DOUBLE PRECISION array, dimension (N)

On entry, elements E(1:N-1) contain the off-diagonal elements

of the bidiagonal matrix whose SVD is desired.

On exit, E is overwritten.

*WORK*

WORK is DOUBLE PRECISION array, dimension (4*N)

*INFO*

INFO is INTEGER

= 0: successful exit

< 0: if INFO = -i, the i-th argument had an illegal value

> 0: the algorithm failed

= 1, a split was marked by a positive value in E

= 2, current block of Z not diagonalized after 100*N

iterations (in inner while loop) On exit D and E

represent a matrix with the same singular values

which the calling subroutine could use to finish the

computation, or even feed back into DLASQ1

= 3, termination criterion of outer while loop not met

(program created more than N unreduced blocks)

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine dlasq2 (integer N, double precision, dimension( * ) Z, integer INFO)¶

**DLASQ2** computes all the eigenvalues of the symmetric
positive definite tridiagonal matrix associated with the qd Array Z to high
relative accuracy. Used by sbdsqr and sstegr.

**Purpose:**

DLASQ2 computes all the eigenvalues of the symmetric positive

definite tridiagonal matrix associated with the qd array Z to high

relative accuracy are computed to high relative accuracy, in the

absence of denormalization, underflow and overflow.

To see the relation of Z to the tridiagonal matrix, let L be a

unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and

let U be an upper bidiagonal matrix with 1's above and diagonal

Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the

symmetric tridiagonal to which it is similar.

Note : DLASQ2 defines a logical variable, IEEE, which is true

on machines which follow ieee-754 floating-point standard in their

handling of infinities and NaNs, and false otherwise. This variable

is passed to DLASQ3.

**Parameters**

*N*

N is INTEGER

The number of rows and columns in the matrix. N >= 0.

*Z*

Z is DOUBLE PRECISION array, dimension ( 4*N )

On entry Z holds the qd array. On exit, entries 1 to N hold

the eigenvalues in decreasing order, Z( 2*N+1 ) holds the

trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If

N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )

holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of

shifts that failed.

*INFO*

INFO is INTEGER

= 0: successful exit

< 0: if the i-th argument is a scalar and had an illegal

value, then INFO = -i, if the i-th argument is an

array and the j-entry had an illegal value, then

INFO = -(i*100+j)

> 0: the algorithm failed

= 1, a split was marked by a positive value in E

= 2, current block of Z not diagonalized after 100*N

iterations (in inner while loop). On exit Z holds

a qd array with the same eigenvalues as the given Z.

= 3, termination criterion of outer while loop not met

(program created more than N unreduced blocks)

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

Local Variables: I0:N0 defines a current unreduced segment of Z.

The shifts are accumulated in SIGMA. Iteration count is in ITER.

Ping-pong is controlled by PP (alternates between 0 and 1).

## subroutine dlasq3 (integer I0, integer N0, double precision, dimension( * ) Z, integer PP, double precision DMIN, double precision SIGMA, double precision DESIG, double precision QMAX, integer NFAIL, integer ITER, integer NDIV, logical IEEE, integer TTYPE, double precision DMIN1, double precision DMIN2, double precision DN, double precision DN1, double precision DN2, double precision G, double precision TAU)¶

**DLASQ3** checks for deflation, computes a shift and calls
dqds. Used by sbdsqr.

**Purpose:**

DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.

In case of failure it changes shifts, and tries again until output

is positive.

**Parameters**

*I0*

I0 is INTEGER

First index.

*N0*

N0 is INTEGER

Last index.

*Z*

Z is DOUBLE PRECISION array, dimension ( 4*N0 )

Z holds the qd array.

*PP*

PP is INTEGER

PP=0 for ping, PP=1 for pong.

PP=2 indicates that flipping was applied to the Z array

and that the initial tests for deflation should not be

performed.

*DMIN*

DMIN is DOUBLE PRECISION

Minimum value of d.

*SIGMA*

SIGMA is DOUBLE PRECISION

Sum of shifts used in current segment.

*DESIG*

DESIG is DOUBLE PRECISION

Lower order part of SIGMA

*QMAX*

QMAX is DOUBLE PRECISION

Maximum value of q.

*NFAIL*

NFAIL is INTEGER

Increment NFAIL by 1 each time the shift was too big.

*ITER*

ITER is INTEGER

Increment ITER by 1 for each iteration.

*NDIV*

NDIV is INTEGER

Increment NDIV by 1 for each division.

*IEEE*

IEEE is LOGICAL

Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).

*TTYPE*

TTYPE is INTEGER

Shift type.

*DMIN1*

DMIN1 is DOUBLE PRECISION

*DMIN2*

DMIN2 is DOUBLE PRECISION

*DN*

DN is DOUBLE PRECISION

*DN1*

DN1 is DOUBLE PRECISION

*DN2*

DN2 is DOUBLE PRECISION

*G*

G is DOUBLE PRECISION

*TAU*

TAU is DOUBLE PRECISION

These are passed as arguments in order to save their values

between calls to DLASQ3.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine dlasq4 (integer I0, integer N0, double precision, dimension( * ) Z, integer PP, integer N0IN, double precision DMIN, double precision DMIN1, double precision DMIN2, double precision DN, double precision DN1, double precision DN2, double precision TAU, integer TTYPE, double precision G)¶

**DLASQ4** computes an approximation to the smallest eigenvalue
using values of d from the previous transform. Used by sbdsqr.

**Purpose:**

DLASQ4 computes an approximation TAU to the smallest eigenvalue

using values of d from the previous transform.

**Parameters**

*I0*

I0 is INTEGER

First index.

*N0*

N0 is INTEGER

Last index.

*Z*

Z is DOUBLE PRECISION array, dimension ( 4*N0 )

Z holds the qd array.

*PP*

PP is INTEGER

PP=0 for ping, PP=1 for pong.

*N0IN*

N0IN is INTEGER

The value of N0 at start of EIGTEST.

*DMIN*

DMIN is DOUBLE PRECISION

Minimum value of d.

*DMIN1*

DMIN1 is DOUBLE PRECISION

Minimum value of d, excluding D( N0 ).

*DMIN2*

DMIN2 is DOUBLE PRECISION

Minimum value of d, excluding D( N0 ) and D( N0-1 ).

*DN*

DN is DOUBLE PRECISION

d(N)

*DN1*

DN1 is DOUBLE PRECISION

d(N-1)

*DN2*

DN2 is DOUBLE PRECISION

d(N-2)

*TAU*

TAU is DOUBLE PRECISION

This is the shift.

*TTYPE*

TTYPE is INTEGER

Shift type.

*G*

G is DOUBLE PRECISION

G is passed as an argument in order to save its value between

calls to DLASQ4.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

CNST1 = 9/16

## subroutine dlasq5 (integer I0, integer N0, double precision, dimension( * ) Z, integer PP, double precision TAU, double precision SIGMA, double precision DMIN, double precision DMIN1, double precision DMIN2, double precision DN, double precision DNM1, double precision DNM2, logical IEEE, double precision EPS)¶

**DLASQ5** computes one dqds transform in ping-pong form. Used
by sbdsqr and sstegr.

**Purpose:**

DLASQ5 computes one dqds transform in ping-pong form, one

version for IEEE machines another for non IEEE machines.

**Parameters**

*I0*

I0 is INTEGER

First index.

*N0*

N0 is INTEGER

Last index.

*Z*

Z is DOUBLE PRECISION array, dimension ( 4*N )

Z holds the qd array. EMIN is stored in Z(4*N0) to avoid

an extra argument.

*PP*

PP is INTEGER

PP=0 for ping, PP=1 for pong.

*TAU*

TAU is DOUBLE PRECISION

This is the shift.

*SIGMA*

SIGMA is DOUBLE PRECISION

This is the accumulated shift up to this step.

*DMIN*

DMIN is DOUBLE PRECISION

Minimum value of d.

*DMIN1*

DMIN1 is DOUBLE PRECISION

Minimum value of d, excluding D( N0 ).

*DMIN2*

DMIN2 is DOUBLE PRECISION

Minimum value of d, excluding D( N0 ) and D( N0-1 ).

*DN*

DN is DOUBLE PRECISION

d(N0), the last value of d.

*DNM1*

DNM1 is DOUBLE PRECISION

d(N0-1).

*DNM2*

DNM2 is DOUBLE PRECISION

d(N0-2).

*IEEE*

IEEE is LOGICAL

Flag for IEEE or non IEEE arithmetic.

*EPS*

EPS is DOUBLE PRECISION

This is the value of epsilon used.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine dlasq6 (integer I0, integer N0, double precision, dimension( * ) Z, integer PP, double precision DMIN, double precision DMIN1, double precision DMIN2, double precision DN, double precision DNM1, double precision DNM2)¶

**DLASQ6** computes one dqd transform in ping-pong form. Used
by sbdsqr and sstegr.

**Purpose:**

DLASQ6 computes one dqd (shift equal to zero) transform in

ping-pong form, with protection against underflow and overflow.

**Parameters**

*I0*

I0 is INTEGER

First index.

*N0*

N0 is INTEGER

Last index.

*Z*

Z is DOUBLE PRECISION array, dimension ( 4*N )

Z holds the qd array. EMIN is stored in Z(4*N0) to avoid

an extra argument.

*PP*

PP is INTEGER

PP=0 for ping, PP=1 for pong.

*DMIN*

DMIN is DOUBLE PRECISION

Minimum value of d.

*DMIN1*

DMIN1 is DOUBLE PRECISION

Minimum value of d, excluding D( N0 ).

*DMIN2*

DMIN2 is DOUBLE PRECISION

Minimum value of d, excluding D( N0 ) and D( N0-1 ).

*DN*

DN is DOUBLE PRECISION

d(N0), the last value of d.

*DNM1*

DNM1 is DOUBLE PRECISION

d(N0-1).

*DNM2*

DNM2 is DOUBLE PRECISION

d(N0-2).

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine dlasrt (character ID, integer N, double precision, dimension( * ) D, integer INFO)¶

**DLASRT** sorts numbers in increasing or decreasing order.

**Purpose:**

Sort the numbers in D in increasing order (if ID = 'I') or

in decreasing order (if ID = 'D' ).

Use Quick Sort, reverting to Insertion sort on arrays of

size <= 20. Dimension of STACK limits N to about 2**32.

**Parameters**

*ID*

ID is CHARACTER*1

= 'I': sort D in increasing order;

= 'D': sort D in decreasing order.

*N*

N is INTEGER

The length of the array D.

*D*

D is DOUBLE PRECISION array, dimension (N)

On entry, the array to be sorted.

On exit, D has been sorted into increasing order

(D(1) <= ... <= D(N) ) or into decreasing order

(D(1) >= ... >= D(N) ), depending on ID.

*INFO*

INFO is INTEGER

= 0: successful exit

< 0: if INFO = -i, the i-th argument had an illegal value

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine dstebz (character RANGE, character ORDER, integer N, double precision VL, double precision VU, integer IL, integer IU, double precision ABSTOL, double precision, dimension( * ) D, double precision, dimension( * ) E, integer M, integer NSPLIT, double precision, dimension( * ) W, integer, dimension( * ) IBLOCK, integer, dimension( * ) ISPLIT, double precision, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)¶

**DSTEBZ**

**Purpose:**

DSTEBZ computes the eigenvalues of a symmetric tridiagonal

matrix T. 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**

*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.

*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.

*N*

N is INTEGER

The order of the tridiagonal matrix T. N >= 0.

*VL*

VL is DOUBLE PRECISION

If RANGE='V', the lower bound 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'.

*VU*

VU is DOUBLE PRECISION

If RANGE='V', the upper bound 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'.

*IL*

IL is INTEGER

If RANGE='I', the index of the

smallest eigenvalue to be returned.

1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.

Not referenced if RANGE = 'A' or 'V'.

*IU*

IU is INTEGER

If RANGE='I', the index of the

largest eigenvalue to be returned.

1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.

Not referenced if RANGE = 'A' or 'V'.

*ABSTOL*

ABSTOL is DOUBLE PRECISION

The absolute tolerance for the eigenvalues. An eigenvalue

(or cluster) is considered to be located if it has been

determined to lie in an interval whose width is ABSTOL or

less. If ABSTOL is less than or equal to zero, then ULP*|T|

will be used, where |T| means the 1-norm of T.

Eigenvalues will be computed most accurately when ABSTOL is

set to twice the underflow threshold 2*DLAMCH('S'), not zero.

*D*

D is DOUBLE PRECISION array, dimension (N)

The n diagonal elements of the tridiagonal matrix T.

*E*

E is DOUBLE PRECISION array, dimension (N-1)

The (n-1) off-diagonal elements of the tridiagonal matrix T.

*M*

M is INTEGER

The actual number of eigenvalues found. 0 <= M <= N.

(See also the description of INFO=2,3.)

*NSPLIT*

NSPLIT is INTEGER

The number of diagonal blocks in the matrix T.

1 <= NSPLIT <= N.

*W*

W is DOUBLE PRECISION array, dimension (N)

On exit, the first M elements of W will contain the

eigenvalues. (DSTEBZ may use the remaining N-M elements as

workspace.)

*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. (DSTEBZ may use the remaining N-M elements as

workspace.)

*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.)

*WORK*

WORK is DOUBLE PRECISION array, dimension (4*N)

*IWORK*

IWORK is INTEGER array, dimension (3*N)

*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:**

RELFAC DOUBLE PRECISION, default = 2.0e0

The relative tolerance. An interval (a,b] lies within

"relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|),

where "ulp" is the machine precision (distance from 1 to

the next larger floating point number.)

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.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine dstedc (character COMPZ, integer N, double precision, dimension( * ) D, double precision, dimension( * ) E, double precision, dimension( ldz, * ) Z, integer LDZ, double precision, dimension( * ) WORK, integer LWORK, integer, dimension( * ) IWORK, integer LIWORK, integer INFO)¶

**DSTEDC**

**Purpose:**

DSTEDC computes all eigenvalues and, optionally, eigenvectors of a

symmetric tridiagonal matrix using the divide and conquer method.

The eigenvectors of a full or band real symmetric matrix can also be

found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this

matrix to tridiagonal form.

This code makes very mild assumptions about floating point

arithmetic. It will work on machines with a guard digit in

add/subtract, or on those binary machines without guard digits

which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.

It could conceivably fail on hexadecimal or decimal machines

without guard digits, but we know of none. See DLAED3 for details.

**Parameters**

*COMPZ*

COMPZ is CHARACTER*1

= 'N': Compute eigenvalues only.

= 'I': Compute eigenvectors of tridiagonal matrix also.

= 'V': Compute eigenvectors of original dense symmetric

matrix also. On entry, Z contains the orthogonal

matrix used to reduce the original matrix to

tridiagonal form.

*N*

N is INTEGER

The dimension of the symmetric tridiagonal matrix. N >= 0.

*D*

D is DOUBLE PRECISION array, dimension (N)

On entry, the diagonal elements of the tridiagonal matrix.

On exit, if INFO = 0, the eigenvalues in ascending order.

*E*

E is DOUBLE PRECISION array, dimension (N-1)

On entry, the subdiagonal elements of the tridiagonal matrix.

On exit, E has been destroyed.

*Z*

Z is DOUBLE PRECISION array, dimension (LDZ,N)

On entry, if COMPZ = 'V', then Z contains the orthogonal

matrix used in the reduction to tridiagonal form.

On exit, if INFO = 0, then if COMPZ = 'V', Z contains the

orthonormal eigenvectors of the original symmetric matrix,

and if COMPZ = 'I', Z contains the orthonormal eigenvectors

of the symmetric tridiagonal matrix.

If COMPZ = 'N', then Z is not referenced.

*LDZ*

LDZ is INTEGER

The leading dimension of the array Z. LDZ >= 1.

If eigenvectors are desired, then LDZ >= max(1,N).

*WORK*

WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))

On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

*LWORK*

LWORK is INTEGER

The dimension of the array WORK.

If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.

If COMPZ = 'V' and N > 1 then LWORK must be at least

( 1 + 3*N + 2*N*lg N + 4*N**2 ),

where lg( N ) = smallest integer k such

that 2**k >= N.

If COMPZ = 'I' and N > 1 then LWORK must be at least

( 1 + 4*N + N**2 ).

Note that for COMPZ = 'I' or 'V', then if N is less than or

equal to the minimum divide size, usually 25, then LWORK need

only be max(1,2*(N-1)).

If LWORK = -1, then a workspace query is assumed; the routine

only calculates the optimal size of the WORK array, returns

this value as the first entry of the WORK array, and no error

message related to LWORK is issued by XERBLA.

*IWORK*

IWORK is INTEGER array, dimension (MAX(1,LIWORK))

On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.

*LIWORK*

LIWORK is INTEGER

The dimension of the array IWORK.

If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.

If COMPZ = 'V' and N > 1 then LIWORK must be at least

( 6 + 6*N + 5*N*lg N ).

If COMPZ = 'I' and N > 1 then LIWORK must be at least

( 3 + 5*N ).

Note that for COMPZ = 'I' or 'V', then if N is less than or

equal to the minimum divide size, usually 25, then LIWORK

need only be 1.

If LIWORK = -1, then a workspace query is assumed; the

routine only calculates the optimal size of the IWORK array,

returns this value as the first entry of the IWORK array, and

no error message related to LIWORK is issued by XERBLA.

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

> 0: The algorithm failed to compute an eigenvalue while

working on the submatrix lying in rows and columns

INFO/(N+1) through mod(INFO,N+1).

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

Modified by Francoise Tisseur, University of Tennessee

## subroutine dsteqr (character COMPZ, integer N, double precision, dimension( * ) D, double precision, dimension( * ) E, double precision, dimension( ldz, * ) Z, integer LDZ, double precision, dimension( * ) WORK, integer INFO)¶

**DSTEQR**

**Purpose:**

DSTEQR computes all eigenvalues and, optionally, eigenvectors of a

symmetric tridiagonal matrix using the implicit QL or QR method.

The eigenvectors of a full or band symmetric matrix can also be found

if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to

tridiagonal form.

**Parameters**

*COMPZ*

COMPZ is CHARACTER*1

= 'N': Compute eigenvalues only.

= 'V': Compute eigenvalues and eigenvectors of the original

symmetric matrix. On entry, Z must contain the

orthogonal matrix used to reduce the original matrix

to tridiagonal form.

= 'I': Compute eigenvalues and eigenvectors of the

tridiagonal matrix. Z is initialized to the identity

matrix.

*N*

N is INTEGER

The order of the matrix. N >= 0.

*D*

D is DOUBLE PRECISION array, dimension (N)

On entry, the diagonal elements of the tridiagonal matrix.

On exit, if INFO = 0, the eigenvalues in ascending order.

*E*

E is DOUBLE PRECISION array, dimension (N-1)

On entry, the (n-1) subdiagonal elements of the tridiagonal

matrix.

On exit, E has been destroyed.

*Z*

Z is DOUBLE PRECISION array, dimension (LDZ, N)

On entry, if COMPZ = 'V', then Z contains the orthogonal

matrix used in the reduction to tridiagonal form.

On exit, if INFO = 0, then if COMPZ = 'V', Z contains the

orthonormal eigenvectors of the original symmetric matrix,

and if COMPZ = 'I', Z contains the orthonormal eigenvectors

of the symmetric tridiagonal matrix.

If COMPZ = 'N', then Z is not referenced.

*LDZ*

LDZ is INTEGER

The leading dimension of the array Z. LDZ >= 1, and if

eigenvectors are desired, then LDZ >= max(1,N).

*WORK*

WORK is DOUBLE PRECISION array, dimension (max(1,2*N-2))

If COMPZ = 'N', then WORK is not referenced.

*INFO*

INFO is INTEGER

= 0: successful exit

< 0: if INFO = -i, the i-th argument had an illegal value

> 0: the algorithm has failed to find all the eigenvalues in

a total of 30*N iterations; if INFO = i, then i

elements of E have not converged to zero; on exit, D

and E contain the elements of a symmetric tridiagonal

matrix which is orthogonally similar to the original

matrix.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine dsterf (integer N, double precision, dimension( * ) D, double precision, dimension( * ) E, integer INFO)¶

**DSTERF**

**Purpose:**

DSTERF computes all eigenvalues of a symmetric tridiagonal matrix

using the Pal-Walker-Kahan variant of the QL or QR algorithm.

**Parameters**

*N*

N is INTEGER

The order of the matrix. N >= 0.

*D*

D is DOUBLE PRECISION array, dimension (N)

On entry, the n diagonal elements of the tridiagonal matrix.

On exit, if INFO = 0, the eigenvalues in ascending order.

*E*

E is DOUBLE PRECISION array, dimension (N-1)

On entry, the (n-1) subdiagonal elements of the tridiagonal

matrix.

On exit, E has been destroyed.

*INFO*

INFO is INTEGER

= 0: successful exit

< 0: if INFO = -i, the i-th argument had an illegal value

> 0: the algorithm failed to find all of the eigenvalues in

a total of 30*N iterations; if INFO = i, then i

elements of E have not converged to zero.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## integer function iladiag (character DIAG)¶

**ILADIAG**

**Purpose:**

This subroutine translated from a character string specifying if a

matrix has unit diagonal or not to the relevant BLAST-specified

integer constant.

ILADIAG returns an INTEGER. If ILADIAG < 0, then the input is not a

character indicating a unit or non-unit diagonal. Otherwise ILADIAG

returns the constant value corresponding to DIAG.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## integer function ilaprec (character PREC)¶

**ILAPREC**

**Purpose:**

This subroutine translated from a character string specifying an

intermediate precision to the relevant BLAST-specified integer

constant.

ILAPREC returns an INTEGER. If ILAPREC < 0, then the input is not a

character indicating a supported intermediate precision. Otherwise

ILAPREC returns the constant value corresponding to PREC.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## integer function ilatrans (character TRANS)¶

**ILATRANS**

**Purpose:**

This subroutine translates from a character string specifying a

transposition operation to the relevant BLAST-specified integer

constant.

ILATRANS returns an INTEGER. If ILATRANS < 0, then the input is not

a character indicating a transposition operator. Otherwise ILATRANS

returns the constant value corresponding to TRANS.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## integer function ilauplo (character UPLO)¶

**ILAUPLO**

**Purpose:**

This subroutine translated from a character string specifying a

upper- or lower-triangular matrix to the relevant BLAST-specified

integer constant.

ILAUPLO returns an INTEGER. If ILAUPLO < 0, then the input is not

a character indicating an upper- or lower-triangular matrix.

Otherwise ILAUPLO returns the constant value corresponding to UPLO.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine sbdsdc (character UPLO, character COMPQ, integer N, real, dimension( * ) D, real, dimension( * ) E, real, dimension( ldu, * ) U, integer LDU, real, dimension( ldvt, * ) VT, integer LDVT, real, dimension( * ) Q, integer, dimension( * ) IQ, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)¶

**SBDSDC**

**Purpose:**

SBDSDC computes the singular value decomposition (SVD) of a real

N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT,

using a divide and conquer method, where S is a diagonal matrix

with non-negative diagonal elements (the singular values of B), and

U and VT are orthogonal matrices of left and right singular vectors,

respectively. SBDSDC can be used to compute all singular values,

and optionally, singular vectors or singular vectors in compact form.

This code makes very mild assumptions about floating point

arithmetic. It will work on machines with a guard digit in

add/subtract, or on those binary machines without guard digits

which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.

It could conceivably fail on hexadecimal or decimal machines

without guard digits, but we know of none. See SLASD3 for details.

The code currently calls SLASDQ if singular values only are desired.

However, it can be slightly modified to compute singular values

using the divide and conquer method.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

= 'U': B is upper bidiagonal.

= 'L': B is lower bidiagonal.

*COMPQ*

COMPQ is CHARACTER*1

Specifies whether singular vectors are to be computed

as follows:

= 'N': Compute singular values only;

= 'P': Compute singular values and compute singular

vectors in compact form;

= 'I': Compute singular values and singular vectors.

*N*

N is INTEGER

The order of the matrix B. N >= 0.

*D*

D is REAL array, dimension (N)

On entry, the n diagonal elements of the bidiagonal matrix B.

On exit, if INFO=0, the singular values of B.

*E*

E is REAL array, dimension (N-1)

On entry, the elements of E contain the offdiagonal

elements of the bidiagonal matrix whose SVD is desired.

On exit, E has been destroyed.

*U*

U is REAL array, dimension (LDU,N)

If COMPQ = 'I', then:

On exit, if INFO = 0, U contains the left singular vectors

of the bidiagonal matrix.

For other values of COMPQ, U is not referenced.

*LDU*

LDU is INTEGER

The leading dimension of the array U. LDU >= 1.

If singular vectors are desired, then LDU >= max( 1, N ).

*VT*

VT is REAL array, dimension (LDVT,N)

If COMPQ = 'I', then:

On exit, if INFO = 0, VT**T contains the right singular

vectors of the bidiagonal matrix.

For other values of COMPQ, VT is not referenced.

*LDVT*

LDVT is INTEGER

The leading dimension of the array VT. LDVT >= 1.

If singular vectors are desired, then LDVT >= max( 1, N ).

*Q*

Q is REAL array, dimension (LDQ)

If COMPQ = 'P', then:

On exit, if INFO = 0, Q and IQ contain the left

and right singular vectors in a compact form,

requiring O(N log N) space instead of 2*N**2.

In particular, Q contains all the REAL data in

LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))

words of memory, where SMLSIZ is returned by ILAENV and

is equal to the maximum size of the subproblems at the

bottom of the computation tree (usually about 25).

For other values of COMPQ, Q is not referenced.

*IQ*

IQ is INTEGER array, dimension (LDIQ)

If COMPQ = 'P', then:

On exit, if INFO = 0, Q and IQ contain the left

and right singular vectors in a compact form,

requiring O(N log N) space instead of 2*N**2.

In particular, IQ contains all INTEGER data in

LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))

words of memory, where SMLSIZ is returned by ILAENV and

is equal to the maximum size of the subproblems at the

bottom of the computation tree (usually about 25).

For other values of COMPQ, IQ is not referenced.

*WORK*

WORK is REAL array, dimension (MAX(1,LWORK))

If COMPQ = 'N' then LWORK >= (4 * N).

If COMPQ = 'P' then LWORK >= (6 * N).

If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).

*IWORK*

IWORK is INTEGER array, dimension (8*N)

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

> 0: The algorithm failed to compute a singular value.

The update process of divide and conquer failed.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine sbdsqr (character UPLO, integer N, integer NCVT, integer NRU, integer NCC, real, dimension( * ) D, real, dimension( * ) E, real, dimension( ldvt, * ) VT, integer LDVT, real, dimension( ldu, * ) U, integer LDU, real, dimension( ldc, * ) C, integer LDC, real, dimension( * ) WORK, integer INFO)¶

**SBDSQR**

**Purpose:**

SBDSQR computes the singular values and, optionally, the right and/or

left singular vectors from the singular value decomposition (SVD) of

a real N-by-N (upper or lower) bidiagonal matrix B using the implicit

zero-shift QR algorithm. The SVD of B has the form

B = Q * S * P**T

where S is the diagonal matrix of singular values, Q is an orthogonal

matrix of left singular vectors, and P is an orthogonal matrix of

right singular vectors. If left singular vectors are requested, this

subroutine actually returns U*Q instead of Q, and, if right singular

vectors are requested, this subroutine returns P**T*VT instead of

P**T, for given real input matrices U and VT. When U and VT are the

orthogonal matrices that reduce a general matrix A to bidiagonal

form: A = U*B*VT, as computed by SGEBRD, then

A = (U*Q) * S * (P**T*VT)

is the SVD of A. Optionally, the subroutine may also compute Q**T*C

for a given real input matrix C.

See "Computing Small Singular Values of Bidiagonal Matrices With

Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,

LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,

no. 5, pp. 873-912, Sept 1990) and

"Accurate singular values and differential qd algorithms," by

B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics

Department, University of California at Berkeley, July 1992

for a detailed description of the algorithm.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

= 'U': B is upper bidiagonal;

= 'L': B is lower bidiagonal.

*N*

N is INTEGER

The order of the matrix B. N >= 0.

*NCVT*

NCVT is INTEGER

The number of columns of the matrix VT. NCVT >= 0.

*NRU*

NRU is INTEGER

The number of rows of the matrix U. NRU >= 0.

*NCC*

NCC is INTEGER

The number of columns of the matrix C. NCC >= 0.

*D*

D is REAL array, dimension (N)

On entry, the n diagonal elements of the bidiagonal matrix B.

On exit, if INFO=0, the singular values of B in decreasing

order.

*E*

E is REAL array, dimension (N-1)

On entry, the N-1 offdiagonal elements of the bidiagonal

matrix B.

On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E

will contain the diagonal and superdiagonal elements of a

bidiagonal matrix orthogonally equivalent to the one given

as input.

*VT*

VT is REAL array, dimension (LDVT, NCVT)

On entry, an N-by-NCVT matrix VT.

On exit, VT is overwritten by P**T * VT.

Not referenced if NCVT = 0.

*LDVT*

LDVT is INTEGER

The leading dimension of the array VT.

LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.

*U*

U is REAL array, dimension (LDU, N)

On entry, an NRU-by-N matrix U.

On exit, U is overwritten by U * Q.

Not referenced if NRU = 0.

*LDU*

LDU is INTEGER

The leading dimension of the array U. LDU >= max(1,NRU).

*C*

C is REAL array, dimension (LDC, NCC)

On entry, an N-by-NCC matrix C.

On exit, C is overwritten by Q**T * C.

Not referenced if NCC = 0.

*LDC*

LDC is INTEGER

The leading dimension of the array C.

LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.

*WORK*

WORK is REAL array, dimension (4*N)

*INFO*

INFO is INTEGER

= 0: successful exit

< 0: If INFO = -i, the i-th argument had an illegal value

> 0:

if NCVT = NRU = NCC = 0,

= 1, a split was marked by a positive value in E

= 2, current block of Z not diagonalized after 30*N

iterations (in inner while loop)

= 3, termination criterion of outer while loop not met

(program created more than N unreduced blocks)

else NCVT = NRU = NCC = 0,

the algorithm did not converge; D and E contain the

elements of a bidiagonal matrix which is orthogonally

similar to the input matrix B; if INFO = i, i

elements of E have not converged to zero.

**Internal Parameters:**

TOLMUL REAL, default = max(10,min(100,EPS**(-1/8)))

TOLMUL controls the convergence criterion of the QR loop.

If it is positive, TOLMUL*EPS is the desired relative

precision in the computed singular values.

If it is negative, abs(TOLMUL*EPS*sigma_max) is the

desired absolute accuracy in the computed singular

values (corresponds to relative accuracy

abs(TOLMUL*EPS) in the largest singular value.

abs(TOLMUL) should be between 1 and 1/EPS, and preferably

between 10 (for fast convergence) and .1/EPS

(for there to be some accuracy in the results).

Default is to lose at either one eighth or 2 of the

available decimal digits in each computed singular value

(whichever is smaller).

MAXITR INTEGER, default = 6

MAXITR controls the maximum number of passes of the

algorithm through its inner loop. The algorithms stops

(and so fails to converge) if the number of passes

through the inner loop exceeds MAXITR*N**2.

**Note:**

Bug report from Cezary Dendek.

On March 23rd 2017, the INTEGER variable MAXIT = MAXITR*N**2 is

removed since it can overflow pretty easily (for N larger or equal

than 18,919). We instead use MAXITDIVN = MAXITR*N.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine sdisna (character JOB, integer M, integer N, real, dimension( * ) D, real, dimension( * ) SEP, integer INFO)¶

**SDISNA**

**Purpose:**

SDISNA computes the reciprocal condition numbers for the eigenvectors

of a real symmetric or complex Hermitian matrix or for the left or

right singular vectors of a general m-by-n matrix. The reciprocal

condition number is the 'gap' between the corresponding eigenvalue or

singular value and the nearest other one.

The bound on the error, measured by angle in radians, in the I-th

computed vector is given by

SLAMCH( 'E' ) * ( ANORM / SEP( I ) )

where ANORM = 2-norm(A) = max( abs( D(j) ) ). SEP(I) is not allowed

to be smaller than SLAMCH( 'E' )*ANORM in order to limit the size of

the error bound.

SDISNA may also be used to compute error bounds for eigenvectors of

the generalized symmetric definite eigenproblem.

**Parameters**

*JOB*

JOB is CHARACTER*1

Specifies for which problem the reciprocal condition numbers

should be computed:

= 'E': the eigenvectors of a symmetric/Hermitian matrix;

= 'L': the left singular vectors of a general matrix;

= 'R': the right singular vectors of a general matrix.

*M*

M is INTEGER

The number of rows of the matrix. M >= 0.

*N*

N is INTEGER

If JOB = 'L' or 'R', the number of columns of the matrix,

in which case N >= 0. Ignored if JOB = 'E'.

*D*

D is REAL array, dimension (M) if JOB = 'E'

dimension (min(M,N)) if JOB = 'L' or 'R'

The eigenvalues (if JOB = 'E') or singular values (if JOB =

'L' or 'R') of the matrix, in either increasing or decreasing

order. If singular values, they must be non-negative.

*SEP*

SEP is REAL array, dimension (M) if JOB = 'E'

dimension (min(M,N)) if JOB = 'L' or 'R'

The reciprocal condition numbers of the vectors.

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine slaed0 (integer ICOMPQ, integer QSIZ, integer N, real, dimension( * ) D, real, dimension( * ) E, real, dimension( ldq, * ) Q, integer LDQ, real, dimension( ldqs, * ) QSTORE, integer LDQS, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)¶

**SLAED0** used by SSTEDC. Computes all eigenvalues and
corresponding eigenvectors of an unreduced symmetric tridiagonal matrix
using the divide and conquer method.

**Purpose:**

SLAED0 computes all eigenvalues and corresponding eigenvectors of a

symmetric tridiagonal matrix using the divide and conquer method.

**Parameters**

*ICOMPQ*

ICOMPQ is INTEGER

= 0: Compute eigenvalues only.

= 1: Compute eigenvectors of original dense symmetric matrix

also. On entry, Q contains the orthogonal matrix used

to reduce the original matrix to tridiagonal form.

= 2: Compute eigenvalues and eigenvectors of tridiagonal

matrix.

*QSIZ*

QSIZ is INTEGER

The dimension of the orthogonal matrix used to reduce

the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.

*N*

N is INTEGER

The dimension of the symmetric tridiagonal matrix. N >= 0.

*D*

D is REAL array, dimension (N)

On entry, the main diagonal of the tridiagonal matrix.

On exit, its eigenvalues.

*E*

E is REAL array, dimension (N-1)

The off-diagonal elements of the tridiagonal matrix.

On exit, E has been destroyed.

*Q*

Q is REAL array, dimension (LDQ, N)

On entry, Q must contain an N-by-N orthogonal matrix.

If ICOMPQ = 0 Q is not referenced.

If ICOMPQ = 1 On entry, Q is a subset of the columns of the

orthogonal matrix used to reduce the full

matrix to tridiagonal form corresponding to

the subset of the full matrix which is being

decomposed at this time.

If ICOMPQ = 2 On entry, Q will be the identity matrix.

On exit, Q contains the eigenvectors of the

tridiagonal matrix.

*LDQ*

LDQ is INTEGER

The leading dimension of the array Q. If eigenvectors are

desired, then LDQ >= max(1,N). In any case, LDQ >= 1.

*QSTORE*

QSTORE is REAL array, dimension (LDQS, N)

Referenced only when ICOMPQ = 1. Used to store parts of

the eigenvector matrix when the updating matrix multiplies

take place.

*LDQS*

LDQS is INTEGER

The leading dimension of the array QSTORE. If ICOMPQ = 1,

then LDQS >= max(1,N). In any case, LDQS >= 1.

*WORK*

WORK is REAL array,

If ICOMPQ = 0 or 1, the dimension of WORK must be at least

1 + 3*N + 2*N*lg N + 3*N**2

( lg( N ) = smallest integer k

such that 2^k >= N )

If ICOMPQ = 2, the dimension of WORK must be at least

4*N + N**2.

*IWORK*

IWORK is INTEGER array,

If ICOMPQ = 0 or 1, the dimension of IWORK must be at least

6 + 6*N + 5*N*lg N.

( lg( N ) = smallest integer k

such that 2^k >= N )

If ICOMPQ = 2, the dimension of IWORK must be at least

3 + 5*N.

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

> 0: The algorithm failed to compute an eigenvalue while

working on the submatrix lying in rows and columns

INFO/(N+1) through mod(INFO,N+1).

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine slaed1 (integer N, real, dimension( * ) D, real, dimension( ldq, * ) Q, integer LDQ, integer, dimension( * ) INDXQ, real RHO, integer CUTPNT, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)¶

**SLAED1** used by SSTEDC. Computes the updated eigensystem of
a diagonal matrix after modification by a rank-one symmetric matrix. Used
when the original matrix is tridiagonal.

**Purpose:**

SLAED1 computes the updated eigensystem of a diagonal

matrix after modification by a rank-one symmetric matrix. This

routine is used only for the eigenproblem which requires all

eigenvalues and eigenvectors of a tridiagonal matrix. SLAED7 handles

the case in which eigenvalues only or eigenvalues and eigenvectors

of a full symmetric matrix (which was reduced to tridiagonal form)

are desired.

T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)

where Z = Q**T*u, u is a vector of length N with ones in the

CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.

The eigenvectors of the original matrix are stored in Q, and the

eigenvalues are in D. The algorithm consists of three stages:

The first stage consists of deflating the size of the problem

when there are multiple eigenvalues or if there is a zero in

the Z vector. For each such occurrence the dimension of the

secular equation problem is reduced by one. This stage is

performed by the routine SLAED2.

The second stage consists of calculating the updated

eigenvalues. This is done by finding the roots of the secular

equation via the routine SLAED4 (as called by SLAED3).

This routine also calculates the eigenvectors of the current

problem.

The final stage consists of computing the updated eigenvectors

directly using the updated eigenvalues. The eigenvectors for

the current problem are multiplied with the eigenvectors from

the overall problem.

**Parameters**

*N*

N is INTEGER

The dimension of the symmetric tridiagonal matrix. N >= 0.

*D*

D is REAL array, dimension (N)

On entry, the eigenvalues of the rank-1-perturbed matrix.

On exit, the eigenvalues of the repaired matrix.

*Q*

Q is REAL array, dimension (LDQ,N)

On entry, the eigenvectors of the rank-1-perturbed matrix.

On exit, the eigenvectors of the repaired tridiagonal matrix.

*LDQ*

LDQ is INTEGER

The leading dimension of the array Q. LDQ >= max(1,N).

*INDXQ*

INDXQ is INTEGER array, dimension (N)

On entry, the permutation which separately sorts the two

subproblems in D into ascending order.

On exit, the permutation which will reintegrate the

subproblems back into sorted order,

i.e. D( INDXQ( I = 1, N ) ) will be in ascending order.

*RHO*

RHO is REAL

The subdiagonal entry used to create the rank-1 modification.

*CUTPNT*

CUTPNT is INTEGER

The location of the last eigenvalue in the leading sub-matrix.

min(1,N) <= CUTPNT <= N/2.

*WORK*

WORK is REAL array, dimension (4*N + N**2)

*IWORK*

IWORK is INTEGER array, dimension (4*N)

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

> 0: if INFO = 1, an eigenvalue did not converge

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

Modified by Francoise Tisseur, University of Tennessee

## subroutine slaed2 (integer K, integer N, integer N1, real, dimension( * ) D, real, dimension( ldq, * ) Q, integer LDQ, integer, dimension( * ) INDXQ, real RHO, real, dimension( * ) Z, real, dimension( * ) DLAMDA, real, dimension( * ) W, real, dimension( * ) Q2, integer, dimension( * ) INDX, integer, dimension( * ) INDXC, integer, dimension( * ) INDXP, integer, dimension( * ) COLTYP, integer INFO)¶

**SLAED2** used by SSTEDC. Merges eigenvalues and deflates
secular equation. Used when the original matrix is tridiagonal.

**Purpose:**

SLAED2 merges the two sets of eigenvalues together into a single

sorted set. Then it tries to deflate the size of the problem.

There are two ways in which deflation can occur: when two or more

eigenvalues are close together or if there is a tiny entry in the

Z vector. For each such occurrence the order of the related secular

equation problem is reduced by one.

**Parameters**

*K*

K is INTEGER

The number of non-deflated eigenvalues, and the order of the

related secular equation. 0 <= K <=N.

*N*

N is INTEGER

The dimension of the symmetric tridiagonal matrix. N >= 0.

*N1*

N1 is INTEGER

The location of the last eigenvalue in the leading sub-matrix.

min(1,N) <= N1 <= N/2.

*D*

D is REAL array, dimension (N)

On entry, D contains the eigenvalues of the two submatrices to

be combined.

On exit, D contains the trailing (N-K) updated eigenvalues

(those which were deflated) sorted into increasing order.

*Q*

Q is REAL array, dimension (LDQ, N)

On entry, Q contains the eigenvectors of two submatrices in

the two square blocks with corners at (1,1), (N1,N1)

and (N1+1, N1+1), (N,N).

On exit, Q contains the trailing (N-K) updated eigenvectors

(those which were deflated) in its last N-K columns.

*LDQ*

LDQ is INTEGER

The leading dimension of the array Q. LDQ >= max(1,N).

*INDXQ*

INDXQ is INTEGER array, dimension (N)

The permutation which separately sorts the two sub-problems

in D into ascending order. Note that elements in the second

half of this permutation must first have N1 added to their

values. Destroyed on exit.

*RHO*

RHO is REAL

On entry, the off-diagonal element associated with the rank-1

cut which originally split the two submatrices which are now

being recombined.

On exit, RHO has been modified to the value required by

SLAED3.

*Z*

Z is REAL array, dimension (N)

On entry, Z contains the updating vector (the last

row of the first sub-eigenvector matrix and the first row of

the second sub-eigenvector matrix).

On exit, the contents of Z have been destroyed by the updating

process.

*DLAMDA*

DLAMDA is REAL array, dimension (N)

A copy of the first K eigenvalues which will be used by

SLAED3 to form the secular equation.

*W*

W is REAL array, dimension (N)

The first k values of the final deflation-altered z-vector

which will be passed to SLAED3.

*Q2*

Q2 is REAL array, dimension (N1**2+(N-N1)**2)

A copy of the first K eigenvectors which will be used by

SLAED3 in a matrix multiply (SGEMM) to solve for the new

eigenvectors.

*INDX*

INDX is INTEGER array, dimension (N)

The permutation used to sort the contents of DLAMDA into

ascending order.

*INDXC*

INDXC is INTEGER array, dimension (N)

The permutation used to arrange the columns of the deflated

Q matrix into three groups: the first group contains non-zero

elements only at and above N1, the second contains

non-zero elements only below N1, and the third is dense.

*INDXP*

INDXP is INTEGER array, dimension (N)

The permutation used to place deflated values of D at the end

of the array. INDXP(1:K) points to the nondeflated D-values

and INDXP(K+1:N) points to the deflated eigenvalues.

*COLTYP*

COLTYP is INTEGER array, dimension (N)

During execution, a label which will indicate which of the

following types a column in the Q2 matrix is:

1 : non-zero in the upper half only;

2 : dense;

3 : non-zero in the lower half only;

4 : deflated.

On exit, COLTYP(i) is the number of columns of type i,

for i=1 to 4 only.

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

Modified by Francoise Tisseur, University of Tennessee

## subroutine slaed3 (integer K, integer N, integer N1, real, dimension( * ) D, real, dimension( ldq, * ) Q, integer LDQ, real RHO, real, dimension( * ) DLAMDA, real, dimension( * ) Q2, integer, dimension( * ) INDX, integer, dimension( * ) CTOT, real, dimension( * ) W, real, dimension( * ) S, integer INFO)¶

**SLAED3** used by SSTEDC. Finds the roots of the secular
equation and updates the eigenvectors. Used when the original matrix is
tridiagonal.

**Purpose:**

SLAED3 finds the roots of the secular equation, as defined by the

values in D, W, and RHO, between 1 and K. It makes the

appropriate calls to SLAED4 and then updates the eigenvectors by

multiplying the matrix of eigenvectors of the pair of eigensystems

being combined by the matrix of eigenvectors of the K-by-K system

which is solved here.

This code makes very mild assumptions about floating point

arithmetic. It will work on machines with a guard digit in

add/subtract, or on those binary machines without guard digits

which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.

It could conceivably fail on hexadecimal or decimal machines

without guard digits, but we know of none.

**Parameters**

*K*

K is INTEGER

The number of terms in the rational function to be solved by

SLAED4. K >= 0.

*N*

N is INTEGER

The number of rows and columns in the Q matrix.

N >= K (deflation may result in N>K).

*N1*

N1 is INTEGER

The location of the last eigenvalue in the leading submatrix.

min(1,N) <= N1 <= N/2.

*D*

D is REAL array, dimension (N)

D(I) contains the updated eigenvalues for

1 <= I <= K.

*Q*

Q is REAL array, dimension (LDQ,N)

Initially the first K columns are used as workspace.

On output the columns 1 to K contain

the updated eigenvectors.

*LDQ*

LDQ is INTEGER

The leading dimension of the array Q. LDQ >= max(1,N).

*RHO*

RHO is REAL

The value of the parameter in the rank one update equation.

RHO >= 0 required.

*DLAMDA*

DLAMDA is REAL array, dimension (K)

The first K elements of this array contain the old roots

of the deflated updating problem. These are the poles

of the secular equation. May be changed on output by

having lowest order bit set to zero on Cray X-MP, Cray Y-MP,

Cray-2, or Cray C-90, as described above.

*Q2*

Q2 is REAL array, dimension (LDQ2*N)

The first K columns of this matrix contain the non-deflated

eigenvectors for the split problem.

*INDX*

INDX is INTEGER array, dimension (N)

The permutation used to arrange the columns of the deflated

Q matrix into three groups (see SLAED2).

The rows of the eigenvectors found by SLAED4 must be likewise

permuted before the matrix multiply can take place.

*CTOT*

CTOT is INTEGER array, dimension (4)

A count of the total number of the various types of columns

in Q, as described in INDX. The fourth column type is any

column which has been deflated.

*W*

W is REAL array, dimension (K)

The first K elements of this array contain the components

of the deflation-adjusted updating vector. Destroyed on

output.

*S*

S is REAL array, dimension (N1 + 1)*K

Will contain the eigenvectors of the repaired matrix which

will be multiplied by the previously accumulated eigenvectors

to update the system.

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

> 0: if INFO = 1, an eigenvalue did not converge

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

Modified by Francoise Tisseur, University of Tennessee

## subroutine slaed4 (integer N, integer I, real, dimension( * ) D, real, dimension( * ) Z, real, dimension( * ) DELTA, real RHO, real DLAM, integer INFO)¶

**SLAED4** used by SSTEDC. Finds a single root of the secular
equation.

**Purpose:**

This subroutine computes the I-th updated eigenvalue of a symmetric

rank-one modification to a diagonal matrix whose elements are

given in the array d, and that

D(i) < D(j) for i < j

and that RHO > 0. This is arranged by the calling routine, and is

no loss in generality. The rank-one modified system is thus

diag( D ) + RHO * Z * Z_transpose.

where we assume the Euclidean norm of Z is 1.

The method consists of approximating the rational functions in the

secular equation by simpler interpolating rational functions.

**Parameters**

*N*

N is INTEGER

The length of all arrays.

*I*

I is INTEGER

The index of the eigenvalue to be computed. 1 <= I <= N.

*D*

D is REAL array, dimension (N)

The original eigenvalues. It is assumed that they are in

order, D(I) < D(J) for I < J.

*Z*

Z is REAL array, dimension (N)

The components of the updating vector.

*DELTA*

DELTA is REAL array, dimension (N)

If N > 2, DELTA contains (D(j) - lambda_I) in its j-th

component. If N = 1, then DELTA(1) = 1. If N = 2, see SLAED5

for detail. The vector DELTA contains the information necessary

to construct the eigenvectors by SLAED3 and SLAED9.

*RHO*

RHO is REAL

The scalar in the symmetric updating formula.

*DLAM*

DLAM is REAL

The computed lambda_I, the I-th updated eigenvalue.

*INFO*

INFO is INTEGER

= 0: successful exit

> 0: if INFO = 1, the updating process failed.

**Internal Parameters:**

Logical variable ORGATI (origin-at-i?) is used for distinguishing

whether D(i) or D(i+1) is treated as the origin.

ORGATI = .true. origin at i

ORGATI = .false. origin at i+1

Logical variable SWTCH3 (switch-for-3-poles?) is for noting

if we are working with THREE poles!

MAXIT is the maximum number of iterations allowed for each

eigenvalue.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine slaed5 (integer I, real, dimension( 2 ) D, real, dimension( 2 ) Z, real, dimension( 2 ) DELTA, real RHO, real DLAM)¶

**SLAED5** used by SSTEDC. Solves the 2-by-2 secular
equation.

**Purpose:**

This subroutine computes the I-th eigenvalue of a symmetric rank-one

modification of a 2-by-2 diagonal matrix

diag( D ) + RHO * Z * transpose(Z) .

The diagonal elements in the array D are assumed to satisfy

D(i) < D(j) for i < j .

We also assume RHO > 0 and that the Euclidean norm of the vector

Z is one.

**Parameters**

*I*

I is INTEGER

The index of the eigenvalue to be computed. I = 1 or I = 2.

*D*

D is REAL array, dimension (2)

The original eigenvalues. We assume D(1) < D(2).

*Z*

Z is REAL array, dimension (2)

The components of the updating vector.

*DELTA*

DELTA is REAL array, dimension (2)

The vector DELTA contains the information necessary

to construct the eigenvectors.

*RHO*

RHO is REAL

The scalar in the symmetric updating formula.

*DLAM*

DLAM is REAL

The computed lambda_I, the I-th updated eigenvalue.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine slaed6 (integer KNITER, logical ORGATI, real RHO, real, dimension( 3 ) D, real, dimension( 3 ) Z, real FINIT, real TAU, integer INFO)¶

**SLAED6** used by SSTEDC. Computes one Newton step in solution
of the secular equation.

**Purpose:**

SLAED6 computes the positive or negative root (closest to the origin)

of

z(1) z(2) z(3)

f(x) = rho + --------- + ---------- + ---------

d(1)-x d(2)-x d(3)-x

It is assumed that

if ORGATI = .true. the root is between d(2) and d(3);

otherwise it is between d(1) and d(2)

This routine will be called by SLAED4 when necessary. In most cases,

the root sought is the smallest in magnitude, though it might not be

in some extremely rare situations.

**Parameters**

*KNITER*

KNITER is INTEGER

Refer to SLAED4 for its significance.

*ORGATI*

ORGATI is LOGICAL

If ORGATI is true, the needed root is between d(2) and

d(3); otherwise it is between d(1) and d(2). See

SLAED4 for further details.

*RHO*

RHO is REAL

Refer to the equation f(x) above.

*D*

D is REAL array, dimension (3)

D satisfies d(1) < d(2) < d(3).

*Z*

Z is REAL array, dimension (3)

Each of the elements in z must be positive.

*FINIT*

FINIT is REAL

The value of f at 0. It is more accurate than the one

evaluated inside this routine (if someone wants to do

so).

*TAU*

TAU is REAL

The root of the equation f(x).

*INFO*

INFO is INTEGER

= 0: successful exit

> 0: if INFO = 1, failure to converge

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

10/02/03: This version has a few statements commented out for thread

safety (machine parameters are computed on each entry). SJH.

05/10/06: Modified from a new version of Ren-Cang Li, use

Gragg-Thornton-Warner cubic convergent scheme for better stability.

**Contributors:**

## subroutine slaed7 (integer ICOMPQ, integer N, integer QSIZ, integer TLVLS, integer CURLVL, integer CURPBM, real, dimension( * ) D, real, dimension( ldq, * ) Q, integer LDQ, integer, dimension( * ) INDXQ, real RHO, integer CUTPNT, real, dimension( * ) QSTORE, integer, dimension( * ) QPTR, integer, dimension( * ) PRMPTR, integer, dimension( * ) PERM, integer, dimension( * ) GIVPTR, integer, dimension( 2, * ) GIVCOL, real, dimension( 2, * ) GIVNUM, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)¶

**SLAED7** used by SSTEDC. Computes the updated eigensystem of
a diagonal matrix after modification by a rank-one symmetric matrix. Used
when the original matrix is dense.

**Purpose:**

SLAED7 computes the updated eigensystem of a diagonal

matrix after modification by a rank-one symmetric matrix. This

routine is used only for the eigenproblem which requires all

eigenvalues and optionally eigenvectors of a dense symmetric matrix

that has been reduced to tridiagonal form. SLAED1 handles

the case in which all eigenvalues and eigenvectors of a symmetric

tridiagonal matrix are desired.

T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)

where Z = Q**Tu, u is a vector of length N with ones in the

CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.

The eigenvectors of the original matrix are stored in Q, and the

eigenvalues are in D. The algorithm consists of three stages:

The first stage consists of deflating the size of the problem

when there are multiple eigenvalues or if there is a zero in

the Z vector. For each such occurrence the dimension of the

secular equation problem is reduced by one. This stage is

performed by the routine SLAED8.

The second stage consists of calculating the updated

eigenvalues. This is done by finding the roots of the secular

equation via the routine SLAED4 (as called by SLAED9).

This routine also calculates the eigenvectors of the current

problem.

The final stage consists of computing the updated eigenvectors

directly using the updated eigenvalues. The eigenvectors for

the current problem are multiplied with the eigenvectors from

the overall problem.

**Parameters**

*ICOMPQ*

ICOMPQ is INTEGER

= 0: Compute eigenvalues only.

= 1: Compute eigenvectors of original dense symmetric matrix

also. On entry, Q contains the orthogonal matrix used

to reduce the original matrix to tridiagonal form.

*N*

N is INTEGER

The dimension of the symmetric tridiagonal matrix. N >= 0.

*QSIZ*

QSIZ is INTEGER

The dimension of the orthogonal matrix used to reduce

the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.

*TLVLS*

TLVLS is INTEGER

The total number of merging levels in the overall divide and

conquer tree.

*CURLVL*

CURLVL is INTEGER

The current level in the overall merge routine,

0 <= CURLVL <= TLVLS.

*CURPBM*

CURPBM is INTEGER

The current problem in the current level in the overall

merge routine (counting from upper left to lower right).

*D*

D is REAL array, dimension (N)

On entry, the eigenvalues of the rank-1-perturbed matrix.

On exit, the eigenvalues of the repaired matrix.

*Q*

Q is REAL array, dimension (LDQ, N)

On entry, the eigenvectors of the rank-1-perturbed matrix.

On exit, the eigenvectors of the repaired tridiagonal matrix.

*LDQ*

LDQ is INTEGER

The leading dimension of the array Q. LDQ >= max(1,N).

*INDXQ*

INDXQ is INTEGER array, dimension (N)

The permutation which will reintegrate the subproblem just

solved back into sorted order, i.e., D( INDXQ( I = 1, N ) )

will be in ascending order.

*RHO*

RHO is REAL

The subdiagonal element used to create the rank-1

modification.

*CUTPNT*

CUTPNT is INTEGER

Contains the location of the last eigenvalue in the leading

sub-matrix. min(1,N) <= CUTPNT <= N.

*QSTORE*

QSTORE is REAL array, dimension (N**2+1)

Stores eigenvectors of submatrices encountered during

divide and conquer, packed together. QPTR points to

beginning of the submatrices.

*QPTR*

QPTR is INTEGER array, dimension (N+2)

List of indices pointing to beginning of submatrices stored

in QSTORE. The submatrices are numbered starting at the

bottom left of the divide and conquer tree, from left to

right and bottom to top.

*PRMPTR*

PRMPTR is INTEGER array, dimension (N lg N)

Contains a list of pointers which indicate where in PERM a

level's permutation is stored. PRMPTR(i+1) - PRMPTR(i)

indicates the size of the permutation and also the size of

the full, non-deflated problem.

*PERM*

PERM is INTEGER array, dimension (N lg N)

Contains the permutations (from deflation and sorting) to be

applied to each eigenblock.

*GIVPTR*

GIVPTR is INTEGER array, dimension (N lg N)

Contains a list of pointers which indicate where in GIVCOL a

level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i)

indicates the number of Givens rotations.

*GIVCOL*

GIVCOL is INTEGER array, dimension (2, N lg N)

Each pair of numbers indicates a pair of columns to take place

in a Givens rotation.

*GIVNUM*

GIVNUM is REAL array, dimension (2, N lg N)

Each number indicates the S value to be used in the

corresponding Givens rotation.

*WORK*

WORK is REAL array, dimension (3*N+2*QSIZ*N)

*IWORK*

IWORK is INTEGER array, dimension (4*N)

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

> 0: if INFO = 1, an eigenvalue did not converge

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine slaed8 (integer ICOMPQ, integer K, integer N, integer QSIZ, real, dimension( * ) D, real, dimension( ldq, * ) Q, integer LDQ, integer, dimension( * ) INDXQ, real RHO, integer CUTPNT, real, dimension( * ) Z, real, dimension( * ) DLAMDA, real, dimension( ldq2, * ) Q2, integer LDQ2, real, dimension( * ) W, integer, dimension( * ) PERM, integer GIVPTR, integer, dimension( 2, * ) GIVCOL, real, dimension( 2, * ) GIVNUM, integer, dimension( * ) INDXP, integer, dimension( * ) INDX, integer INFO)¶

**SLAED8** used by SSTEDC. Merges eigenvalues and deflates
secular equation. Used when the original matrix is dense.

**Purpose:**

SLAED8 merges the two sets of eigenvalues together into a single

sorted set. Then it tries to deflate the size of the problem.

There are two ways in which deflation can occur: when two or more

eigenvalues are close together or if there is a tiny element in the

Z vector. For each such occurrence the order of the related secular

equation problem is reduced by one.

**Parameters**

*ICOMPQ*

ICOMPQ is INTEGER

= 0: Compute eigenvalues only.

= 1: Compute eigenvectors of original dense symmetric matrix

also. On entry, Q contains the orthogonal matrix used

to reduce the original matrix to tridiagonal form.

*K*

K is INTEGER

The number of non-deflated eigenvalues, and the order of the

related secular equation.

*N*

N is INTEGER

The dimension of the symmetric tridiagonal matrix. N >= 0.

*QSIZ*

QSIZ is INTEGER

The dimension of the orthogonal matrix used to reduce

the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.

*D*

D is REAL array, dimension (N)

On entry, the eigenvalues of the two submatrices to be

combined. On exit, the trailing (N-K) updated eigenvalues

(those which were deflated) sorted into increasing order.

*Q*

Q is REAL array, dimension (LDQ,N)

If ICOMPQ = 0, Q is not referenced. Otherwise,

on entry, Q contains the eigenvectors of the partially solved

system which has been previously updated in matrix

multiplies with other partially solved eigensystems.

On exit, Q contains the trailing (N-K) updated eigenvectors

(those which were deflated) in its last N-K columns.

*LDQ*

LDQ is INTEGER

The leading dimension of the array Q. LDQ >= max(1,N).

*INDXQ*

INDXQ is INTEGER array, dimension (N)

The permutation which separately sorts the two sub-problems

in D into ascending order. Note that elements in the second

half of this permutation must first have CUTPNT added to

their values in order to be accurate.

*RHO*

RHO is REAL

On entry, the off-diagonal element associated with the rank-1

cut which originally split the two submatrices which are now

being recombined.

On exit, RHO has been modified to the value required by

SLAED3.

*CUTPNT*

CUTPNT is INTEGER

The location of the last eigenvalue in the leading

sub-matrix. min(1,N) <= CUTPNT <= N.

*Z*

Z is REAL array, dimension (N)

On entry, Z contains the updating vector (the last row of

the first sub-eigenvector matrix and the first row of the

second sub-eigenvector matrix).

On exit, the contents of Z are destroyed by the updating

process.

*DLAMDA*

DLAMDA is REAL array, dimension (N)

A copy of the first K eigenvalues which will be used by

SLAED3 to form the secular equation.

*Q2*

Q2 is REAL array, dimension (LDQ2,N)

If ICOMPQ = 0, Q2 is not referenced. Otherwise,

a copy of the first K eigenvectors which will be used by

SLAED7 in a matrix multiply (SGEMM) to update the new

eigenvectors.

*LDQ2*

LDQ2 is INTEGER

The leading dimension of the array Q2. LDQ2 >= max(1,N).

*W*

W is REAL array, dimension (N)

The first k values of the final deflation-altered z-vector and

will be passed to SLAED3.

*PERM*

PERM is INTEGER array, dimension (N)

The permutations (from deflation and sorting) to be applied

to each eigenblock.

*GIVPTR*

GIVPTR is INTEGER

The number of Givens rotations which took place in this

subproblem.

*GIVCOL*

GIVCOL is INTEGER array, dimension (2, N)

Each pair of numbers indicates a pair of columns to take place

in a Givens rotation.

*GIVNUM*

GIVNUM is REAL array, dimension (2, N)

Each number indicates the S value to be used in the

corresponding Givens rotation.

*INDXP*

INDXP is INTEGER array, dimension (N)

The permutation used to place deflated values of D at the end

of the array. INDXP(1:K) points to the nondeflated D-values

and INDXP(K+1:N) points to the deflated eigenvalues.

*INDX*

INDX is INTEGER array, dimension (N)

The permutation used to sort the contents of D into ascending

order.

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine slaed9 (integer K, integer KSTART, integer KSTOP, integer N, real, dimension( * ) D, real, dimension( ldq, * ) Q, integer LDQ, real RHO, real, dimension( * ) DLAMDA, real, dimension( * ) W, real, dimension( lds, * ) S, integer LDS, integer INFO)¶

**SLAED9** used by SSTEDC. Finds the roots of the secular
equation and updates the eigenvectors. Used when the original matrix is
dense.

**Purpose:**

SLAED9 finds the roots of the secular equation, as defined by the

values in D, Z, and RHO, between KSTART and KSTOP. It makes the

appropriate calls to SLAED4 and then stores the new matrix of

eigenvectors for use in calculating the next level of Z vectors.

**Parameters**

*K*

K is INTEGER

The number of terms in the rational function to be solved by

SLAED4. K >= 0.

*KSTART*

KSTART is INTEGER

*KSTOP*

KSTOP is INTEGER

The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP

are to be computed. 1 <= KSTART <= KSTOP <= K.

*N*

N is INTEGER

The number of rows and columns in the Q matrix.

N >= K (delation may result in N > K).

*D*

D is REAL array, dimension (N)

D(I) contains the updated eigenvalues

for KSTART <= I <= KSTOP.

*Q*

Q is REAL array, dimension (LDQ,N)

*LDQ*

LDQ is INTEGER

The leading dimension of the array Q. LDQ >= max( 1, N ).

*RHO*

RHO is REAL

The value of the parameter in the rank one update equation.

RHO >= 0 required.

*DLAMDA*

DLAMDA is REAL array, dimension (K)

The first K elements of this array contain the old roots

of the deflated updating problem. These are the poles

of the secular equation.

*W*

W is REAL array, dimension (K)

The first K elements of this array contain the components

of the deflation-adjusted updating vector.

*S*

S is REAL array, dimension (LDS, K)

Will contain the eigenvectors of the repaired matrix which

will be stored for subsequent Z vector calculation and

multiplied by the previously accumulated eigenvectors

to update the system.

*LDS*

LDS is INTEGER

The leading dimension of S. LDS >= max( 1, K ).

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

> 0: if INFO = 1, an eigenvalue did not converge

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine slaeda (integer N, integer TLVLS, integer CURLVL, integer CURPBM, integer, dimension( * ) PRMPTR, integer, dimension( * ) PERM, integer, dimension( * ) GIVPTR, integer, dimension( 2, * ) GIVCOL, real, dimension( 2, * ) GIVNUM, real, dimension( * ) Q, integer, dimension( * ) QPTR, real, dimension( * ) Z, real, dimension( * ) ZTEMP, integer INFO)¶

**SLAEDA** used by SSTEDC. Computes the Z vector determining
the rank-one modification of the diagonal matrix. Used when the original
matrix is dense.

**Purpose:**

SLAEDA computes the Z vector corresponding to the merge step in the

CURLVLth step of the merge process with TLVLS steps for the CURPBMth

problem.

**Parameters**

*N*

N is INTEGER

The dimension of the symmetric tridiagonal matrix. N >= 0.

*TLVLS*

TLVLS is INTEGER

The total number of merging levels in the overall divide and

conquer tree.

*CURLVL*

CURLVL is INTEGER

The current level in the overall merge routine,

0 <= curlvl <= tlvls.

*CURPBM*

CURPBM is INTEGER

The current problem in the current level in the overall

merge routine (counting from upper left to lower right).

*PRMPTR*

PRMPTR is INTEGER array, dimension (N lg N)

Contains a list of pointers which indicate where in PERM a

level's permutation is stored. PRMPTR(i+1) - PRMPTR(i)

indicates the size of the permutation and incidentally the

size of the full, non-deflated problem.

*PERM*

PERM is INTEGER array, dimension (N lg N)

Contains the permutations (from deflation and sorting) to be

applied to each eigenblock.

*GIVPTR*

GIVPTR is INTEGER array, dimension (N lg N)

Contains a list of pointers which indicate where in GIVCOL a

level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i)

indicates the number of Givens rotations.

*GIVCOL*

GIVCOL is INTEGER array, dimension (2, N lg N)

Each pair of numbers indicates a pair of columns to take place

in a Givens rotation.

*GIVNUM*

GIVNUM is REAL array, dimension (2, N lg N)

Each number indicates the S value to be used in the

corresponding Givens rotation.

*Q*

Q is REAL array, dimension (N**2)

Contains the square eigenblocks from previous levels, the

starting positions for blocks are given by QPTR.

*QPTR*

QPTR is INTEGER array, dimension (N+2)

Contains a list of pointers which indicate where in Q an

eigenblock is stored. SQRT( QPTR(i+1) - QPTR(i) ) indicates

the size of the block.

*Z*

Z is REAL array, dimension (N)

On output this vector contains the updating vector (the last

row of the first sub-eigenvector matrix and the first row of

the second sub-eigenvector matrix).

*ZTEMP*

ZTEMP is REAL array, dimension (N)

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine slagtf (integer N, real, dimension( * ) A, real LAMBDA, real, dimension( * ) B, real, dimension( * ) C, real TOL, real, dimension( * ) D, integer, dimension( * ) IN, integer INFO)¶

**SLAGTF** computes an LU factorization of a matrix
T-λI, where T is a general tridiagonal matrix, and λ a scalar,
using partial pivoting with row interchanges.

**Purpose:**

SLAGTF factorizes the matrix (T - lambda*I), where T is an n by n

tridiagonal matrix and lambda is a scalar, as

T - lambda*I = PLU,

where P is a permutation matrix, L is a unit lower tridiagonal matrix

with at most one non-zero sub-diagonal elements per column and U is

an upper triangular matrix with at most two non-zero super-diagonal

elements per column.

The factorization is obtained by Gaussian elimination with partial

pivoting and implicit row scaling.

The parameter LAMBDA is included in the routine so that SLAGTF may

be used, in conjunction with SLAGTS, to obtain eigenvectors of T by

inverse iteration.

**Parameters**

*N*

N is INTEGER

The order of the matrix T.

*A*

A is REAL array, dimension (N)

On entry, A must contain the diagonal elements of T.

On exit, A is overwritten by the n diagonal elements of the

upper triangular matrix U of the factorization of T.

*LAMBDA*

LAMBDA is REAL

On entry, the scalar lambda.

*B*

B is REAL array, dimension (N-1)

On entry, B must contain the (n-1) super-diagonal elements of

T.

On exit, B is overwritten by the (n-1) super-diagonal

elements of the matrix U of the factorization of T.

*C*

C is REAL array, dimension (N-1)

On entry, C must contain the (n-1) sub-diagonal elements of

T.

On exit, C is overwritten by the (n-1) sub-diagonal elements

of the matrix L of the factorization of T.

*TOL*

TOL is REAL

On entry, a relative tolerance used to indicate whether or

not the matrix (T - lambda*I) is nearly singular. TOL should

normally be chose as approximately the largest relative error

in the elements of T. For example, if the elements of T are

correct to about 4 significant figures, then TOL should be

set to about 5*10**(-4). If TOL is supplied as less than eps,

where eps is the relative machine precision, then the value

eps is used in place of TOL.

*D*

D is REAL array, dimension (N-2)

On exit, D is overwritten by the (n-2) second super-diagonal

elements of the matrix U of the factorization of T.

*IN*

IN is INTEGER array, dimension (N)

On exit, IN contains details of the permutation matrix P. If

an interchange occurred at the kth step of the elimination,

then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)

returns the smallest positive integer j such that

abs( u(j,j) ) <= norm( (T - lambda*I)(j) )*TOL,

where norm( A(j) ) denotes the sum of the absolute values of

the jth row of the matrix A. If no such j exists then IN(n)

is returned as zero. If IN(n) is returned as positive, then a

diagonal element of U is small, indicating that

(T - lambda*I) is singular or nearly singular,

*INFO*

INFO is INTEGER

= 0: successful exit

< 0: if INFO = -k, the kth argument had an illegal value

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine slamrg (integer N1, integer N2, real, dimension( * ) A, integer STRD1, integer STRD2, integer, dimension( * ) INDEX)¶

**SLAMRG** creates a permutation list to merge the entries of
two independently sorted sets into a single set sorted in ascending
order.

**Purpose:**

SLAMRG will create a permutation list which will merge the elements

of A (which is composed of two independently sorted sets) into a

single set which is sorted in ascending order.

**Parameters**

*N1*

N1 is INTEGER

*N2*

N2 is INTEGER

These arguments contain the respective lengths of the two

sorted lists to be merged.

*A*

A is REAL array, dimension (N1+N2)

The first N1 elements of A contain a list of numbers which

are sorted in either ascending or descending order. Likewise

for the final N2 elements.

*STRD1*

STRD1 is INTEGER

*STRD2*

STRD2 is INTEGER

These are the strides to be taken through the array A.

Allowable strides are 1 and -1. They indicate whether a

subset of A is sorted in ascending (STRDx = 1) or descending

(STRDx = -1) order.

*INDEX*

INDEX is INTEGER array, dimension (N1+N2)

On exit this array will contain a permutation such that

if B( I ) = A( INDEX( I ) ) for I=1,N1+N2, then B will be

sorted in ascending order.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine slartgs (real X, real Y, real SIGMA, real CS, real SN)¶

**SLARTGS** generates a plane rotation designed to introduce a
bulge in implicit QR iteration for the bidiagonal SVD problem.

**Purpose:**

SLARTGS generates a plane rotation designed to introduce a bulge in

Golub-Reinsch-style implicit QR iteration for the bidiagonal SVD

problem. X and Y are the top-row entries, and SIGMA is the shift.

The computed CS and SN define a plane rotation satisfying

[ CS SN ] . [ X^2 - SIGMA ] = [ R ],

[ -SN CS ] [ X * Y ] [ 0 ]

with R nonnegative. If X^2 - SIGMA and X * Y are 0, then the

rotation is by PI/2.

**Parameters**

*X*

X is REAL

The (1,1) entry of an upper bidiagonal matrix.

*Y*

Y is REAL

The (1,2) entry of an upper bidiagonal matrix.

*SIGMA*

SIGMA is REAL

The shift.

*CS*

CS is REAL

The cosine of the rotation.

*SN*

SN is REAL

The sine of the rotation.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine slasq1 (integer N, real, dimension( * ) D, real, dimension( * ) E, real, dimension( * ) WORK, integer INFO)¶

**SLASQ1** computes the singular values of a real square
bidiagonal matrix. Used by sbdsqr.

**Purpose:**

SLASQ1 computes the singular values of a real N-by-N bidiagonal

matrix with diagonal D and off-diagonal E. The singular values

are computed to high relative accuracy, in the absence of

denormalization, underflow and overflow. The algorithm was first

presented in

"Accurate singular values and differential qd algorithms" by K. V.

Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,

1994,

and the present implementation is described in "An implementation of

the dqds Algorithm (Positive Case)", LAPACK Working Note.

**Parameters**

*N*

N is INTEGER

The number of rows and columns in the matrix. N >= 0.

*D*

D is REAL array, dimension (N)

On entry, D contains the diagonal elements of the

bidiagonal matrix whose SVD is desired. On normal exit,

D contains the singular values in decreasing order.

*E*

E is REAL array, dimension (N)

On entry, elements E(1:N-1) contain the off-diagonal elements

of the bidiagonal matrix whose SVD is desired.

On exit, E is overwritten.

*WORK*

WORK is REAL array, dimension (4*N)

*INFO*

INFO is INTEGER

= 0: successful exit

< 0: if INFO = -i, the i-th argument had an illegal value

> 0: the algorithm failed

= 1, a split was marked by a positive value in E

= 2, current block of Z not diagonalized after 100*N

iterations (in inner while loop) On exit D and E

represent a matrix with the same singular values

which the calling subroutine could use to finish the

computation, or even feed back into SLASQ1

= 3, termination criterion of outer while loop not met

(program created more than N unreduced blocks)

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine slasq2 (integer N, real, dimension( * ) Z, integer INFO)¶

**SLASQ2** computes all the eigenvalues of the symmetric
positive definite tridiagonal matrix associated with the qd Array Z to high
relative accuracy. Used by sbdsqr and sstegr.

**Purpose:**

SLASQ2 computes all the eigenvalues of the symmetric positive

definite tridiagonal matrix associated with the qd array Z to high

relative accuracy are computed to high relative accuracy, in the

absence of denormalization, underflow and overflow.

To see the relation of Z to the tridiagonal matrix, let L be a

unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and

let U be an upper bidiagonal matrix with 1's above and diagonal

Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the

symmetric tridiagonal to which it is similar.

Note : SLASQ2 defines a logical variable, IEEE, which is true

on machines which follow ieee-754 floating-point standard in their

handling of infinities and NaNs, and false otherwise. This variable

is passed to SLASQ3.

**Parameters**

*N*

N is INTEGER

The number of rows and columns in the matrix. N >= 0.

*Z*

Z is REAL array, dimension ( 4*N )

On entry Z holds the qd array. On exit, entries 1 to N hold

the eigenvalues in decreasing order, Z( 2*N+1 ) holds the

trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If

N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )

holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of

shifts that failed.

*INFO*

INFO is INTEGER

= 0: successful exit

< 0: if the i-th argument is a scalar and had an illegal

value, then INFO = -i, if the i-th argument is an

array and the j-entry had an illegal value, then

INFO = -(i*100+j)

> 0: the algorithm failed

= 1, a split was marked by a positive value in E

= 2, current block of Z not diagonalized after 100*N

iterations (in inner while loop). On exit Z holds

a qd array with the same eigenvalues as the given Z.

= 3, termination criterion of outer while loop not met

(program created more than N unreduced blocks)

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

Local Variables: I0:N0 defines a current unreduced segment of Z.

The shifts are accumulated in SIGMA. Iteration count is in ITER.

Ping-pong is controlled by PP (alternates between 0 and 1).

## subroutine slasq3 (integer I0, integer N0, real, dimension( * ) Z, integer PP, real DMIN, real SIGMA, real DESIG, real QMAX, integer NFAIL, integer ITER, integer NDIV, logical IEEE, integer TTYPE, real DMIN1, real DMIN2, real DN, real DN1, real DN2, real G, real TAU)¶

**SLASQ3** checks for deflation, computes a shift and calls
dqds. Used by sbdsqr.

**Purpose:**

SLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.

In case of failure it changes shifts, and tries again until output

is positive.

**Parameters**

*I0*

I0 is INTEGER

First index.

*N0*

N0 is INTEGER

Last index.

*Z*

Z is REAL array, dimension ( 4*N0 )

Z holds the qd array.

*PP*

PP is INTEGER

PP=0 for ping, PP=1 for pong.

PP=2 indicates that flipping was applied to the Z array

and that the initial tests for deflation should not be

performed.

*DMIN*

DMIN is REAL

Minimum value of d.

*SIGMA*

SIGMA is REAL

Sum of shifts used in current segment.

*DESIG*

DESIG is REAL

Lower order part of SIGMA

*QMAX*

QMAX is REAL

Maximum value of q.

*NFAIL*

NFAIL is INTEGER

Increment NFAIL by 1 each time the shift was too big.

*ITER*

ITER is INTEGER

Increment ITER by 1 for each iteration.

*NDIV*

NDIV is INTEGER

Increment NDIV by 1 for each division.

*IEEE*

IEEE is LOGICAL

Flag for IEEE or non IEEE arithmetic (passed to SLASQ5).

*TTYPE*

TTYPE is INTEGER

Shift type.

*DMIN1*

DMIN1 is REAL

*DMIN2*

DMIN2 is REAL

*DN*

DN is REAL

*DN1*

DN1 is REAL

*DN2*

DN2 is REAL

*G*

G is REAL

*TAU*

TAU is REAL

These are passed as arguments in order to save their values

between calls to SLASQ3.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine slasq4 (integer I0, integer N0, real, dimension( * ) Z, integer PP, integer N0IN, real DMIN, real DMIN1, real DMIN2, real DN, real DN1, real DN2, real TAU, integer TTYPE, real G)¶

**SLASQ4** computes an approximation to the smallest eigenvalue
using values of d from the previous transform. Used by sbdsqr.

**Purpose:**

SLASQ4 computes an approximation TAU to the smallest eigenvalue

using values of d from the previous transform.

**Parameters**

*I0*

I0 is INTEGER

First index.

*N0*

N0 is INTEGER

Last index.

*Z*

Z is REAL array, dimension ( 4*N0 )

Z holds the qd array.

*PP*

PP is INTEGER

PP=0 for ping, PP=1 for pong.

*N0IN*

N0IN is INTEGER

The value of N0 at start of EIGTEST.

*DMIN*

DMIN is REAL

Minimum value of d.

*DMIN1*

DMIN1 is REAL

Minimum value of d, excluding D( N0 ).

*DMIN2*

DMIN2 is REAL

Minimum value of d, excluding D( N0 ) and D( N0-1 ).

*DN*

DN is REAL

d(N)

*DN1*

DN1 is REAL

d(N-1)

*DN2*

DN2 is REAL

d(N-2)

*TAU*

TAU is REAL

This is the shift.

*TTYPE*

TTYPE is INTEGER

Shift type.

*G*

G is REAL

G is passed as an argument in order to save its value between

calls to SLASQ4.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

CNST1 = 9/16

## subroutine slasq5 (integer I0, integer N0, real, dimension( * ) Z, integer PP, real TAU, real SIGMA, real DMIN, real DMIN1, real DMIN2, real DN, real DNM1, real DNM2, logical IEEE, real EPS)¶

** SLASQ5 computes one dqds transform in ping-pong form. Used by
sbdsqr and sstegr. **

**Purpose:**

SLASQ5 computes one dqds transform in ping-pong form, one

version for IEEE machines another for non IEEE machines.

**Parameters**

*I0*

I0 is INTEGER

First index.

*N0*

N0 is INTEGER

Last index.

*Z*

Z is REAL array, dimension ( 4*N )

Z holds the qd array. EMIN is stored in Z(4*N0) to avoid

an extra argument.

*PP*

PP is INTEGER

PP=0 for ping, PP=1 for pong.

*TAU*

TAU is REAL

This is the shift.

*SIGMA*

SIGMA is REAL

This is the accumulated shift up to this step.

*DMIN*

DMIN is REAL

Minimum value of d.

*DMIN1*

DMIN1 is REAL

Minimum value of d, excluding D( N0 ).

*DMIN2*

DMIN2 is REAL

Minimum value of d, excluding D( N0 ) and D( N0-1 ).

*DN*

DN is REAL

d(N0), the last value of d.

*DNM1*

DNM1 is REAL

d(N0-1).

*DNM2*

DNM2 is REAL

d(N0-2).

*IEEE*

IEEE is LOGICAL

Flag for IEEE or non IEEE arithmetic.

*EPS*

EPS is REAL

This is the value of epsilon used.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine slasq6 (integer I0, integer N0, real, dimension( * ) Z, integer PP, real DMIN, real DMIN1, real DMIN2, real DN, real DNM1, real DNM2)¶

**SLASQ6** computes one dqd transform in ping-pong form. Used
by sbdsqr and sstegr.

**Purpose:**

SLASQ6 computes one dqd (shift equal to zero) transform in

ping-pong form, with protection against underflow and overflow.

**Parameters**

*I0*

I0 is INTEGER

First index.

*N0*

N0 is INTEGER

Last index.

*Z*

Z is REAL array, dimension ( 4*N )

Z holds the qd array. EMIN is stored in Z(4*N0) to avoid

an extra argument.

*PP*

PP is INTEGER

PP=0 for ping, PP=1 for pong.

*DMIN*

DMIN is REAL

Minimum value of d.

*DMIN1*

DMIN1 is REAL

Minimum value of d, excluding D( N0 ).

*DMIN2*

DMIN2 is REAL

Minimum value of d, excluding D( N0 ) and D( N0-1 ).

*DN*

DN is REAL

d(N0), the last value of d.

*DNM1*

DNM1 is REAL

d(N0-1).

*DNM2*

DNM2 is REAL

d(N0-2).

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine slasrt (character ID, integer N, real, dimension( * ) D, integer INFO)¶

**SLASRT** sorts numbers in increasing or decreasing order.

**Purpose:**

Sort the numbers in D in increasing order (if ID = 'I') or

in decreasing order (if ID = 'D' ).

Use Quick Sort, reverting to Insertion sort on arrays of

size <= 20. Dimension of STACK limits N to about 2**32.

**Parameters**

*ID*

ID is CHARACTER*1

= 'I': sort D in increasing order;

= 'D': sort D in decreasing order.

*N*

N is INTEGER

The length of the array D.

*D*

D is REAL array, dimension (N)

On entry, the array to be sorted.

On exit, D has been sorted into increasing order

(D(1) <= ... <= D(N) ) or into decreasing order

(D(1) >= ... >= D(N) ), depending on ID.

*INFO*

INFO is INTEGER

= 0: successful exit

< 0: if INFO = -i, the i-th argument had an illegal value

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine spttrf (integer N, real, dimension( * ) D, real, dimension( * ) E, integer INFO)¶

**SPTTRF**

**Purpose:**

SPTTRF computes the L*D*L**T factorization of a real symmetric

positive definite tridiagonal matrix A. The factorization may also

be regarded as having the form A = U**T*D*U.

**Parameters**

*N*

N is INTEGER

The order of the matrix A. N >= 0.

*D*

D is REAL array, dimension (N)

On entry, the n diagonal elements of the tridiagonal matrix

A. On exit, the n diagonal elements of the diagonal matrix

D from the L*D*L**T factorization of A.

*E*

E is REAL array, dimension (N-1)

On entry, the (n-1) subdiagonal elements of the tridiagonal

matrix A. On exit, the (n-1) subdiagonal elements of the

unit bidiagonal factor L from the L*D*L**T factorization of A.

E can also be regarded as the superdiagonal of the unit

bidiagonal factor U from the U**T*D*U factorization of A.

*INFO*

INFO is INTEGER

= 0: successful exit

< 0: if INFO = -k, the k-th argument had an illegal value

> 0: if INFO = k, the leading minor of order k is not

positive definite; if k < N, the factorization could not

be completed, while if k = N, the factorization was

completed, but D(N) <= 0.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine sstebz (character RANGE, character ORDER, integer N, real VL, real VU, integer IL, integer IU, real ABSTOL, real, dimension( * ) D, real, dimension( * ) E, integer M, integer NSPLIT, real, dimension( * ) W, integer, dimension( * ) IBLOCK, integer, dimension( * ) ISPLIT, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)¶

**SSTEBZ**

**Purpose:**

SSTEBZ computes the eigenvalues of a symmetric tridiagonal

matrix T. 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**

*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.

*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.

*N*

N is INTEGER

The order of the tridiagonal matrix T. N >= 0.

*VL*

VL is REAL

If RANGE='V', the lower bound 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'.

*VU*

VU is REAL

If RANGE='V', the upper bound 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'.

*IL*

IL is INTEGER

If RANGE='I', the index of the

smallest eigenvalue to be returned.

1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.

Not referenced if RANGE = 'A' or 'V'.

*IU*

IU is INTEGER

If RANGE='I', the index of the

largest eigenvalue to be returned.

1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.

Not referenced if RANGE = 'A' or 'V'.

*ABSTOL*

ABSTOL is REAL

The absolute tolerance for the eigenvalues. An eigenvalue

(or cluster) is considered to be located if it has been

determined to lie in an interval whose width is ABSTOL or

less. If ABSTOL is less than or equal to zero, then ULP*|T|

will be used, where |T| means the 1-norm of T.

Eigenvalues will be computed most accurately when ABSTOL is

set to twice the underflow threshold 2*SLAMCH('S'), not zero.

*D*

D is REAL array, dimension (N)

The n diagonal elements of the tridiagonal matrix T.

*E*

E is REAL array, dimension (N-1)

The (n-1) off-diagonal elements of the tridiagonal matrix T.

*M*

M is INTEGER

The actual number of eigenvalues found. 0 <= M <= N.

(See also the description of INFO=2,3.)

*NSPLIT*

NSPLIT is INTEGER

The number of diagonal blocks in the matrix T.

1 <= NSPLIT <= N.

*W*

W is REAL array, dimension (N)

On exit, the first M elements of W will contain the

eigenvalues. (SSTEBZ may use the remaining N-M elements as

workspace.)

*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. (SSTEBZ may use the remaining N-M elements as

workspace.)

*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.)

*WORK*

WORK is REAL array, dimension (4*N)

*IWORK*

IWORK is INTEGER array, dimension (3*N)

*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:**

RELFAC REAL, default = 2.0e0

The relative tolerance. An interval (a,b] lies within

"relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|),

where "ulp" is the machine precision (distance from 1 to

the next larger floating point number.)

FUDGE REAL, 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.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine sstedc (character COMPZ, integer N, real, dimension( * ) D, real, dimension( * ) E, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer LWORK, integer, dimension( * ) IWORK, integer LIWORK, integer INFO)¶

**SSTEDC**

**Purpose:**

SSTEDC computes all eigenvalues and, optionally, eigenvectors of a

symmetric tridiagonal matrix using the divide and conquer method.

The eigenvectors of a full or band real symmetric matrix can also be

found if SSYTRD or SSPTRD or SSBTRD has been used to reduce this

matrix to tridiagonal form.

This code makes very mild assumptions about floating point

arithmetic. It will work on machines with a guard digit in

add/subtract, or on those binary machines without guard digits

which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.

It could conceivably fail on hexadecimal or decimal machines

without guard digits, but we know of none. See SLAED3 for details.

**Parameters**

*COMPZ*

COMPZ is CHARACTER*1

= 'N': Compute eigenvalues only.

= 'I': Compute eigenvectors of tridiagonal matrix also.

= 'V': Compute eigenvectors of original dense symmetric

matrix also. On entry, Z contains the orthogonal

matrix used to reduce the original matrix to

tridiagonal form.

*N*

N is INTEGER

The dimension of the symmetric tridiagonal matrix. N >= 0.

*D*

D is REAL array, dimension (N)

On entry, the diagonal elements of the tridiagonal matrix.

On exit, if INFO = 0, the eigenvalues in ascending order.

*E*

E is REAL array, dimension (N-1)

On entry, the subdiagonal elements of the tridiagonal matrix.

On exit, E has been destroyed.

*Z*

Z is REAL array, dimension (LDZ,N)

On entry, if COMPZ = 'V', then Z contains the orthogonal

matrix used in the reduction to tridiagonal form.

On exit, if INFO = 0, then if COMPZ = 'V', Z contains the

orthonormal eigenvectors of the original symmetric matrix,

and if COMPZ = 'I', Z contains the orthonormal eigenvectors

of the symmetric tridiagonal matrix.

If COMPZ = 'N', then Z is not referenced.

*LDZ*

LDZ is INTEGER

The leading dimension of the array Z. LDZ >= 1.

If eigenvectors are desired, then LDZ >= max(1,N).

*WORK*

WORK is REAL array, dimension (MAX(1,LWORK))

On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

*LWORK*

LWORK is INTEGER

The dimension of the array WORK.

If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.

If COMPZ = 'V' and N > 1 then LWORK must be at least

( 1 + 3*N + 2*N*lg N + 4*N**2 ),

where lg( N ) = smallest integer k such

that 2**k >= N.

If COMPZ = 'I' and N > 1 then LWORK must be at least

( 1 + 4*N + N**2 ).

Note that for COMPZ = 'I' or 'V', then if N is less than or

equal to the minimum divide size, usually 25, then LWORK need

only be max(1,2*(N-1)).

If LWORK = -1, then a workspace query is assumed; the routine

only calculates the optimal size of the WORK array, returns

this value as the first entry of the WORK array, and no error

message related to LWORK is issued by XERBLA.

*IWORK*

IWORK is INTEGER array, dimension (MAX(1,LIWORK))

On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.

*LIWORK*

LIWORK is INTEGER

The dimension of the array IWORK.

If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.

If COMPZ = 'V' and N > 1 then LIWORK must be at least

( 6 + 6*N + 5*N*lg N ).

If COMPZ = 'I' and N > 1 then LIWORK must be at least

( 3 + 5*N ).

Note that for COMPZ = 'I' or 'V', then if N is less than or

equal to the minimum divide size, usually 25, then LIWORK

need only be 1.

If LIWORK = -1, then a workspace query is assumed; the

routine only calculates the optimal size of the IWORK array,

returns this value as the first entry of the IWORK array, and

no error message related to LIWORK is issued by XERBLA.

*INFO*

INFO is INTEGER

= 0: successful exit.

< 0: if INFO = -i, the i-th argument had an illegal value.

> 0: The algorithm failed to compute an eigenvalue while

working on the submatrix lying in rows and columns

INFO/(N+1) through mod(INFO,N+1).

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

Modified by Francoise Tisseur, University of Tennessee

## subroutine ssteqr (character COMPZ, integer N, real, dimension( * ) D, real, dimension( * ) E, real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer INFO)¶

**SSTEQR**

**Purpose:**

SSTEQR computes all eigenvalues and, optionally, eigenvectors of a

symmetric tridiagonal matrix using the implicit QL or QR method.

The eigenvectors of a full or band symmetric matrix can also be found

if SSYTRD or SSPTRD or SSBTRD has been used to reduce this matrix to

tridiagonal form.

**Parameters**

*COMPZ*

COMPZ is CHARACTER*1

= 'N': Compute eigenvalues only.

= 'V': Compute eigenvalues and eigenvectors of the original

symmetric matrix. On entry, Z must contain the

orthogonal matrix used to reduce the original matrix

to tridiagonal form.

= 'I': Compute eigenvalues and eigenvectors of the

tridiagonal matrix. Z is initialized to the identity

matrix.

*N*

N is INTEGER

The order of the matrix. N >= 0.

*D*

D is REAL array, dimension (N)

On entry, the diagonal elements of the tridiagonal matrix.

On exit, if INFO = 0, the eigenvalues in ascending order.

*E*

E is REAL array, dimension (N-1)

On entry, the (n-1) subdiagonal elements of the tridiagonal

matrix.

On exit, E has been destroyed.

*Z*

Z is REAL array, dimension (LDZ, N)

On entry, if COMPZ = 'V', then Z contains the orthogonal

matrix used in the reduction to tridiagonal form.

On exit, if INFO = 0, then if COMPZ = 'V', Z contains the

orthonormal eigenvectors of the original symmetric matrix,

and if COMPZ = 'I', Z contains the orthonormal eigenvectors

of the symmetric tridiagonal matrix.

If COMPZ = 'N', then Z is not referenced.

*LDZ*

LDZ is INTEGER

The leading dimension of the array Z. LDZ >= 1, and if

eigenvectors are desired, then LDZ >= max(1,N).

*WORK*

WORK is REAL array, dimension (max(1,2*N-2))

If COMPZ = 'N', then WORK is not referenced.

*INFO*

INFO is INTEGER

= 0: successful exit

< 0: if INFO = -i, the i-th argument had an illegal value

> 0: the algorithm has failed to find all the eigenvalues in

a total of 30*N iterations; if INFO = i, then i

elements of E have not converged to zero; on exit, D

and E contain the elements of a symmetric tridiagonal

matrix which is orthogonally similar to the original

matrix.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine ssterf (integer N, real, dimension( * ) D, real, dimension( * ) E, integer INFO)¶

**SSTERF**

**Purpose:**

SSTERF computes all eigenvalues of a symmetric tridiagonal matrix

using the Pal-Walker-Kahan variant of the QL or QR algorithm.

**Parameters**

*N*

N is INTEGER

The order of the matrix. N >= 0.

*D*

D is REAL array, dimension (N)

On entry, the n diagonal elements of the tridiagonal matrix.

On exit, if INFO = 0, the eigenvalues in ascending order.

*E*

E is REAL array, dimension (N-1)

On entry, the (n-1) subdiagonal elements of the tridiagonal

matrix.

On exit, E has been destroyed.

*INFO*

INFO is INTEGER

= 0: successful exit

< 0: if INFO = -i, the i-th argument had an illegal value

> 0: the algorithm failed to find all of the eigenvalues in

a total of 30*N iterations; if INFO = i, then i

elements of E have not converged to zero.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

# Author¶

Generated automatically by Doxygen for LAPACK from the source code.

Wed Sep 1 2021 | Version 3.10.0 |