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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

 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.

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.

  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.

  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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLAED0 computes all eigenvalues and corresponding eigenvectors of a


 symmetric tridiagonal matrix using the divide and conquer method.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

 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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

 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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

 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.

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.

  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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA

 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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  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.

Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA

 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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

 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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

 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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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)

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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)

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLASQ4 computes an approximation TAU to the smallest eigenvalue


 using values of d from the previous transform.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  CNST1 = 9/16

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


 version for IEEE machines another for non IEEE machines.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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


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

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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.

  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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DSTERF computes all eigenvalues of a symmetric tridiagonal matrix


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

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

 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.

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.

  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.

  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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 SLAED0 computes all eigenvalues and corresponding eigenvectors of a


 symmetric tridiagonal matrix using the divide and conquer method.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

 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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

 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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

 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.

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.

  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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA

 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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  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.

Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA

 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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

 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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

 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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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)

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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)

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 SLASQ4 computes an approximation TAU to the smallest eigenvalue


 using values of d from the previous transform.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  CNST1 = 9/16

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


 version for IEEE machines another for non IEEE machines.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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


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

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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.

  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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

 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.

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 SSTERF computes all eigenvalues of a symmetric tridiagonal matrix


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

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.