Provided by: liblapack-doc_3.11.0-2build1_all bug

NAME

       auxOTHERcomputational - auxiliary Computational routines

SYNOPSIS

   Functions
       character *1 function chla_transtype (TRANS)
           CHLA_TRANSTYPE
       subroutine dbdsdc (UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
           DBDSDC
       subroutine dbdsqr (UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
           DBDSQR
       subroutine ddisna (JOB, M, N, D, SEP, INFO)
           DDISNA
       subroutine dlaed0 (ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, WORK, IWORK, INFO)
           DLAED0 used by DSTEDC. Computes all eigenvalues and corresponding eigenvectors of an
           unreduced symmetric tridiagonal matrix using the divide and conquer method.
       subroutine dlaed1 (N, D, Q, LDQ, INDXQ, RHO, CUTPNT, WORK, IWORK, INFO)
           DLAED1 used by DSTEDC. Computes the updated eigensystem of a diagonal matrix after
           modification by a rank-one symmetric matrix. Used when the original matrix is
           tridiagonal.
       subroutine dlaed2 (K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMDA, W, Q2, INDX, INDXC, INDXP,
           COLTYP, INFO)
           DLAED2 used by DSTEDC. Merges eigenvalues and deflates secular equation. Used when the
           original matrix is tridiagonal.
       subroutine dlaed3 (K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX, CTOT, W, S, INFO)
           DLAED3 used by DSTEDC. Finds the roots of the secular equation and updates the
           eigenvectors. Used when the original matrix is tridiagonal.
       subroutine dlaed4 (N, I, D, Z, DELTA, RHO, DLAM, INFO)
           DLAED4 used by DSTEDC. Finds a single root of the secular equation.
       subroutine dlaed5 (I, D, Z, DELTA, RHO, DLAM)
           DLAED5 used by DSTEDC. Solves the 2-by-2 secular equation.
       subroutine dlaed6 (KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO)
           DLAED6 used by DSTEDC. Computes one Newton step in solution of the secular equation.
       subroutine dlaed7 (ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q, LDQ, INDXQ, RHO, CUTPNT,
           QSTORE, QPTR, PRMPTR, PERM, GIVPTR, GIVCOL, GIVNUM, WORK, IWORK, INFO)
           DLAED7 used by DSTEDC. Computes the updated eigensystem of a diagonal matrix after
           modification by a rank-one symmetric matrix. Used when the original matrix is dense.
       subroutine dlaed8 (ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, CUTPNT, Z, DLAMDA, Q2, LDQ2,
           W, PERM, GIVPTR, GIVCOL, GIVNUM, INDXP, INDX, INFO)
           DLAED8 used by DSTEDC. Merges eigenvalues and deflates secular equation. Used when the
           original matrix is dense.
       subroutine dlaed9 (K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W, S, LDS, INFO)
           DLAED9 used by DSTEDC. Finds the roots of the secular equation and updates the
           eigenvectors. Used when the original matrix is dense.
       subroutine dlaeda (N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR, GIVCOL, GIVNUM, Q,
           QPTR, Z, ZTEMP, INFO)
           DLAEDA used by DSTEDC. Computes the Z vector determining the rank-one modification of
           the diagonal matrix. Used when the original matrix is dense.
       subroutine dlagtf (N, A, LAMBDA, B, C, TOL, D, IN, INFO)
           DLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal
           matrix, and λ a scalar, using partial pivoting with row interchanges.
       subroutine dlamrg (N1, N2, A, DTRD1, DTRD2, INDEX)
           DLAMRG creates a permutation list to merge the entries of two independently sorted
           sets into a single set sorted in ascending order.
       subroutine dlartgs (X, Y, SIGMA, CS, SN)
           DLARTGS generates a plane rotation designed to introduce a bulge in implicit QR
           iteration for the bidiagonal SVD problem.
       subroutine dlasq1 (N, D, E, WORK, INFO)
           DLASQ1 computes the singular values of a real square bidiagonal matrix. Used by
           sbdsqr.
       subroutine dlasq2 (N, Z, INFO)
           DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal
           matrix associated with the qd Array Z to high relative accuracy. Used by sbdsqr and
           sstegr.
       subroutine dlasq3 (I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, ITER, NDIV, IEEE,
           TTYPE, DMIN1, DMIN2, DN, DN1, DN2, G, TAU)
           DLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr.
       subroutine dlasq4 (I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1, DN2, TAU, TTYPE, G)
           DLASQ4 computes an approximation to the smallest eigenvalue using values of d from the
           previous transform. Used by sbdsqr.
       subroutine dlasq5 (I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, IEEE,
           EPS)
           DLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.
       subroutine dlasq6 (I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2)
           DLASQ6 computes one dqd transform in ping-pong form. Used by sbdsqr and sstegr.
       subroutine dlasrt (ID, N, D, INFO)
           DLASRT sorts numbers in increasing or decreasing order.
       subroutine dstebz (RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK,
           ISPLIT, WORK, IWORK, INFO)
           DSTEBZ
       subroutine dstedc (COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
           DSTEDC
       subroutine dsteqr (COMPZ, N, D, E, Z, LDZ, WORK, INFO)
           DSTEQR
       subroutine dsterf (N, D, E, INFO)
           DSTERF
       integer function iladiag (DIAG)
           ILADIAG
       integer function ilaprec (PREC)
           ILAPREC
       integer function ilatrans (TRANS)
           ILATRANS
       integer function ilauplo (UPLO)
           ILAUPLO
       subroutine sbdsdc (UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
           SBDSDC
       subroutine sbdsqr (UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
           SBDSQR
       subroutine sdisna (JOB, M, N, D, SEP, INFO)
           SDISNA
       subroutine slaed0 (ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, WORK, IWORK, INFO)
           SLAED0 used by SSTEDC. Computes all eigenvalues and corresponding eigenvectors of an
           unreduced symmetric tridiagonal matrix using the divide and conquer method.
       subroutine slaed1 (N, D, Q, LDQ, INDXQ, RHO, CUTPNT, WORK, IWORK, INFO)
           SLAED1 used by SSTEDC. Computes the updated eigensystem of a diagonal matrix after
           modification by a rank-one symmetric matrix. Used when the original matrix is
           tridiagonal.
       subroutine slaed2 (K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMDA, W, Q2, INDX, INDXC, INDXP,
           COLTYP, INFO)
           SLAED2 used by SSTEDC. Merges eigenvalues and deflates secular equation. Used when the
           original matrix is tridiagonal.
       subroutine slaed3 (K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX, CTOT, W, S, INFO)
           SLAED3 used by SSTEDC. Finds the roots of the secular equation and updates the
           eigenvectors. Used when the original matrix is tridiagonal.
       subroutine slaed4 (N, I, D, Z, DELTA, RHO, DLAM, INFO)
           SLAED4 used by SSTEDC. Finds a single root of the secular equation.
       subroutine slaed5 (I, D, Z, DELTA, RHO, DLAM)
           SLAED5 used by SSTEDC. Solves the 2-by-2 secular equation.
       subroutine slaed6 (KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO)
           SLAED6 used by SSTEDC. Computes one Newton step in solution of the secular equation.
       subroutine slaed7 (ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q, LDQ, INDXQ, RHO, CUTPNT,
           QSTORE, QPTR, PRMPTR, PERM, GIVPTR, GIVCOL, GIVNUM, WORK, IWORK, INFO)
           SLAED7 used by SSTEDC. Computes the updated eigensystem of a diagonal matrix after
           modification by a rank-one symmetric matrix. Used when the original matrix is dense.
       subroutine slaed8 (ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, CUTPNT, Z, DLAMDA, Q2, LDQ2,
           W, PERM, GIVPTR, GIVCOL, GIVNUM, INDXP, INDX, INFO)
           SLAED8 used by SSTEDC. Merges eigenvalues and deflates secular equation. Used when the
           original matrix is dense.
       subroutine slaed9 (K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W, S, LDS, INFO)
           SLAED9 used by SSTEDC. Finds the roots of the secular equation and updates the
           eigenvectors. Used when the original matrix is dense.
       subroutine slaeda (N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR, GIVCOL, GIVNUM, Q,
           QPTR, Z, ZTEMP, INFO)
           SLAEDA used by SSTEDC. Computes the Z vector determining the rank-one modification of
           the diagonal matrix. Used when the original matrix is dense.
       subroutine slagtf (N, A, LAMBDA, B, C, TOL, D, IN, INFO)
           SLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal
           matrix, and λ a scalar, using partial pivoting with row interchanges.
       subroutine slamrg (N1, N2, A, STRD1, STRD2, INDEX)
           SLAMRG creates a permutation list to merge the entries of two independently sorted
           sets into a single set sorted in ascending order.
       subroutine slartgs (X, Y, SIGMA, CS, SN)
           SLARTGS generates a plane rotation designed to introduce a bulge in implicit QR
           iteration for the bidiagonal SVD problem.
       subroutine slasq1 (N, D, E, WORK, INFO)
           SLASQ1 computes the singular values of a real square bidiagonal matrix. Used by
           sbdsqr.
       subroutine slasq2 (N, Z, INFO)
           SLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal
           matrix associated with the qd Array Z to high relative accuracy. Used by sbdsqr and
           sstegr.
       subroutine slasq3 (I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, ITER, NDIV, IEEE,
           TTYPE, DMIN1, DMIN2, DN, DN1, DN2, G, TAU)
           SLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr.
       subroutine slasq4 (I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1, DN2, TAU, TTYPE, G)
           SLASQ4 computes an approximation to the smallest eigenvalue using values of d from the
           previous transform. Used by sbdsqr.
       subroutine slasq5 (I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, IEEE,
           EPS)
            SLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.
       subroutine slasq6 (I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2)
           SLASQ6 computes one dqd transform in ping-pong form. Used by sbdsqr and sstegr.
       subroutine slasrt (ID, N, D, INFO)
           SLASRT sorts numbers in increasing or decreasing order.
       subroutine spttrf (N, D, E, INFO)
           SPTTRF
       subroutine sstebz (RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK,
           ISPLIT, WORK, IWORK, INFO)
           SSTEBZ
       subroutine sstedc (COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
           SSTEDC
       subroutine ssteqr (COMPZ, N, D, E, Z, LDZ, WORK, INFO)
           SSTEQR
       subroutine ssterf (N, D, E, INFO)
           SSTERF

Detailed Description

       This is the group of auxiliary Computational routines

Function Documentation

   character*1 function chla_transtype (integer TRANS)
       CHLA_TRANSTYPE

       Purpose:

            This subroutine translates from a BLAST-specified integer constant to
            the character string specifying a transposition operation.

            CHLA_TRANSTYPE returns an CHARACTER*1.  If CHLA_TRANSTYPE is 'X',
            then input is not an integer indicating a transposition operator.
            Otherwise CHLA_TRANSTYPE returns the constant value corresponding to
            TRANS.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

       Purpose:

            DBDSDC computes the singular value decomposition (SVD) of a real
            N-by-N (upper or lower) bidiagonal matrix B:  B = U * S * VT,
            using a divide and conquer method, where S is a diagonal matrix
            with non-negative diagonal elements (the singular values of B), and
            U and VT are orthogonal matrices of left and right singular vectors,
            respectively. DBDSDC can be used to compute all singular values,
            and optionally, singular vectors or singular vectors in compact form.

            This code makes very mild assumptions about floating point
            arithmetic. It will work on machines with a guard digit in
            add/subtract, or on those binary machines without guard digits
            which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
            It could conceivably fail on hexadecimal or decimal machines
            without guard digits, but we know of none.  See DLASD3 for details.

            The code currently calls DLASDQ if singular values only are desired.
            However, it can be slightly modified to compute singular values
            using the divide and conquer method.

       Parameters
           UPLO

                     UPLO is CHARACTER*1
                     = 'U':  B is upper bidiagonal.
                     = 'L':  B is lower bidiagonal.

           COMPQ

                     COMPQ is CHARACTER*1
                     Specifies whether singular vectors are to be computed
                     as follows:
                     = 'N':  Compute singular values only;
                     = 'P':  Compute singular values and compute singular
                             vectors in compact form;
                     = 'I':  Compute singular values and singular vectors.

           N

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

           D

                     D is DOUBLE PRECISION array, dimension (N)
                     On entry, the n diagonal elements of the bidiagonal matrix B.
                     On exit, if INFO=0, the singular values of B.

           E

                     E is DOUBLE PRECISION array, dimension (N-1)
                     On entry, the elements of E contain the offdiagonal
                     elements of the bidiagonal matrix whose SVD is desired.
                     On exit, E has been destroyed.

           U

                     U is DOUBLE PRECISION array, dimension (LDU,N)
                     If  COMPQ = 'I', then:
                        On exit, if INFO = 0, U contains the left singular vectors
                        of the bidiagonal matrix.
                     For other values of COMPQ, U is not referenced.

           LDU

                     LDU is INTEGER
                     The leading dimension of the array U.  LDU >= 1.
                     If singular vectors are desired, then LDU >= max( 1, N ).

           VT

                     VT is DOUBLE PRECISION array, dimension (LDVT,N)
                     If  COMPQ = 'I', then:
                        On exit, if INFO = 0, VT**T contains the right singular
                        vectors of the bidiagonal matrix.
                     For other values of COMPQ, VT is not referenced.

           LDVT

                     LDVT is INTEGER
                     The leading dimension of the array VT.  LDVT >= 1.
                     If singular vectors are desired, then LDVT >= max( 1, N ).

           Q

                     Q is DOUBLE PRECISION array, dimension (LDQ)
                     If  COMPQ = 'P', then:
                        On exit, if INFO = 0, Q and IQ contain the left
                        and right singular vectors in a compact form,
                        requiring O(N log N) space instead of 2*N**2.
                        In particular, Q contains all the DOUBLE PRECISION data in
                        LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
                        words of memory, where SMLSIZ is returned by ILAENV and
                        is equal to the maximum size of the subproblems at the
                        bottom of the computation tree (usually about 25).
                     For other values of COMPQ, Q is not referenced.

           IQ

                     IQ is INTEGER array, dimension (LDIQ)
                     If  COMPQ = 'P', then:
                        On exit, if INFO = 0, Q and IQ contain the left
                        and right singular vectors in a compact form,
                        requiring O(N log N) space instead of 2*N**2.
                        In particular, IQ contains all INTEGER data in
                        LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
                        words of memory, where SMLSIZ is returned by ILAENV and
                        is equal to the maximum size of the subproblems at the
                        bottom of the computation tree (usually about 25).
                     For other values of COMPQ, IQ is not referenced.

           WORK

                     WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
                     If COMPQ = 'N' then LWORK >= (4 * N).
                     If COMPQ = 'P' then LWORK >= (6 * N).
                     If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).

           IWORK

                     IWORK is INTEGER array, dimension (8*N)

           INFO

                     INFO is INTEGER
                     = 0:  successful exit.
                     < 0:  if INFO = -i, the i-th argument had an illegal value.
                     > 0:  The algorithm failed to compute a singular value.
                           The update process of divide and conquer failed.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

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

       Purpose:

            DBDSQR computes the singular values and, optionally, the right and/or
            left singular vectors from the singular value decomposition (SVD) of
            a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
            zero-shift QR algorithm.  The SVD of B has the form

               B = Q * S * P**T

            where S is the diagonal matrix of singular values, Q is an orthogonal
            matrix of left singular vectors, and P is an orthogonal matrix of
            right singular vectors.  If left singular vectors are requested, this
            subroutine actually returns U*Q instead of Q, and, if right singular
            vectors are requested, this subroutine returns P**T*VT instead of
            P**T, for given real input matrices U and VT.  When U and VT are the
            orthogonal matrices that reduce a general matrix A to bidiagonal
            form:  A = U*B*VT, as computed by DGEBRD, then

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

            is the SVD of A.  Optionally, the subroutine may also compute Q**T*C
            for a given real input matrix C.

            See 'Computing  Small Singular Values of Bidiagonal Matrices With
            Guaranteed High Relative Accuracy,' by J. Demmel and W. Kahan,
            LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
            no. 5, pp. 873-912, Sept 1990) and
            'Accurate singular values and differential qd algorithms,' by
            B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
            Department, University of California at Berkeley, July 1992
            for a detailed description of the algorithm.

       Parameters
           UPLO

                     UPLO is CHARACTER*1
                     = 'U':  B is upper bidiagonal;
                     = 'L':  B is lower bidiagonal.

           N

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

           NCVT

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

           NRU

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

           NCC

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

           D

                     D is DOUBLE PRECISION array, dimension (N)
                     On entry, the n diagonal elements of the bidiagonal matrix B.
                     On exit, if INFO=0, the singular values of B in decreasing
                     order.

           E

                     E is DOUBLE PRECISION array, dimension (N-1)
                     On entry, the N-1 offdiagonal elements of the bidiagonal
                     matrix B.
                     On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
                     will contain the diagonal and superdiagonal elements of a
                     bidiagonal matrix orthogonally equivalent to the one given
                     as input.

           VT

                     VT is DOUBLE PRECISION array, dimension (LDVT, NCVT)
                     On entry, an N-by-NCVT matrix VT.
                     On exit, VT is overwritten by P**T * VT.
                     Not referenced if NCVT = 0.

           LDVT

                     LDVT is INTEGER
                     The leading dimension of the array VT.
                     LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.

           U

                     U is DOUBLE PRECISION array, dimension (LDU, N)
                     On entry, an NRU-by-N matrix U.
                     On exit, U is overwritten by U * Q.
                     Not referenced if NRU = 0.

           LDU

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

           C

                     C is DOUBLE PRECISION array, dimension (LDC, NCC)
                     On entry, an N-by-NCC matrix C.
                     On exit, C is overwritten by Q**T * C.
                     Not referenced if NCC = 0.

           LDC

                     LDC is INTEGER
                     The leading dimension of the array C.
                     LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.

           WORK

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

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  If INFO = -i, the i-th argument had an illegal value
                     > 0:
                        if NCVT = NRU = NCC = 0,
                           = 1, a split was marked by a positive value in E
                           = 2, current block of Z not diagonalized after 30*N
                                iterations (in inner while loop)
                           = 3, termination criterion of outer while loop not met
                                (program created more than N unreduced blocks)
                        else NCVT = NRU = NCC = 0,
                              the algorithm did not converge; D and E contain the
                              elements of a bidiagonal matrix which is orthogonally
                              similar to the input matrix B;  if INFO = i, i
                              elements of E have not converged to zero.

       Internal Parameters:

             TOLMUL  DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))
                     TOLMUL controls the convergence criterion of the QR loop.
                     If it is positive, TOLMUL*EPS is the desired relative
                        precision in the computed singular values.
                     If it is negative, abs(TOLMUL*EPS*sigma_max) is the
                        desired absolute accuracy in the computed singular
                        values (corresponds to relative accuracy
                        abs(TOLMUL*EPS) in the largest singular value.
                     abs(TOLMUL) should be between 1 and 1/EPS, and preferably
                        between 10 (for fast convergence) and .1/EPS
                        (for there to be some accuracy in the results).
                     Default is to lose at either one eighth or 2 of the
                        available decimal digits in each computed singular value
                        (whichever is smaller).

             MAXITR  INTEGER, default = 6
                     MAXITR controls the maximum number of passes of the
                     algorithm through its inner loop. The algorithms stops
                     (and so fails to converge) if the number of passes
                     through the inner loop exceeds MAXITR*N**2.

       Note:

             Bug report from Cezary Dendek.
             On March 23rd 2017, the INTEGER variable MAXIT = MAXITR*N**2 is
             removed since it can overflow pretty easily (for N larger or equal
             than 18,919). We instead use MAXITDIVN = MAXITR*N.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

       Purpose:

            DDISNA computes the reciprocal condition numbers for the eigenvectors
            of a real symmetric or complex Hermitian matrix or for the left or
            right singular vectors of a general m-by-n matrix. The reciprocal
            condition number is the 'gap' between the corresponding eigenvalue or
            singular value and the nearest other one.

            The bound on the error, measured by angle in radians, in the I-th
            computed vector is given by

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

            where ANORM = 2-norm(A) = max( abs( D(j) ) ).  SEP(I) is not allowed
            to be smaller than DLAMCH( 'E' )*ANORM in order to limit the size of
            the error bound.

            DDISNA may also be used to compute error bounds for eigenvectors of
            the generalized symmetric definite eigenproblem.

       Parameters
           JOB

                     JOB is CHARACTER*1
                     Specifies for which problem the reciprocal condition numbers
                     should be computed:
                     = 'E':  the eigenvectors of a symmetric/Hermitian matrix;
                     = 'L':  the left singular vectors of a general matrix;
                     = 'R':  the right singular vectors of a general matrix.

           M

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

           N

                     N is INTEGER
                     If JOB = 'L' or 'R', the number of columns of the matrix,
                     in which case N >= 0. Ignored if JOB = 'E'.

           D

                     D is DOUBLE PRECISION array, dimension (M) if JOB = 'E'
                                         dimension (min(M,N)) if JOB = 'L' or 'R'
                     The eigenvalues (if JOB = 'E') or singular values (if JOB =
                     'L' or 'R') of the matrix, in either increasing or decreasing
                     order. If singular values, they must be non-negative.

           SEP

                     SEP is DOUBLE PRECISION array, dimension (M) if JOB = 'E'
                                          dimension (min(M,N)) if JOB = 'L' or 'R'
                     The reciprocal condition numbers of the vectors.

           INFO

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine dlaed0 (integer ICOMPQ, integer QSIZ, integer N, double precision, dimension( * )
       D, double precision, dimension( * ) E, double precision, dimension( ldq, * ) Q, integer
       LDQ, double precision, dimension( ldqs, * ) QSTORE, integer LDQS, double precision,
       dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)
       DLAED0 used by DSTEDC. Computes all eigenvalues and corresponding eigenvectors of an
       unreduced symmetric tridiagonal matrix using the divide and conquer method.

       Purpose:

            DLAED0 computes all eigenvalues and corresponding eigenvectors of a
            symmetric tridiagonal matrix using the divide and conquer method.

       Parameters
           ICOMPQ

                     ICOMPQ is INTEGER
                     = 0:  Compute eigenvalues only.
                     = 1:  Compute eigenvectors of original dense symmetric matrix
                           also.  On entry, Q contains the orthogonal matrix used
                           to reduce the original matrix to tridiagonal form.
                     = 2:  Compute eigenvalues and eigenvectors of tridiagonal
                           matrix.

           QSIZ

                     QSIZ is INTEGER
                    The dimension of the orthogonal matrix used to reduce
                    the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.

           N

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

           D

                     D is DOUBLE PRECISION array, dimension (N)
                    On entry, the main diagonal of the tridiagonal matrix.
                    On exit, its eigenvalues.

           E

                     E is DOUBLE PRECISION array, dimension (N-1)
                    The off-diagonal elements of the tridiagonal matrix.
                    On exit, E has been destroyed.

           Q

                     Q is DOUBLE PRECISION array, dimension (LDQ, N)
                    On entry, Q must contain an N-by-N orthogonal matrix.
                    If ICOMPQ = 0    Q is not referenced.
                    If ICOMPQ = 1    On entry, Q is a subset of the columns of the
                                     orthogonal matrix used to reduce the full
                                     matrix to tridiagonal form corresponding to
                                     the subset of the full matrix which is being
                                     decomposed at this time.
                    If ICOMPQ = 2    On entry, Q will be the identity matrix.
                                     On exit, Q contains the eigenvectors of the
                                     tridiagonal matrix.

           LDQ

                     LDQ is INTEGER
                    The leading dimension of the array Q.  If eigenvectors are
                    desired, then  LDQ >= max(1,N).  In any case,  LDQ >= 1.

           QSTORE

                     QSTORE is DOUBLE PRECISION array, dimension (LDQS, N)
                    Referenced only when ICOMPQ = 1.  Used to store parts of
                    the eigenvector matrix when the updating matrix multiplies
                    take place.

           LDQS

                     LDQS is INTEGER
                    The leading dimension of the array QSTORE.  If ICOMPQ = 1,
                    then  LDQS >= max(1,N).  In any case,  LDQS >= 1.

           WORK

                     WORK is DOUBLE PRECISION array,
                    If ICOMPQ = 0 or 1, the dimension of WORK must be at least
                                1 + 3*N + 2*N*lg N + 3*N**2
                                ( lg( N ) = smallest integer k
                                            such that 2^k >= N )
                    If ICOMPQ = 2, the dimension of WORK must be at least
                                4*N + N**2.

           IWORK

                     IWORK is INTEGER array,
                    If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
                                   6 + 6*N + 5*N*lg N.
                                   ( lg( N ) = smallest integer k
                                               such that 2^k >= N )
                    If ICOMPQ = 2, the dimension of IWORK must be at least
                                   3 + 5*N.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit.
                     < 0:  if INFO = -i, the i-th argument had an illegal value.
                     > 0:  The algorithm failed to compute an eigenvalue while
                           working on the submatrix lying in rows and columns
                           INFO/(N+1) through mod(INFO,N+1).

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine dlaed1 (integer N, double precision, dimension( * ) D, double precision, dimension(
       ldq, * ) Q, integer LDQ, integer, dimension( * ) INDXQ, double precision RHO, integer
       CUTPNT, double precision, dimension( * ) WORK, integer, dimension( * ) IWORK, integer
       INFO)
       DLAED1 used by DSTEDC. Computes the updated eigensystem of a diagonal matrix after
       modification by a rank-one symmetric matrix. Used when the original matrix is tridiagonal.

       Purpose:

            DLAED1 computes the updated eigensystem of a diagonal
            matrix after modification by a rank-one symmetric matrix.  This
            routine is used only for the eigenproblem which requires all
            eigenvalues and eigenvectors of a tridiagonal matrix.  DLAED7 handles
            the case in which eigenvalues only or eigenvalues and eigenvectors
            of a full symmetric matrix (which was reduced to tridiagonal form)
            are desired.

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

               where Z = Q**T*u, u is a vector of length N with ones in the
               CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.

               The eigenvectors of the original matrix are stored in Q, and the
               eigenvalues are in D.  The algorithm consists of three stages:

                  The first stage consists of deflating the size of the problem
                  when there are multiple eigenvalues or if there is a zero in
                  the Z vector.  For each such occurrence the dimension of the
                  secular equation problem is reduced by one.  This stage is
                  performed by the routine DLAED2.

                  The second stage consists of calculating the updated
                  eigenvalues. This is done by finding the roots of the secular
                  equation via the routine DLAED4 (as called by DLAED3).
                  This routine also calculates the eigenvectors of the current
                  problem.

                  The final stage consists of computing the updated eigenvectors
                  directly using the updated eigenvalues.  The eigenvectors for
                  the current problem are multiplied with the eigenvectors from
                  the overall problem.

       Parameters
           N

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

           D

                     D is DOUBLE PRECISION array, dimension (N)
                    On entry, the eigenvalues of the rank-1-perturbed matrix.
                    On exit, the eigenvalues of the repaired matrix.

           Q

                     Q is DOUBLE PRECISION array, dimension (LDQ,N)
                    On entry, the eigenvectors of the rank-1-perturbed matrix.
                    On exit, the eigenvectors of the repaired tridiagonal matrix.

           LDQ

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

           INDXQ

                     INDXQ is INTEGER array, dimension (N)
                    On entry, the permutation which separately sorts the two
                    subproblems in D into ascending order.
                    On exit, the permutation which will reintegrate the
                    subproblems back into sorted order,
                    i.e. D( INDXQ( I = 1, N ) ) will be in ascending order.

           RHO

                     RHO is DOUBLE PRECISION
                    The subdiagonal entry used to create the rank-1 modification.

           CUTPNT

                     CUTPNT is INTEGER
                    The location of the last eigenvalue in the leading sub-matrix.
                    min(1,N) <= CUTPNT <= N/2.

           WORK

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

           IWORK

                     IWORK is INTEGER array, dimension (4*N)

           INFO

                     INFO is INTEGER
                     = 0:  successful exit.
                     < 0:  if INFO = -i, the i-th argument had an illegal value.
                     > 0:  if INFO = 1, an eigenvalue did not converge

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine dlaed2 (integer K, integer N, integer N1, double precision, dimension( * ) D,
       double precision, dimension( ldq, * ) Q, integer LDQ, integer, dimension( * ) INDXQ,
       double precision RHO, double precision, dimension( * ) Z, double precision, dimension( * )
       DLAMDA, double precision, dimension( * ) W, double precision, dimension( * ) Q2, integer,
       dimension( * ) INDX, integer, dimension( * ) INDXC, integer, dimension( * ) INDXP,
       integer, dimension( * ) COLTYP, integer INFO)
       DLAED2 used by DSTEDC. Merges eigenvalues and deflates secular equation. Used when the
       original matrix is tridiagonal.

       Purpose:

            DLAED2 merges the two sets of eigenvalues together into a single
            sorted set.  Then it tries to deflate the size of the problem.
            There are two ways in which deflation can occur:  when two or more
            eigenvalues are close together or if there is a tiny entry in the
            Z vector.  For each such occurrence the order of the related secular
            equation problem is reduced by one.

       Parameters
           K

                     K is INTEGER
                    The number of non-deflated eigenvalues, and the order of the
                    related secular equation. 0 <= K <=N.

           N

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

           N1

                     N1 is INTEGER
                    The location of the last eigenvalue in the leading sub-matrix.
                    min(1,N) <= N1 <= N/2.

           D

                     D is DOUBLE PRECISION array, dimension (N)
                    On entry, D contains the eigenvalues of the two submatrices to
                    be combined.
                    On exit, D contains the trailing (N-K) updated eigenvalues
                    (those which were deflated) sorted into increasing order.

           Q

                     Q is DOUBLE PRECISION array, dimension (LDQ, N)
                    On entry, Q contains the eigenvectors of two submatrices in
                    the two square blocks with corners at (1,1), (N1,N1)
                    and (N1+1, N1+1), (N,N).
                    On exit, Q contains the trailing (N-K) updated eigenvectors
                    (those which were deflated) in its last N-K columns.

           LDQ

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

           INDXQ

                     INDXQ is INTEGER array, dimension (N)
                    The permutation which separately sorts the two sub-problems
                    in D into ascending order.  Note that elements in the second
                    half of this permutation must first have N1 added to their
                    values. Destroyed on exit.

           RHO

                     RHO is DOUBLE PRECISION
                    On entry, the off-diagonal element associated with the rank-1
                    cut which originally split the two submatrices which are now
                    being recombined.
                    On exit, RHO has been modified to the value required by
                    DLAED3.

           Z

                     Z is DOUBLE PRECISION array, dimension (N)
                    On entry, Z contains the updating vector (the last
                    row of the first sub-eigenvector matrix and the first row of
                    the second sub-eigenvector matrix).
                    On exit, the contents of Z have been destroyed by the updating
                    process.

           DLAMDA

                     DLAMDA is DOUBLE PRECISION array, dimension (N)
                    A copy of the first K eigenvalues which will be used by
                    DLAED3 to form the secular equation.

           W

                     W is DOUBLE PRECISION array, dimension (N)
                    The first k values of the final deflation-altered z-vector
                    which will be passed to DLAED3.

           Q2

                     Q2 is DOUBLE PRECISION array, dimension (N1**2+(N-N1)**2)
                    A copy of the first K eigenvectors which will be used by
                    DLAED3 in a matrix multiply (DGEMM) to solve for the new
                    eigenvectors.

           INDX

                     INDX is INTEGER array, dimension (N)
                    The permutation used to sort the contents of DLAMDA into
                    ascending order.

           INDXC

                     INDXC is INTEGER array, dimension (N)
                    The permutation used to arrange the columns of the deflated
                    Q matrix into three groups:  the first group contains non-zero
                    elements only at and above N1, the second contains
                    non-zero elements only below N1, and the third is dense.

           INDXP

                     INDXP is INTEGER array, dimension (N)
                    The permutation used to place deflated values of D at the end
                    of the array.  INDXP(1:K) points to the nondeflated D-values
                    and INDXP(K+1:N) points to the deflated eigenvalues.

           COLTYP

                     COLTYP is INTEGER array, dimension (N)
                    During execution, a label which will indicate which of the
                    following types a column in the Q2 matrix is:
                    1 : non-zero in the upper half only;
                    2 : dense;
                    3 : non-zero in the lower half only;
                    4 : deflated.
                    On exit, COLTYP(i) is the number of columns of type i,
                    for i=1 to 4 only.

           INFO

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine dlaed3 (integer K, integer N, integer N1, double precision, dimension( * ) D,
       double precision, dimension( ldq, * ) Q, integer LDQ, double precision RHO, double
       precision, dimension( * ) DLAMDA, double precision, dimension( * ) Q2, integer, dimension(
       * ) INDX, integer, dimension( * ) CTOT, double precision, dimension( * ) W, double
       precision, dimension( * ) S, integer INFO)
       DLAED3 used by DSTEDC. Finds the roots of the secular equation and updates the
       eigenvectors. Used when the original matrix is tridiagonal.

       Purpose:

            DLAED3 finds the roots of the secular equation, as defined by the
            values in D, W, and RHO, between 1 and K.  It makes the
            appropriate calls to DLAED4 and then updates the eigenvectors by
            multiplying the matrix of eigenvectors of the pair of eigensystems
            being combined by the matrix of eigenvectors of the K-by-K system
            which is solved here.

            This code makes very mild assumptions about floating point
            arithmetic. It will work on machines with a guard digit in
            add/subtract, or on those binary machines without guard digits
            which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
            It could conceivably fail on hexadecimal or decimal machines
            without guard digits, but we know of none.

       Parameters
           K

                     K is INTEGER
                     The number of terms in the rational function to be solved by
                     DLAED4.  K >= 0.

           N

                     N is INTEGER
                     The number of rows and columns in the Q matrix.
                     N >= K (deflation may result in N>K).

           N1

                     N1 is INTEGER
                     The location of the last eigenvalue in the leading submatrix.
                     min(1,N) <= N1 <= N/2.

           D

                     D is DOUBLE PRECISION array, dimension (N)
                     D(I) contains the updated eigenvalues for
                     1 <= I <= K.

           Q

                     Q is DOUBLE PRECISION array, dimension (LDQ,N)
                     Initially the first K columns are used as workspace.
                     On output the columns 1 to K contain
                     the updated eigenvectors.

           LDQ

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

           RHO

                     RHO is DOUBLE PRECISION
                     The value of the parameter in the rank one update equation.
                     RHO >= 0 required.

           DLAMDA

                     DLAMDA is DOUBLE PRECISION array, dimension (K)
                     The first K elements of this array contain the old roots
                     of the deflated updating problem.  These are the poles
                     of the secular equation. May be changed on output by
                     having lowest order bit set to zero on Cray X-MP, Cray Y-MP,
                     Cray-2, or Cray C-90, as described above.

           Q2

                     Q2 is DOUBLE PRECISION array, dimension (LDQ2*N)
                     The first K columns of this matrix contain the non-deflated
                     eigenvectors for the split problem.

           INDX

                     INDX is INTEGER array, dimension (N)
                     The permutation used to arrange the columns of the deflated
                     Q matrix into three groups (see DLAED2).
                     The rows of the eigenvectors found by DLAED4 must be likewise
                     permuted before the matrix multiply can take place.

           CTOT

                     CTOT is INTEGER array, dimension (4)
                     A count of the total number of the various types of columns
                     in Q, as described in INDX.  The fourth column type is any
                     column which has been deflated.

           W

                     W is DOUBLE PRECISION array, dimension (K)
                     The first K elements of this array contain the components
                     of the deflation-adjusted updating vector. Destroyed on
                     output.

           S

                     S is DOUBLE PRECISION array, dimension (N1 + 1)*K
                     Will contain the eigenvectors of the repaired matrix which
                     will be multiplied by the previously accumulated eigenvectors
                     to update the system.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit.
                     < 0:  if INFO = -i, the i-th argument had an illegal value.
                     > 0:  if INFO = 1, an eigenvalue did not converge

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine dlaed4 (integer N, integer I, double precision, dimension( * ) D, double precision,
       dimension( * ) Z, double precision, dimension( * ) DELTA, double precision RHO, double
       precision DLAM, integer INFO)
       DLAED4 used by DSTEDC. Finds a single root of the secular equation.

       Purpose:

            This subroutine computes the I-th updated eigenvalue of a symmetric
            rank-one modification to a diagonal matrix whose elements are
            given in the array d, and that

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

            and that RHO > 0.  This is arranged by the calling routine, and is
            no loss in generality.  The rank-one modified system is thus

                       diag( D )  +  RHO * Z * Z_transpose.

            where we assume the Euclidean norm of Z is 1.

            The method consists of approximating the rational functions in the
            secular equation by simpler interpolating rational functions.

       Parameters
           N

                     N is INTEGER
                    The length of all arrays.

           I

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

           D

                     D is DOUBLE PRECISION array, dimension (N)
                    The original eigenvalues.  It is assumed that they are in
                    order, D(I) < D(J)  for I < J.

           Z

                     Z is DOUBLE PRECISION array, dimension (N)
                    The components of the updating vector.

           DELTA

                     DELTA is DOUBLE PRECISION array, dimension (N)
                    If N > 2, DELTA contains (D(j) - lambda_I) in its  j-th
                    component.  If N = 1, then DELTA(1) = 1. If N = 2, see DLAED5
                    for detail. The vector DELTA contains the information necessary
                    to construct the eigenvectors by DLAED3 and DLAED9.

           RHO

                     RHO is DOUBLE PRECISION
                    The scalar in the symmetric updating formula.

           DLAM

                     DLAM is DOUBLE PRECISION
                    The computed lambda_I, the I-th updated eigenvalue.

           INFO

                     INFO is INTEGER
                    = 0:  successful exit
                    > 0:  if INFO = 1, the updating process failed.

       Internal Parameters:

             Logical variable ORGATI (origin-at-i?) is used for distinguishing
             whether D(i) or D(i+1) is treated as the origin.

                       ORGATI = .true.    origin at i
                       ORGATI = .false.   origin at i+1

              Logical variable SWTCH3 (switch-for-3-poles?) is for noting
              if we are working with THREE poles!

              MAXIT is the maximum number of iterations allowed for each
              eigenvalue.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine dlaed5 (integer I, double precision, dimension( 2 ) D, double precision, dimension(
       2 ) Z, double precision, dimension( 2 ) DELTA, double precision RHO, double precision
       DLAM)
       DLAED5 used by DSTEDC. Solves the 2-by-2 secular equation.

       Purpose:

            This subroutine computes the I-th eigenvalue of a symmetric rank-one
            modification of a 2-by-2 diagonal matrix

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

            The diagonal elements in the array D are assumed to satisfy

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

            We also assume RHO > 0 and that the Euclidean norm of the vector
            Z is one.

       Parameters
           I

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

           D

                     D is DOUBLE PRECISION array, dimension (2)
                    The original eigenvalues.  We assume D(1) < D(2).

           Z

                     Z is DOUBLE PRECISION array, dimension (2)
                    The components of the updating vector.

           DELTA

                     DELTA is DOUBLE PRECISION array, dimension (2)
                    The vector DELTA contains the information necessary
                    to construct the eigenvectors.

           RHO

                     RHO is DOUBLE PRECISION
                    The scalar in the symmetric updating formula.

           DLAM

                     DLAM is DOUBLE PRECISION
                    The computed lambda_I, the I-th updated eigenvalue.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine dlaed6 (integer KNITER, logical ORGATI, double precision RHO, double precision,
       dimension( 3 ) D, double precision, dimension( 3 ) Z, double precision FINIT, double
       precision TAU, integer INFO)
       DLAED6 used by DSTEDC. Computes one Newton step in solution of the secular equation.

       Purpose:

            DLAED6 computes the positive or negative root (closest to the origin)
            of
                             z(1)        z(2)        z(3)
            f(x) =   rho + --------- + ---------- + ---------
                            d(1)-x      d(2)-x      d(3)-x

            It is assumed that

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

            This routine will be called by DLAED4 when necessary. In most cases,
            the root sought is the smallest in magnitude, though it might not be
            in some extremely rare situations.

       Parameters
           KNITER

                     KNITER is INTEGER
                          Refer to DLAED4 for its significance.

           ORGATI

                     ORGATI is LOGICAL
                          If ORGATI is true, the needed root is between d(2) and
                          d(3); otherwise it is between d(1) and d(2).  See
                          DLAED4 for further details.

           RHO

                     RHO is DOUBLE PRECISION
                          Refer to the equation f(x) above.

           D

                     D is DOUBLE PRECISION array, dimension (3)
                          D satisfies d(1) < d(2) < d(3).

           Z

                     Z is DOUBLE PRECISION array, dimension (3)
                          Each of the elements in z must be positive.

           FINIT

                     FINIT is DOUBLE PRECISION
                          The value of f at 0. It is more accurate than the one
                          evaluated inside this routine (if someone wants to do
                          so).

           TAU

                     TAU is DOUBLE PRECISION
                          The root of the equation f(x).

           INFO

                     INFO is INTEGER
                          = 0: successful exit
                          > 0: if INFO = 1, failure to converge

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Further Details:

             10/02/03: This version has a few statements commented out for thread
             safety (machine parameters are computed on each entry). SJH.

             05/10/06: Modified from a new version of Ren-Cang Li, use
                Gragg-Thornton-Warner cubic convergent scheme for better stability.

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

   subroutine dlaed7 (integer ICOMPQ, integer N, integer QSIZ, integer TLVLS, integer CURLVL,
       integer CURPBM, double precision, dimension( * ) D, double precision, dimension( ldq, * )
       Q, integer LDQ, integer, dimension( * ) INDXQ, double precision RHO, integer CUTPNT,
       double precision, dimension( * ) QSTORE, integer, dimension( * ) QPTR, integer, dimension(
       * ) PRMPTR, integer, dimension( * ) PERM, integer, dimension( * ) GIVPTR, integer,
       dimension( 2, * ) GIVCOL, double precision, dimension( 2, * ) GIVNUM, double precision,
       dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)
       DLAED7 used by DSTEDC. Computes the updated eigensystem of a diagonal matrix after
       modification by a rank-one symmetric matrix. Used when the original matrix is dense.

       Purpose:

            DLAED7 computes the updated eigensystem of a diagonal
            matrix after modification by a rank-one symmetric matrix. This
            routine is used only for the eigenproblem which requires all
            eigenvalues and optionally eigenvectors of a dense symmetric matrix
            that has been reduced to tridiagonal form.  DLAED1 handles
            the case in which all eigenvalues and eigenvectors of a symmetric
            tridiagonal matrix are desired.

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

               where Z = Q**Tu, u is a vector of length N with ones in the
               CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.

               The eigenvectors of the original matrix are stored in Q, and the
               eigenvalues are in D.  The algorithm consists of three stages:

                  The first stage consists of deflating the size of the problem
                  when there are multiple eigenvalues or if there is a zero in
                  the Z vector.  For each such occurrence the dimension of the
                  secular equation problem is reduced by one.  This stage is
                  performed by the routine DLAED8.

                  The second stage consists of calculating the updated
                  eigenvalues. This is done by finding the roots of the secular
                  equation via the routine DLAED4 (as called by DLAED9).
                  This routine also calculates the eigenvectors of the current
                  problem.

                  The final stage consists of computing the updated eigenvectors
                  directly using the updated eigenvalues.  The eigenvectors for
                  the current problem are multiplied with the eigenvectors from
                  the overall problem.

       Parameters
           ICOMPQ

                     ICOMPQ is INTEGER
                     = 0:  Compute eigenvalues only.
                     = 1:  Compute eigenvectors of original dense symmetric matrix
                           also.  On entry, Q contains the orthogonal matrix used
                           to reduce the original matrix to tridiagonal form.

           N

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

           QSIZ

                     QSIZ is INTEGER
                    The dimension of the orthogonal matrix used to reduce
                    the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.

           TLVLS

                     TLVLS is INTEGER
                    The total number of merging levels in the overall divide and
                    conquer tree.

           CURLVL

                     CURLVL is INTEGER
                    The current level in the overall merge routine,
                    0 <= CURLVL <= TLVLS.

           CURPBM

                     CURPBM is INTEGER
                    The current problem in the current level in the overall
                    merge routine (counting from upper left to lower right).

           D

                     D is DOUBLE PRECISION array, dimension (N)
                    On entry, the eigenvalues of the rank-1-perturbed matrix.
                    On exit, the eigenvalues of the repaired matrix.

           Q

                     Q is DOUBLE PRECISION array, dimension (LDQ, N)
                    On entry, the eigenvectors of the rank-1-perturbed matrix.
                    On exit, the eigenvectors of the repaired tridiagonal matrix.

           LDQ

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

           INDXQ

                     INDXQ is INTEGER array, dimension (N)
                    The permutation which will reintegrate the subproblem just
                    solved back into sorted order, i.e., D( INDXQ( I = 1, N ) )
                    will be in ascending order.

           RHO

                     RHO is DOUBLE PRECISION
                    The subdiagonal element used to create the rank-1
                    modification.

           CUTPNT

                     CUTPNT is INTEGER
                    Contains the location of the last eigenvalue in the leading
                    sub-matrix.  min(1,N) <= CUTPNT <= N.

           QSTORE

                     QSTORE is DOUBLE PRECISION array, dimension (N**2+1)
                    Stores eigenvectors of submatrices encountered during
                    divide and conquer, packed together. QPTR points to
                    beginning of the submatrices.

           QPTR

                     QPTR is INTEGER array, dimension (N+2)
                    List of indices pointing to beginning of submatrices stored
                    in QSTORE. The submatrices are numbered starting at the
                    bottom left of the divide and conquer tree, from left to
                    right and bottom to top.

           PRMPTR

                     PRMPTR is INTEGER array, dimension (N lg N)
                    Contains a list of pointers which indicate where in PERM a
                    level's permutation is stored.  PRMPTR(i+1) - PRMPTR(i)
                    indicates the size of the permutation and also the size of
                    the full, non-deflated problem.

           PERM

                     PERM is INTEGER array, dimension (N lg N)
                    Contains the permutations (from deflation and sorting) to be
                    applied to each eigenblock.

           GIVPTR

                     GIVPTR is INTEGER array, dimension (N lg N)
                    Contains a list of pointers which indicate where in GIVCOL a
                    level's Givens rotations are stored.  GIVPTR(i+1) - GIVPTR(i)
                    indicates the number of Givens rotations.

           GIVCOL

                     GIVCOL is INTEGER array, dimension (2, N lg N)
                    Each pair of numbers indicates a pair of columns to take place
                    in a Givens rotation.

           GIVNUM

                     GIVNUM is DOUBLE PRECISION array, dimension (2, N lg N)
                    Each number indicates the S value to be used in the
                    corresponding Givens rotation.

           WORK

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

           IWORK

                     IWORK is INTEGER array, dimension (4*N)

           INFO

                     INFO is INTEGER
                     = 0:  successful exit.
                     < 0:  if INFO = -i, the i-th argument had an illegal value.
                     > 0:  if INFO = 1, an eigenvalue did not converge

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine dlaed8 (integer ICOMPQ, integer K, integer N, integer QSIZ, double precision,
       dimension( * ) D, double precision, dimension( ldq, * ) Q, integer LDQ, integer,
       dimension( * ) INDXQ, double precision RHO, integer CUTPNT, double precision, dimension( *
       ) Z, double precision, dimension( * ) DLAMDA, double precision, dimension( ldq2, * ) Q2,
       integer LDQ2, double precision, dimension( * ) W, integer, dimension( * ) PERM, integer
       GIVPTR, integer, dimension( 2, * ) GIVCOL, double precision, dimension( 2, * ) GIVNUM,
       integer, dimension( * ) INDXP, integer, dimension( * ) INDX, integer INFO)
       DLAED8 used by DSTEDC. Merges eigenvalues and deflates secular equation. Used when the
       original matrix is dense.

       Purpose:

            DLAED8 merges the two sets of eigenvalues together into a single
            sorted set.  Then it tries to deflate the size of the problem.
            There are two ways in which deflation can occur:  when two or more
            eigenvalues are close together or if there is a tiny element in the
            Z vector.  For each such occurrence the order of the related secular
            equation problem is reduced by one.

       Parameters
           ICOMPQ

                     ICOMPQ is INTEGER
                     = 0:  Compute eigenvalues only.
                     = 1:  Compute eigenvectors of original dense symmetric matrix
                           also.  On entry, Q contains the orthogonal matrix used
                           to reduce the original matrix to tridiagonal form.

           K

                     K is INTEGER
                    The number of non-deflated eigenvalues, and the order of the
                    related secular equation.

           N

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

           QSIZ

                     QSIZ is INTEGER
                    The dimension of the orthogonal matrix used to reduce
                    the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.

           D

                     D is DOUBLE PRECISION array, dimension (N)
                    On entry, the eigenvalues of the two submatrices to be
                    combined.  On exit, the trailing (N-K) updated eigenvalues
                    (those which were deflated) sorted into increasing order.

           Q

                     Q is DOUBLE PRECISION array, dimension (LDQ,N)
                    If ICOMPQ = 0, Q is not referenced.  Otherwise,
                    on entry, Q contains the eigenvectors of the partially solved
                    system which has been previously updated in matrix
                    multiplies with other partially solved eigensystems.
                    On exit, Q contains the trailing (N-K) updated eigenvectors
                    (those which were deflated) in its last N-K columns.

           LDQ

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

           INDXQ

                     INDXQ is INTEGER array, dimension (N)
                    The permutation which separately sorts the two sub-problems
                    in D into ascending order.  Note that elements in the second
                    half of this permutation must first have CUTPNT added to
                    their values in order to be accurate.

           RHO

                     RHO is DOUBLE PRECISION
                    On entry, the off-diagonal element associated with the rank-1
                    cut which originally split the two submatrices which are now
                    being recombined.
                    On exit, RHO has been modified to the value required by
                    DLAED3.

           CUTPNT

                     CUTPNT is INTEGER
                    The location of the last eigenvalue in the leading
                    sub-matrix.  min(1,N) <= CUTPNT <= N.

           Z

                     Z is DOUBLE PRECISION array, dimension (N)
                    On entry, Z contains the updating vector (the last row of
                    the first sub-eigenvector matrix and the first row of the
                    second sub-eigenvector matrix).
                    On exit, the contents of Z are destroyed by the updating
                    process.

           DLAMDA

                     DLAMDA is DOUBLE PRECISION array, dimension (N)
                    A copy of the first K eigenvalues which will be used by
                    DLAED3 to form the secular equation.

           Q2

                     Q2 is DOUBLE PRECISION array, dimension (LDQ2,N)
                    If ICOMPQ = 0, Q2 is not referenced.  Otherwise,
                    a copy of the first K eigenvectors which will be used by
                    DLAED7 in a matrix multiply (DGEMM) to update the new
                    eigenvectors.

           LDQ2

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

           W

                     W is DOUBLE PRECISION array, dimension (N)
                    The first k values of the final deflation-altered z-vector and
                    will be passed to DLAED3.

           PERM

                     PERM is INTEGER array, dimension (N)
                    The permutations (from deflation and sorting) to be applied
                    to each eigenblock.

           GIVPTR

                     GIVPTR is INTEGER
                    The number of Givens rotations which took place in this
                    subproblem.

           GIVCOL

                     GIVCOL is INTEGER array, dimension (2, N)
                    Each pair of numbers indicates a pair of columns to take place
                    in a Givens rotation.

           GIVNUM

                     GIVNUM is DOUBLE PRECISION array, dimension (2, N)
                    Each number indicates the S value to be used in the
                    corresponding Givens rotation.

           INDXP

                     INDXP is INTEGER array, dimension (N)
                    The permutation used to place deflated values of D at the end
                    of the array.  INDXP(1:K) points to the nondeflated D-values
                    and INDXP(K+1:N) points to the deflated eigenvalues.

           INDX

                     INDX is INTEGER array, dimension (N)
                    The permutation used to sort the contents of D into ascending
                    order.

           INFO

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine dlaed9 (integer K, integer KSTART, integer KSTOP, integer N, double precision,
       dimension( * ) D, double precision, dimension( ldq, * ) Q, integer LDQ, double precision
       RHO, double precision, dimension( * ) DLAMDA, double precision, dimension( * ) W, double
       precision, dimension( lds, * ) S, integer LDS, integer INFO)
       DLAED9 used by DSTEDC. Finds the roots of the secular equation and updates the
       eigenvectors. Used when the original matrix is dense.

       Purpose:

            DLAED9 finds the roots of the secular equation, as defined by the
            values in D, Z, and RHO, between KSTART and KSTOP.  It makes the
            appropriate calls to DLAED4 and then stores the new matrix of
            eigenvectors for use in calculating the next level of Z vectors.

       Parameters
           K

                     K is INTEGER
                     The number of terms in the rational function to be solved by
                     DLAED4.  K >= 0.

           KSTART

                     KSTART is INTEGER

           KSTOP

                     KSTOP is INTEGER
                     The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP
                     are to be computed.  1 <= KSTART <= KSTOP <= K.

           N

                     N is INTEGER
                     The number of rows and columns in the Q matrix.
                     N >= K (delation may result in N > K).

           D

                     D is DOUBLE PRECISION array, dimension (N)
                     D(I) contains the updated eigenvalues
                     for KSTART <= I <= KSTOP.

           Q

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

           LDQ

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

           RHO

                     RHO is DOUBLE PRECISION
                     The value of the parameter in the rank one update equation.
                     RHO >= 0 required.

           DLAMDA

                     DLAMDA is DOUBLE PRECISION array, dimension (K)
                     The first K elements of this array contain the old roots
                     of the deflated updating problem.  These are the poles
                     of the secular equation.

           W

                     W is DOUBLE PRECISION array, dimension (K)
                     The first K elements of this array contain the components
                     of the deflation-adjusted updating vector.

           S

                     S is DOUBLE PRECISION array, dimension (LDS, K)
                     Will contain the eigenvectors of the repaired matrix which
                     will be stored for subsequent Z vector calculation and
                     multiplied by the previously accumulated eigenvectors
                     to update the system.

           LDS

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

           INFO

                     INFO is INTEGER
                     = 0:  successful exit.
                     < 0:  if INFO = -i, the i-th argument had an illegal value.
                     > 0:  if INFO = 1, an eigenvalue did not converge

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine dlaeda (integer N, integer TLVLS, integer CURLVL, integer CURPBM, integer,
       dimension( * ) PRMPTR, integer, dimension( * ) PERM, integer, dimension( * ) GIVPTR,
       integer, dimension( 2, * ) GIVCOL, double precision, dimension( 2, * ) GIVNUM, double
       precision, dimension( * ) Q, integer, dimension( * ) QPTR, double precision, dimension( *
       ) Z, double precision, dimension( * ) ZTEMP, integer INFO)
       DLAEDA used by DSTEDC. Computes the Z vector determining the rank-one modification of the
       diagonal matrix. Used when the original matrix is dense.

       Purpose:

            DLAEDA computes the Z vector corresponding to the merge step in the
            CURLVLth step of the merge process with TLVLS steps for the CURPBMth
            problem.

       Parameters
           N

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

           TLVLS

                     TLVLS is INTEGER
                    The total number of merging levels in the overall divide and
                    conquer tree.

           CURLVL

                     CURLVL is INTEGER
                    The current level in the overall merge routine,
                    0 <= curlvl <= tlvls.

           CURPBM

                     CURPBM is INTEGER
                    The current problem in the current level in the overall
                    merge routine (counting from upper left to lower right).

           PRMPTR

                     PRMPTR is INTEGER array, dimension (N lg N)
                    Contains a list of pointers which indicate where in PERM a
                    level's permutation is stored.  PRMPTR(i+1) - PRMPTR(i)
                    indicates the size of the permutation and incidentally the
                    size of the full, non-deflated problem.

           PERM

                     PERM is INTEGER array, dimension (N lg N)
                    Contains the permutations (from deflation and sorting) to be
                    applied to each eigenblock.

           GIVPTR

                     GIVPTR is INTEGER array, dimension (N lg N)
                    Contains a list of pointers which indicate where in GIVCOL a
                    level's Givens rotations are stored.  GIVPTR(i+1) - GIVPTR(i)
                    indicates the number of Givens rotations.

           GIVCOL

                     GIVCOL is INTEGER array, dimension (2, N lg N)
                    Each pair of numbers indicates a pair of columns to take place
                    in a Givens rotation.

           GIVNUM

                     GIVNUM is DOUBLE PRECISION array, dimension (2, N lg N)
                    Each number indicates the S value to be used in the
                    corresponding Givens rotation.

           Q

                     Q is DOUBLE PRECISION array, dimension (N**2)
                    Contains the square eigenblocks from previous levels, the
                    starting positions for blocks are given by QPTR.

           QPTR

                     QPTR is INTEGER array, dimension (N+2)
                    Contains a list of pointers which indicate where in Q an
                    eigenblock is stored.  SQRT( QPTR(i+1) - QPTR(i) ) indicates
                    the size of the block.

           Z

                     Z is DOUBLE PRECISION array, dimension (N)
                    On output this vector contains the updating vector (the last
                    row of the first sub-eigenvector matrix and the first row of
                    the second sub-eigenvector matrix).

           ZTEMP

                     ZTEMP is DOUBLE PRECISION array, dimension (N)

           INFO

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine dlagtf (integer N, double precision, dimension( * ) A, double precision LAMBDA,
       double precision, dimension( * ) B, double precision, dimension( * ) C, double precision
       TOL, double precision, dimension( * ) D, integer, dimension( * ) IN, integer INFO)
       DLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal
       matrix, and λ a scalar, using partial pivoting with row interchanges.

       Purpose:

            DLAGTF factorizes the matrix (T - lambda*I), where T is an n by n
            tridiagonal matrix and lambda is a scalar, as

               T - lambda*I = PLU,

            where P is a permutation matrix, L is a unit lower tridiagonal matrix
            with at most one non-zero sub-diagonal elements per column and U is
            an upper triangular matrix with at most two non-zero super-diagonal
            elements per column.

            The factorization is obtained by Gaussian elimination with partial
            pivoting and implicit row scaling.

            The parameter LAMBDA is included in the routine so that DLAGTF may
            be used, in conjunction with DLAGTS, to obtain eigenvectors of T by
            inverse iteration.

       Parameters
           N

                     N is INTEGER
                     The order of the matrix T.

           A

                     A is DOUBLE PRECISION array, dimension (N)
                     On entry, A must contain the diagonal elements of T.

                     On exit, A is overwritten by the n diagonal elements of the
                     upper triangular matrix U of the factorization of T.

           LAMBDA

                     LAMBDA is DOUBLE PRECISION
                     On entry, the scalar lambda.

           B

                     B is DOUBLE PRECISION array, dimension (N-1)
                     On entry, B must contain the (n-1) super-diagonal elements of
                     T.

                     On exit, B is overwritten by the (n-1) super-diagonal
                     elements of the matrix U of the factorization of T.

           C

                     C is DOUBLE PRECISION array, dimension (N-1)
                     On entry, C must contain the (n-1) sub-diagonal elements of
                     T.

                     On exit, C is overwritten by the (n-1) sub-diagonal elements
                     of the matrix L of the factorization of T.

           TOL

                     TOL is DOUBLE PRECISION
                     On entry, a relative tolerance used to indicate whether or
                     not the matrix (T - lambda*I) is nearly singular. TOL should
                     normally be chose as approximately the largest relative error
                     in the elements of T. For example, if the elements of T are
                     correct to about 4 significant figures, then TOL should be
                     set to about 5*10**(-4). If TOL is supplied as less than eps,
                     where eps is the relative machine precision, then the value
                     eps is used in place of TOL.

           D

                     D is DOUBLE PRECISION array, dimension (N-2)
                     On exit, D is overwritten by the (n-2) second super-diagonal
                     elements of the matrix U of the factorization of T.

           IN

                     IN is INTEGER array, dimension (N)
                     On exit, IN contains details of the permutation matrix P. If
                     an interchange occurred at the kth step of the elimination,
                     then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)
                     returns the smallest positive integer j such that

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

                     where norm( A(j) ) denotes the sum of the absolute values of
                     the jth row of the matrix A. If no such j exists then IN(n)
                     is returned as zero. If IN(n) is returned as positive, then a
                     diagonal element of U is small, indicating that
                     (T - lambda*I) is singular or nearly singular,

           INFO

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine dlamrg (integer N1, integer N2, double precision, dimension( * ) A, integer DTRD1,
       integer DTRD2, integer, dimension( * ) INDEX)
       DLAMRG creates a permutation list to merge the entries of two independently sorted sets
       into a single set sorted in ascending order.

       Purpose:

            DLAMRG will create a permutation list which will merge the elements
            of A (which is composed of two independently sorted sets) into a
            single set which is sorted in ascending order.

       Parameters
           N1

                     N1 is INTEGER

           N2

                     N2 is INTEGER
                    These arguments contain the respective lengths of the two
                    sorted lists to be merged.

           A

                     A is DOUBLE PRECISION array, dimension (N1+N2)
                    The first N1 elements of A contain a list of numbers which
                    are sorted in either ascending or descending order.  Likewise
                    for the final N2 elements.

           DTRD1

                     DTRD1 is INTEGER

           DTRD2

                     DTRD2 is INTEGER
                    These are the strides to be taken through the array A.
                    Allowable strides are 1 and -1.  They indicate whether a
                    subset of A is sorted in ascending (DTRDx = 1) or descending
                    (DTRDx = -1) order.

           INDEX

                     INDEX is INTEGER array, dimension (N1+N2)
                    On exit this array will contain a permutation such that
                    if B( I ) = A( INDEX( I ) ) for I=1,N1+N2, then B will be
                    sorted in ascending order.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine dlartgs (double precision X, double precision Y, double precision SIGMA, double
       precision CS, double precision SN)
       DLARTGS generates a plane rotation designed to introduce a bulge in implicit QR iteration
       for the bidiagonal SVD problem.

       Purpose:

            DLARTGS generates a plane rotation designed to introduce a bulge in
            Golub-Reinsch-style implicit QR iteration for the bidiagonal SVD
            problem. X and Y are the top-row entries, and SIGMA is the shift.
            The computed CS and SN define a plane rotation satisfying

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

            with R nonnegative.  If X^2 - SIGMA and X * Y are 0, then the
            rotation is by PI/2.

       Parameters
           X

                     X is DOUBLE PRECISION
                     The (1,1) entry of an upper bidiagonal matrix.

           Y

                     Y is DOUBLE PRECISION
                     The (1,2) entry of an upper bidiagonal matrix.

           SIGMA

                     SIGMA is DOUBLE PRECISION
                     The shift.

           CS

                     CS is DOUBLE PRECISION
                     The cosine of the rotation.

           SN

                     SN is DOUBLE PRECISION
                     The sine of the rotation.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine dlasq1 (integer N, double precision, dimension( * ) D, double precision, dimension(
       * ) E, double precision, dimension( * ) WORK, integer INFO)
       DLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.

       Purpose:

            DLASQ1 computes the singular values of a real N-by-N bidiagonal
            matrix with diagonal D and off-diagonal E. The singular values
            are computed to high relative accuracy, in the absence of
            denormalization, underflow and overflow. The algorithm was first
            presented in

            'Accurate singular values and differential qd algorithms' by K. V.
            Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
            1994,

            and the present implementation is described in 'An implementation of
            the dqds Algorithm (Positive Case)', LAPACK Working Note.

       Parameters
           N

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

           D

                     D is DOUBLE PRECISION array, dimension (N)
                   On entry, D contains the diagonal elements of the
                   bidiagonal matrix whose SVD is desired. On normal exit,
                   D contains the singular values in decreasing order.

           E

                     E is DOUBLE PRECISION array, dimension (N)
                   On entry, elements E(1:N-1) contain the off-diagonal elements
                   of the bidiagonal matrix whose SVD is desired.
                   On exit, E is overwritten.

           WORK

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

           INFO

                     INFO is INTEGER
                   = 0: successful exit
                   < 0: if INFO = -i, the i-th argument had an illegal value
                   > 0: the algorithm failed
                        = 1, a split was marked by a positive value in E
                        = 2, current block of Z not diagonalized after 100*N
                             iterations (in inner while loop)  On exit D and E
                             represent a matrix with the same singular values
                             which the calling subroutine could use to finish the
                             computation, or even feed back into DLASQ1
                        = 3, termination criterion of outer while loop not met
                             (program created more than N unreduced blocks)

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine dlasq2 (integer N, double precision, dimension( * ) Z, integer INFO)
       DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix
       associated with the qd Array Z to high relative accuracy. Used by sbdsqr and sstegr.

       Purpose:

            DLASQ2 computes all the eigenvalues of the symmetric positive
            definite tridiagonal matrix associated with the qd array Z to high
            relative accuracy are computed to high relative accuracy, in the
            absence of denormalization, underflow and overflow.

            To see the relation of Z to the tridiagonal matrix, let L be a
            unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
            let U be an upper bidiagonal matrix with 1's above and diagonal
            Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
            symmetric tridiagonal to which it is similar.

            Note : DLASQ2 defines a logical variable, IEEE, which is true
            on machines which follow ieee-754 floating-point standard in their
            handling of infinities and NaNs, and false otherwise. This variable
            is passed to DLASQ3.

       Parameters
           N

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

           Z

                     Z is DOUBLE PRECISION array, dimension ( 4*N )
                   On entry Z holds the qd array. On exit, entries 1 to N hold
                   the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
                   trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
                   N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
                   holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
                   shifts that failed.

           INFO

                     INFO is INTEGER
                   = 0: successful exit
                   < 0: if the i-th argument is a scalar and had an illegal
                        value, then INFO = -i, if the i-th argument is an
                        array and the j-entry had an illegal value, then
                        INFO = -(i*100+j)
                   > 0: the algorithm failed
                         = 1, a split was marked by a positive value in E
                         = 2, current block of Z not diagonalized after 100*N
                              iterations (in inner while loop).  On exit Z holds
                              a qd array with the same eigenvalues as the given Z.
                         = 3, termination criterion of outer while loop not met
                              (program created more than N unreduced blocks)

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Further Details:

             Local Variables: I0:N0 defines a current unreduced segment of Z.
             The shifts are accumulated in SIGMA. Iteration count is in ITER.
             Ping-pong is controlled by PP (alternates between 0 and 1).

   subroutine dlasq3 (integer I0, integer N0, double precision, dimension( * ) Z, integer PP,
       double precision DMIN, double precision SIGMA, double precision DESIG, double precision
       QMAX, integer NFAIL, integer ITER, integer NDIV, logical IEEE, integer TTYPE, double
       precision DMIN1, double precision DMIN2, double precision DN, double precision DN1, double
       precision DN2, double precision G, double precision TAU)
       DLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr.

       Purpose:

            DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.
            In case of failure it changes shifts, and tries again until output
            is positive.

       Parameters
           I0

                     I0 is INTEGER
                    First index.

           N0

                     N0 is INTEGER
                    Last index.

           Z

                     Z is DOUBLE PRECISION array, dimension ( 4*N0 )
                    Z holds the qd array.

           PP

                     PP is INTEGER
                    PP=0 for ping, PP=1 for pong.
                    PP=2 indicates that flipping was applied to the Z array
                    and that the initial tests for deflation should not be
                    performed.

           DMIN

                     DMIN is DOUBLE PRECISION
                    Minimum value of d.

           SIGMA

                     SIGMA is DOUBLE PRECISION
                    Sum of shifts used in current segment.

           DESIG

                     DESIG is DOUBLE PRECISION
                    Lower order part of SIGMA

           QMAX

                     QMAX is DOUBLE PRECISION
                    Maximum value of q.

           NFAIL

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

           ITER

                     ITER is INTEGER
                    Increment ITER by 1 for each iteration.

           NDIV

                     NDIV is INTEGER
                    Increment NDIV by 1 for each division.

           IEEE

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

           TTYPE

                     TTYPE is INTEGER
                    Shift type.

           DMIN1

                     DMIN1 is DOUBLE PRECISION

           DMIN2

                     DMIN2 is DOUBLE PRECISION

           DN

                     DN is DOUBLE PRECISION

           DN1

                     DN1 is DOUBLE PRECISION

           DN2

                     DN2 is DOUBLE PRECISION

           G

                     G is DOUBLE PRECISION

           TAU

                     TAU is DOUBLE PRECISION

                    These are passed as arguments in order to save their values
                    between calls to DLASQ3.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine dlasq4 (integer I0, integer N0, double precision, dimension( * ) Z, integer PP,
       integer N0IN, double precision DMIN, double precision DMIN1, double precision DMIN2,
       double precision DN, double precision DN1, double precision DN2, double precision TAU,
       integer TTYPE, double precision G)
       DLASQ4 computes an approximation to the smallest eigenvalue using values of d from the
       previous transform. Used by sbdsqr.

       Purpose:

            DLASQ4 computes an approximation TAU to the smallest eigenvalue
            using values of d from the previous transform.

       Parameters
           I0

                     I0 is INTEGER
                   First index.

           N0

                     N0 is INTEGER
                   Last index.

           Z

                     Z is DOUBLE PRECISION array, dimension ( 4*N0 )
                   Z holds the qd array.

           PP

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

           N0IN

                     N0IN is INTEGER
                   The value of N0 at start of EIGTEST.

           DMIN

                     DMIN is DOUBLE PRECISION
                   Minimum value of d.

           DMIN1

                     DMIN1 is DOUBLE PRECISION
                   Minimum value of d, excluding D( N0 ).

           DMIN2

                     DMIN2 is DOUBLE PRECISION
                   Minimum value of d, excluding D( N0 ) and D( N0-1 ).

           DN

                     DN is DOUBLE PRECISION
                   d(N)

           DN1

                     DN1 is DOUBLE PRECISION
                   d(N-1)

           DN2

                     DN2 is DOUBLE PRECISION
                   d(N-2)

           TAU

                     TAU is DOUBLE PRECISION
                   This is the shift.

           TTYPE

                     TTYPE is INTEGER
                   Shift type.

           G

                     G is DOUBLE PRECISION
                   G is passed as an argument in order to save its value between
                   calls to DLASQ4.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Further Details:

             CNST1 = 9/16

   subroutine dlasq5 (integer I0, integer N0, double precision, dimension( * ) Z, integer PP,
       double precision TAU, double precision SIGMA, double precision DMIN, double precision
       DMIN1, double precision DMIN2, double precision DN, double precision DNM1, double
       precision DNM2, logical IEEE, double precision EPS)
       DLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.

       Purpose:

            DLASQ5 computes one dqds transform in ping-pong form, one
            version for IEEE machines another for non IEEE machines.

       Parameters
           I0

                     I0 is INTEGER
                   First index.

           N0

                     N0 is INTEGER
                   Last index.

           Z

                     Z is DOUBLE PRECISION array, dimension ( 4*N )
                   Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
                   an extra argument.

           PP

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

           TAU

                     TAU is DOUBLE PRECISION
                   This is the shift.

           SIGMA

                     SIGMA is DOUBLE PRECISION
                   This is the accumulated shift up to this step.

           DMIN

                     DMIN is DOUBLE PRECISION
                   Minimum value of d.

           DMIN1

                     DMIN1 is DOUBLE PRECISION
                   Minimum value of d, excluding D( N0 ).

           DMIN2

                     DMIN2 is DOUBLE PRECISION
                   Minimum value of d, excluding D( N0 ) and D( N0-1 ).

           DN

                     DN is DOUBLE PRECISION
                   d(N0), the last value of d.

           DNM1

                     DNM1 is DOUBLE PRECISION
                   d(N0-1).

           DNM2

                     DNM2 is DOUBLE PRECISION
                   d(N0-2).

           IEEE

                     IEEE is LOGICAL
                   Flag for IEEE or non IEEE arithmetic.

           EPS

                     EPS is DOUBLE PRECISION
                   This is the value of epsilon used.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine dlasq6 (integer I0, integer N0, double precision, dimension( * ) Z, integer PP,
       double precision DMIN, double precision DMIN1, double precision DMIN2, double precision
       DN, double precision DNM1, double precision DNM2)
       DLASQ6 computes one dqd transform in ping-pong form. Used by sbdsqr and sstegr.

       Purpose:

            DLASQ6 computes one dqd (shift equal to zero) transform in
            ping-pong form, with protection against underflow and overflow.

       Parameters
           I0

                     I0 is INTEGER
                   First index.

           N0

                     N0 is INTEGER
                   Last index.

           Z

                     Z is DOUBLE PRECISION array, dimension ( 4*N )
                   Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
                   an extra argument.

           PP

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

           DMIN

                     DMIN is DOUBLE PRECISION
                   Minimum value of d.

           DMIN1

                     DMIN1 is DOUBLE PRECISION
                   Minimum value of d, excluding D( N0 ).

           DMIN2

                     DMIN2 is DOUBLE PRECISION
                   Minimum value of d, excluding D( N0 ) and D( N0-1 ).

           DN

                     DN is DOUBLE PRECISION
                   d(N0), the last value of d.

           DNM1

                     DNM1 is DOUBLE PRECISION
                   d(N0-1).

           DNM2

                     DNM2 is DOUBLE PRECISION
                   d(N0-2).

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine dlasrt (character ID, integer N, double precision, dimension( * ) D, integer INFO)
       DLASRT sorts numbers in increasing or decreasing order.

       Purpose:

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

            Use Quick Sort, reverting to Insertion sort on arrays of
            size <= 20. Dimension of STACK limits N to about 2**32.

       Parameters
           ID

                     ID is CHARACTER*1
                     = 'I': sort D in increasing order;
                     = 'D': sort D in decreasing order.

           N

                     N is INTEGER
                     The length of the array D.

           D

                     D is DOUBLE PRECISION array, dimension (N)
                     On entry, the array to be sorted.
                     On exit, D has been sorted into increasing order
                     (D(1) <= ... <= D(N) ) or into decreasing order
                     (D(1) >= ... >= D(N) ), depending on ID.

           INFO

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

       Purpose:

            DSTEBZ computes the eigenvalues of a symmetric tridiagonal
            matrix T.  The user may ask for all eigenvalues, all eigenvalues
            in the half-open interval (VL, VU], or the IL-th through IU-th
            eigenvalues.

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

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

       Parameters
           RANGE

                     RANGE is CHARACTER*1
                     = 'A': ('All')   all eigenvalues will be found.
                     = 'V': ('Value') all eigenvalues in the half-open interval
                                      (VL, VU] will be found.
                     = 'I': ('Index') the IL-th through IU-th eigenvalues (of the
                                      entire matrix) will be found.

           ORDER

                     ORDER is CHARACTER*1
                     = 'B': ('By Block') the eigenvalues will be grouped by
                                         split-off block (see IBLOCK, ISPLIT) and
                                         ordered from smallest to largest within
                                         the block.
                     = 'E': ('Entire matrix')
                                         the eigenvalues for the entire matrix
                                         will be ordered from smallest to
                                         largest.

           N

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

           VL

                     VL is DOUBLE PRECISION

                     If RANGE='V', the lower bound of the interval to
                     be searched for eigenvalues.  Eigenvalues less than or equal
                     to VL, or greater than VU, will not be returned.  VL < VU.
                     Not referenced if RANGE = 'A' or 'I'.

           VU

                     VU is DOUBLE PRECISION

                     If RANGE='V', the upper bound of the interval to
                     be searched for eigenvalues.  Eigenvalues less than or equal
                     to VL, or greater than VU, will not be returned.  VL < VU.
                     Not referenced if RANGE = 'A' or 'I'.

           IL

                     IL is INTEGER

                     If RANGE='I', the index of the
                     smallest eigenvalue to be returned.
                     1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
                     Not referenced if RANGE = 'A' or 'V'.

           IU

                     IU is INTEGER

                     If RANGE='I', the index of the
                     largest eigenvalue to be returned.
                     1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
                     Not referenced if RANGE = 'A' or 'V'.

           ABSTOL

                     ABSTOL is DOUBLE PRECISION
                     The absolute tolerance for the eigenvalues.  An eigenvalue
                     (or cluster) is considered to be located if it has been
                     determined to lie in an interval whose width is ABSTOL or
                     less.  If ABSTOL is less than or equal to zero, then ULP*|T|
                     will be used, where |T| means the 1-norm of T.

                     Eigenvalues will be computed most accurately when ABSTOL is
                     set to twice the underflow threshold 2*DLAMCH('S'), not zero.

           D

                     D is DOUBLE PRECISION array, dimension (N)
                     The n diagonal elements of the tridiagonal matrix T.

           E

                     E is DOUBLE PRECISION array, dimension (N-1)
                     The (n-1) off-diagonal elements of the tridiagonal matrix T.

           M

                     M is INTEGER
                     The actual number of eigenvalues found. 0 <= M <= N.
                     (See also the description of INFO=2,3.)

           NSPLIT

                     NSPLIT is INTEGER
                     The number of diagonal blocks in the matrix T.
                     1 <= NSPLIT <= N.

           W

                     W is DOUBLE PRECISION array, dimension (N)
                     On exit, the first M elements of W will contain the
                     eigenvalues.  (DSTEBZ may use the remaining N-M elements as
                     workspace.)

           IBLOCK

                     IBLOCK is INTEGER array, dimension (N)
                     At each row/column j where E(j) is zero or small, the
                     matrix T is considered to split into a block diagonal
                     matrix.  On exit, if INFO = 0, IBLOCK(i) specifies to which
                     block (from 1 to the number of blocks) the eigenvalue W(i)
                     belongs.  (DSTEBZ may use the remaining N-M elements as
                     workspace.)

           ISPLIT

                     ISPLIT is INTEGER array, dimension (N)
                     The splitting points, at which T breaks up into submatrices.
                     The first submatrix consists of rows/columns 1 to ISPLIT(1),
                     the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
                     etc., and the NSPLIT-th consists of rows/columns
                     ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
                     (Only the first NSPLIT elements will actually be used, but
                     since the user cannot know a priori what value NSPLIT will
                     have, N words must be reserved for ISPLIT.)

           WORK

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

           IWORK

                     IWORK is INTEGER array, dimension (3*N)

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  if INFO = -i, the i-th argument had an illegal value
                     > 0:  some or all of the eigenvalues failed to converge or
                           were not computed:
                           =1 or 3: Bisection failed to converge for some
                                   eigenvalues; these eigenvalues are flagged by a
                                   negative block number.  The effect is that the
                                   eigenvalues may not be as accurate as the
                                   absolute and relative tolerances.  This is
                                   generally caused by unexpectedly inaccurate
                                   arithmetic.
                           =2 or 3: RANGE='I' only: Not all of the eigenvalues
                                   IL:IU were found.
                                   Effect: M < IU+1-IL
                                   Cause:  non-monotonic arithmetic, causing the
                                           Sturm sequence to be non-monotonic.
                                   Cure:   recalculate, using RANGE='A', and pick
                                           out eigenvalues IL:IU.  In some cases,
                                           increasing the PARAMETER 'FUDGE' may
                                           make things work.
                           = 4:    RANGE='I', and the Gershgorin interval
                                   initially used was too small.  No eigenvalues
                                   were computed.
                                   Probable cause: your machine has sloppy
                                                   floating-point arithmetic.
                                   Cure: Increase the PARAMETER 'FUDGE',
                                         recompile, and try again.

       Internal Parameters:

             RELFAC  DOUBLE PRECISION, default = 2.0e0
                     The relative tolerance.  An interval (a,b] lies within
                     'relative tolerance' if  b-a < RELFAC*ulp*max(|a|,|b|),
                     where 'ulp' is the machine precision (distance from 1 to
                     the next larger floating point number.)

             FUDGE   DOUBLE PRECISION, default = 2
                     A 'fudge factor' to widen the Gershgorin intervals.  Ideally,
                     a value of 1 should work, but on machines with sloppy
                     arithmetic, this needs to be larger.  The default for
                     publicly released versions should be large enough to handle
                     the worst machine around.  Note that this has no effect
                     on accuracy of the solution.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

       Purpose:

            DSTEDC computes all eigenvalues and, optionally, eigenvectors of a
            symmetric tridiagonal matrix using the divide and conquer method.
            The eigenvectors of a full or band real symmetric matrix can also be
            found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this
            matrix to tridiagonal form.

            This code makes very mild assumptions about floating point
            arithmetic. It will work on machines with a guard digit in
            add/subtract, or on those binary machines without guard digits
            which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
            It could conceivably fail on hexadecimal or decimal machines
            without guard digits, but we know of none.  See DLAED3 for details.

       Parameters
           COMPZ

                     COMPZ is CHARACTER*1
                     = 'N':  Compute eigenvalues only.
                     = 'I':  Compute eigenvectors of tridiagonal matrix also.
                     = 'V':  Compute eigenvectors of original dense symmetric
                             matrix also.  On entry, Z contains the orthogonal
                             matrix used to reduce the original matrix to
                             tridiagonal form.

           N

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

           D

                     D is DOUBLE PRECISION array, dimension (N)
                     On entry, the diagonal elements of the tridiagonal matrix.
                     On exit, if INFO = 0, the eigenvalues in ascending order.

           E

                     E is DOUBLE PRECISION array, dimension (N-1)
                     On entry, the subdiagonal elements of the tridiagonal matrix.
                     On exit, E has been destroyed.

           Z

                     Z is DOUBLE PRECISION array, dimension (LDZ,N)
                     On entry, if COMPZ = 'V', then Z contains the orthogonal
                     matrix used in the reduction to tridiagonal form.
                     On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
                     orthonormal eigenvectors of the original symmetric matrix,
                     and if COMPZ = 'I', Z contains the orthonormal eigenvectors
                     of the symmetric tridiagonal matrix.
                     If  COMPZ = 'N', then Z is not referenced.

           LDZ

                     LDZ is INTEGER
                     The leading dimension of the array Z.  LDZ >= 1.
                     If eigenvectors are desired, then LDZ >= max(1,N).

           WORK

                     WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
                     On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

           LWORK

                     LWORK is INTEGER
                     The dimension of the array WORK.
                     If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
                     If COMPZ = 'V' and N > 1 then LWORK must be at least
                                    ( 1 + 3*N + 2*N*lg N + 4*N**2 ),
                                    where lg( N ) = smallest integer k such
                                    that 2**k >= N.
                     If COMPZ = 'I' and N > 1 then LWORK must be at least
                                    ( 1 + 4*N + N**2 ).
                     Note that for COMPZ = 'I' or 'V', then if N is less than or
                     equal to the minimum divide size, usually 25, then LWORK need
                     only be max(1,2*(N-1)).

                     If LWORK = -1, then a workspace query is assumed; the routine
                     only calculates the optimal size of the WORK array, returns
                     this value as the first entry of the WORK array, and no error
                     message related to LWORK is issued by XERBLA.

           IWORK

                     IWORK is INTEGER array, dimension (MAX(1,LIWORK))
                     On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.

           LIWORK

                     LIWORK is INTEGER
                     The dimension of the array IWORK.
                     If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
                     If COMPZ = 'V' and N > 1 then LIWORK must be at least
                                    ( 6 + 6*N + 5*N*lg N ).
                     If COMPZ = 'I' and N > 1 then LIWORK must be at least
                                    ( 3 + 5*N ).
                     Note that for COMPZ = 'I' or 'V', then if N is less than or
                     equal to the minimum divide size, usually 25, then LIWORK
                     need only be 1.

                     If LIWORK = -1, then a workspace query is assumed; the
                     routine only calculates the optimal size of the IWORK array,
                     returns this value as the first entry of the IWORK array, and
                     no error message related to LIWORK is issued by XERBLA.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit.
                     < 0:  if INFO = -i, the i-th argument had an illegal value.
                     > 0:  The algorithm failed to compute an eigenvalue while
                           working on the submatrix lying in rows and columns
                           INFO/(N+1) through mod(INFO,N+1).

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

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

       Purpose:

            DSTEQR computes all eigenvalues and, optionally, eigenvectors of a
            symmetric tridiagonal matrix using the implicit QL or QR method.
            The eigenvectors of a full or band symmetric matrix can also be found
            if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to
            tridiagonal form.

       Parameters
           COMPZ

                     COMPZ is CHARACTER*1
                     = 'N':  Compute eigenvalues only.
                     = 'V':  Compute eigenvalues and eigenvectors of the original
                             symmetric matrix.  On entry, Z must contain the
                             orthogonal matrix used to reduce the original matrix
                             to tridiagonal form.
                     = 'I':  Compute eigenvalues and eigenvectors of the
                             tridiagonal matrix.  Z is initialized to the identity
                             matrix.

           N

                     N is INTEGER
                     The order of the matrix.  N >= 0.

           D

                     D is DOUBLE PRECISION array, dimension (N)
                     On entry, the diagonal elements of the tridiagonal matrix.
                     On exit, if INFO = 0, the eigenvalues in ascending order.

           E

                     E is DOUBLE PRECISION array, dimension (N-1)
                     On entry, the (n-1) subdiagonal elements of the tridiagonal
                     matrix.
                     On exit, E has been destroyed.

           Z

                     Z is DOUBLE PRECISION array, dimension (LDZ, N)
                     On entry, if  COMPZ = 'V', then Z contains the orthogonal
                     matrix used in the reduction to tridiagonal form.
                     On exit, if INFO = 0, then if  COMPZ = 'V', Z contains the
                     orthonormal eigenvectors of the original symmetric matrix,
                     and if COMPZ = 'I', Z contains the orthonormal eigenvectors
                     of the symmetric tridiagonal matrix.
                     If COMPZ = 'N', then Z is not referenced.

           LDZ

                     LDZ is INTEGER
                     The leading dimension of the array Z.  LDZ >= 1, and if
                     eigenvectors are desired, then  LDZ >= max(1,N).

           WORK

                     WORK is DOUBLE PRECISION array, dimension (max(1,2*N-2))
                     If COMPZ = 'N', then WORK is not referenced.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  if INFO = -i, the i-th argument had an illegal value
                     > 0:  the algorithm has failed to find all the eigenvalues in
                           a total of 30*N iterations; if INFO = i, then i
                           elements of E have not converged to zero; on exit, D
                           and E contain the elements of a symmetric tridiagonal
                           matrix which is orthogonally similar to the original
                           matrix.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

       Purpose:

            DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
            using the Pal-Walker-Kahan variant of the QL or QR algorithm.

       Parameters
           N

                     N is INTEGER
                     The order of the matrix.  N >= 0.

           D

                     D is DOUBLE PRECISION array, dimension (N)
                     On entry, the n diagonal elements of the tridiagonal matrix.
                     On exit, if INFO = 0, the eigenvalues in ascending order.

           E

                     E is DOUBLE PRECISION array, dimension (N-1)
                     On entry, the (n-1) subdiagonal elements of the tridiagonal
                     matrix.
                     On exit, E has been destroyed.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  if INFO = -i, the i-th argument had an illegal value
                     > 0:  the algorithm failed to find all of the eigenvalues in
                           a total of 30*N iterations; if INFO = i, then i
                           elements of E have not converged to zero.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   integer function iladiag (character DIAG)
       ILADIAG

       Purpose:

            This subroutine translated from a character string specifying if a
            matrix has unit diagonal or not to the relevant BLAST-specified
            integer constant.

            ILADIAG returns an INTEGER.  If ILADIAG < 0, then the input is not a
            character indicating a unit or non-unit diagonal.  Otherwise ILADIAG
            returns the constant value corresponding to DIAG.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   integer function ilaprec (character PREC)
       ILAPREC

       Purpose:

            This subroutine translated from a character string specifying an
            intermediate precision to the relevant BLAST-specified integer
            constant.

            ILAPREC returns an INTEGER.  If ILAPREC < 0, then the input is not a
            character indicating a supported intermediate precision.  Otherwise
            ILAPREC returns the constant value corresponding to PREC.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   integer function ilatrans (character TRANS)
       ILATRANS

       Purpose:

            This subroutine translates from a character string specifying a
            transposition operation to the relevant BLAST-specified integer
            constant.

            ILATRANS returns an INTEGER.  If ILATRANS < 0, then the input is not
            a character indicating a transposition operator.  Otherwise ILATRANS
            returns the constant value corresponding to TRANS.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   integer function ilauplo (character UPLO)
       ILAUPLO

       Purpose:

            This subroutine translated from a character string specifying a
            upper- or lower-triangular matrix to the relevant BLAST-specified
            integer constant.

            ILAUPLO returns an INTEGER.  If ILAUPLO < 0, then the input is not
            a character indicating an upper- or lower-triangular matrix.
            Otherwise ILAUPLO returns the constant value corresponding to UPLO.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

       Purpose:

            SBDSDC computes the singular value decomposition (SVD) of a real
            N-by-N (upper or lower) bidiagonal matrix B:  B = U * S * VT,
            using a divide and conquer method, where S is a diagonal matrix
            with non-negative diagonal elements (the singular values of B), and
            U and VT are orthogonal matrices of left and right singular vectors,
            respectively. SBDSDC can be used to compute all singular values,
            and optionally, singular vectors or singular vectors in compact form.

            This code makes very mild assumptions about floating point
            arithmetic. It will work on machines with a guard digit in
            add/subtract, or on those binary machines without guard digits
            which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
            It could conceivably fail on hexadecimal or decimal machines
            without guard digits, but we know of none.  See SLASD3 for details.

            The code currently calls SLASDQ if singular values only are desired.
            However, it can be slightly modified to compute singular values
            using the divide and conquer method.

       Parameters
           UPLO

                     UPLO is CHARACTER*1
                     = 'U':  B is upper bidiagonal.
                     = 'L':  B is lower bidiagonal.

           COMPQ

                     COMPQ is CHARACTER*1
                     Specifies whether singular vectors are to be computed
                     as follows:
                     = 'N':  Compute singular values only;
                     = 'P':  Compute singular values and compute singular
                             vectors in compact form;
                     = 'I':  Compute singular values and singular vectors.

           N

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

           D

                     D is REAL array, dimension (N)
                     On entry, the n diagonal elements of the bidiagonal matrix B.
                     On exit, if INFO=0, the singular values of B.

           E

                     E is REAL array, dimension (N-1)
                     On entry, the elements of E contain the offdiagonal
                     elements of the bidiagonal matrix whose SVD is desired.
                     On exit, E has been destroyed.

           U

                     U is REAL array, dimension (LDU,N)
                     If  COMPQ = 'I', then:
                        On exit, if INFO = 0, U contains the left singular vectors
                        of the bidiagonal matrix.
                     For other values of COMPQ, U is not referenced.

           LDU

                     LDU is INTEGER
                     The leading dimension of the array U.  LDU >= 1.
                     If singular vectors are desired, then LDU >= max( 1, N ).

           VT

                     VT is REAL array, dimension (LDVT,N)
                     If  COMPQ = 'I', then:
                        On exit, if INFO = 0, VT**T contains the right singular
                        vectors of the bidiagonal matrix.
                     For other values of COMPQ, VT is not referenced.

           LDVT

                     LDVT is INTEGER
                     The leading dimension of the array VT.  LDVT >= 1.
                     If singular vectors are desired, then LDVT >= max( 1, N ).

           Q

                     Q is REAL array, dimension (LDQ)
                     If  COMPQ = 'P', then:
                        On exit, if INFO = 0, Q and IQ contain the left
                        and right singular vectors in a compact form,
                        requiring O(N log N) space instead of 2*N**2.
                        In particular, Q contains all the REAL data in
                        LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
                        words of memory, where SMLSIZ is returned by ILAENV and
                        is equal to the maximum size of the subproblems at the
                        bottom of the computation tree (usually about 25).
                     For other values of COMPQ, Q is not referenced.

           IQ

                     IQ is INTEGER array, dimension (LDIQ)
                     If  COMPQ = 'P', then:
                        On exit, if INFO = 0, Q and IQ contain the left
                        and right singular vectors in a compact form,
                        requiring O(N log N) space instead of 2*N**2.
                        In particular, IQ contains all INTEGER data in
                        LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
                        words of memory, where SMLSIZ is returned by ILAENV and
                        is equal to the maximum size of the subproblems at the
                        bottom of the computation tree (usually about 25).
                     For other values of COMPQ, IQ is not referenced.

           WORK

                     WORK is REAL array, dimension (MAX(1,LWORK))
                     If COMPQ = 'N' then LWORK >= (4 * N).
                     If COMPQ = 'P' then LWORK >= (6 * N).
                     If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).

           IWORK

                     IWORK is INTEGER array, dimension (8*N)

           INFO

                     INFO is INTEGER
                     = 0:  successful exit.
                     < 0:  if INFO = -i, the i-th argument had an illegal value.
                     > 0:  The algorithm failed to compute a singular value.
                           The update process of divide and conquer failed.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

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

       Purpose:

            SBDSQR computes the singular values and, optionally, the right and/or
            left singular vectors from the singular value decomposition (SVD) of
            a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
            zero-shift QR algorithm.  The SVD of B has the form

               B = Q * S * P**T

            where S is the diagonal matrix of singular values, Q is an orthogonal
            matrix of left singular vectors, and P is an orthogonal matrix of
            right singular vectors.  If left singular vectors are requested, this
            subroutine actually returns U*Q instead of Q, and, if right singular
            vectors are requested, this subroutine returns P**T*VT instead of
            P**T, for given real input matrices U and VT.  When U and VT are the
            orthogonal matrices that reduce a general matrix A to bidiagonal
            form:  A = U*B*VT, as computed by SGEBRD, then

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

            is the SVD of A.  Optionally, the subroutine may also compute Q**T*C
            for a given real input matrix C.

            See 'Computing  Small Singular Values of Bidiagonal Matrices With
            Guaranteed High Relative Accuracy,' by J. Demmel and W. Kahan,
            LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
            no. 5, pp. 873-912, Sept 1990) and
            'Accurate singular values and differential qd algorithms,' by
            B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
            Department, University of California at Berkeley, July 1992
            for a detailed description of the algorithm.

       Parameters
           UPLO

                     UPLO is CHARACTER*1
                     = 'U':  B is upper bidiagonal;
                     = 'L':  B is lower bidiagonal.

           N

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

           NCVT

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

           NRU

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

           NCC

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

           D

                     D is REAL array, dimension (N)
                     On entry, the n diagonal elements of the bidiagonal matrix B.
                     On exit, if INFO=0, the singular values of B in decreasing
                     order.

           E

                     E is REAL array, dimension (N-1)
                     On entry, the N-1 offdiagonal elements of the bidiagonal
                     matrix B.
                     On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
                     will contain the diagonal and superdiagonal elements of a
                     bidiagonal matrix orthogonally equivalent to the one given
                     as input.

           VT

                     VT is REAL array, dimension (LDVT, NCVT)
                     On entry, an N-by-NCVT matrix VT.
                     On exit, VT is overwritten by P**T * VT.
                     Not referenced if NCVT = 0.

           LDVT

                     LDVT is INTEGER
                     The leading dimension of the array VT.
                     LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.

           U

                     U is REAL array, dimension (LDU, N)
                     On entry, an NRU-by-N matrix U.
                     On exit, U is overwritten by U * Q.
                     Not referenced if NRU = 0.

           LDU

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

           C

                     C is REAL array, dimension (LDC, NCC)
                     On entry, an N-by-NCC matrix C.
                     On exit, C is overwritten by Q**T * C.
                     Not referenced if NCC = 0.

           LDC

                     LDC is INTEGER
                     The leading dimension of the array C.
                     LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.

           WORK

                     WORK is REAL array, dimension (4*N)

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  If INFO = -i, the i-th argument had an illegal value
                     > 0:
                        if NCVT = NRU = NCC = 0,
                           = 1, a split was marked by a positive value in E
                           = 2, current block of Z not diagonalized after 30*N
                                iterations (in inner while loop)
                           = 3, termination criterion of outer while loop not met
                                (program created more than N unreduced blocks)
                        else NCVT = NRU = NCC = 0,
                              the algorithm did not converge; D and E contain the
                              elements of a bidiagonal matrix which is orthogonally
                              similar to the input matrix B;  if INFO = i, i
                              elements of E have not converged to zero.

       Internal Parameters:

             TOLMUL  REAL, default = max(10,min(100,EPS**(-1/8)))
                     TOLMUL controls the convergence criterion of the QR loop.
                     If it is positive, TOLMUL*EPS is the desired relative
                        precision in the computed singular values.
                     If it is negative, abs(TOLMUL*EPS*sigma_max) is the
                        desired absolute accuracy in the computed singular
                        values (corresponds to relative accuracy
                        abs(TOLMUL*EPS) in the largest singular value.
                     abs(TOLMUL) should be between 1 and 1/EPS, and preferably
                        between 10 (for fast convergence) and .1/EPS
                        (for there to be some accuracy in the results).
                     Default is to lose at either one eighth or 2 of the
                        available decimal digits in each computed singular value
                        (whichever is smaller).

             MAXITR  INTEGER, default = 6
                     MAXITR controls the maximum number of passes of the
                     algorithm through its inner loop. The algorithms stops
                     (and so fails to converge) if the number of passes
                     through the inner loop exceeds MAXITR*N**2.

       Note:

             Bug report from Cezary Dendek.
             On March 23rd 2017, the INTEGER variable MAXIT = MAXITR*N**2 is
             removed since it can overflow pretty easily (for N larger or equal
             than 18,919). We instead use MAXITDIVN = MAXITR*N.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

       Purpose:

            SDISNA computes the reciprocal condition numbers for the eigenvectors
            of a real symmetric or complex Hermitian matrix or for the left or
            right singular vectors of a general m-by-n matrix. The reciprocal
            condition number is the 'gap' between the corresponding eigenvalue or
            singular value and the nearest other one.

            The bound on the error, measured by angle in radians, in the I-th
            computed vector is given by

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

            where ANORM = 2-norm(A) = max( abs( D(j) ) ).  SEP(I) is not allowed
            to be smaller than SLAMCH( 'E' )*ANORM in order to limit the size of
            the error bound.

            SDISNA may also be used to compute error bounds for eigenvectors of
            the generalized symmetric definite eigenproblem.

       Parameters
           JOB

                     JOB is CHARACTER*1
                     Specifies for which problem the reciprocal condition numbers
                     should be computed:
                     = 'E':  the eigenvectors of a symmetric/Hermitian matrix;
                     = 'L':  the left singular vectors of a general matrix;
                     = 'R':  the right singular vectors of a general matrix.

           M

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

           N

                     N is INTEGER
                     If JOB = 'L' or 'R', the number of columns of the matrix,
                     in which case N >= 0. Ignored if JOB = 'E'.

           D

                     D is REAL array, dimension (M) if JOB = 'E'
                                         dimension (min(M,N)) if JOB = 'L' or 'R'
                     The eigenvalues (if JOB = 'E') or singular values (if JOB =
                     'L' or 'R') of the matrix, in either increasing or decreasing
                     order. If singular values, they must be non-negative.

           SEP

                     SEP is REAL array, dimension (M) if JOB = 'E'
                                          dimension (min(M,N)) if JOB = 'L' or 'R'
                     The reciprocal condition numbers of the vectors.

           INFO

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine slaed0 (integer ICOMPQ, integer QSIZ, integer N, real, dimension( * ) D, real,
       dimension( * ) E, real, dimension( ldq, * ) Q, integer LDQ, real, dimension( ldqs, * )
       QSTORE, integer LDQS, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer
       INFO)
       SLAED0 used by SSTEDC. Computes all eigenvalues and corresponding eigenvectors of an
       unreduced symmetric tridiagonal matrix using the divide and conquer method.

       Purpose:

            SLAED0 computes all eigenvalues and corresponding eigenvectors of a
            symmetric tridiagonal matrix using the divide and conquer method.

       Parameters
           ICOMPQ

                     ICOMPQ is INTEGER
                     = 0:  Compute eigenvalues only.
                     = 1:  Compute eigenvectors of original dense symmetric matrix
                           also.  On entry, Q contains the orthogonal matrix used
                           to reduce the original matrix to tridiagonal form.
                     = 2:  Compute eigenvalues and eigenvectors of tridiagonal
                           matrix.

           QSIZ

                     QSIZ is INTEGER
                    The dimension of the orthogonal matrix used to reduce
                    the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.

           N

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

           D

                     D is REAL array, dimension (N)
                    On entry, the main diagonal of the tridiagonal matrix.
                    On exit, its eigenvalues.

           E

                     E is REAL array, dimension (N-1)
                    The off-diagonal elements of the tridiagonal matrix.
                    On exit, E has been destroyed.

           Q

                     Q is REAL array, dimension (LDQ, N)
                    On entry, Q must contain an N-by-N orthogonal matrix.
                    If ICOMPQ = 0    Q is not referenced.
                    If ICOMPQ = 1    On entry, Q is a subset of the columns of the
                                     orthogonal matrix used to reduce the full
                                     matrix to tridiagonal form corresponding to
                                     the subset of the full matrix which is being
                                     decomposed at this time.
                    If ICOMPQ = 2    On entry, Q will be the identity matrix.
                                     On exit, Q contains the eigenvectors of the
                                     tridiagonal matrix.

           LDQ

                     LDQ is INTEGER
                    The leading dimension of the array Q.  If eigenvectors are
                    desired, then  LDQ >= max(1,N).  In any case,  LDQ >= 1.

           QSTORE

                     QSTORE is REAL array, dimension (LDQS, N)
                    Referenced only when ICOMPQ = 1.  Used to store parts of
                    the eigenvector matrix when the updating matrix multiplies
                    take place.

           LDQS

                     LDQS is INTEGER
                    The leading dimension of the array QSTORE.  If ICOMPQ = 1,
                    then  LDQS >= max(1,N).  In any case,  LDQS >= 1.

           WORK

                     WORK is REAL array,
                    If ICOMPQ = 0 or 1, the dimension of WORK must be at least
                                1 + 3*N + 2*N*lg N + 3*N**2
                                ( lg( N ) = smallest integer k
                                            such that 2^k >= N )
                    If ICOMPQ = 2, the dimension of WORK must be at least
                                4*N + N**2.

           IWORK

                     IWORK is INTEGER array,
                    If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
                                   6 + 6*N + 5*N*lg N.
                                   ( lg( N ) = smallest integer k
                                               such that 2^k >= N )
                    If ICOMPQ = 2, the dimension of IWORK must be at least
                                   3 + 5*N.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit.
                     < 0:  if INFO = -i, the i-th argument had an illegal value.
                     > 0:  The algorithm failed to compute an eigenvalue while
                           working on the submatrix lying in rows and columns
                           INFO/(N+1) through mod(INFO,N+1).

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine slaed1 (integer N, real, dimension( * ) D, real, dimension( ldq, * ) Q, integer
       LDQ, integer, dimension( * ) INDXQ, real RHO, integer CUTPNT, real, dimension( * ) WORK,
       integer, dimension( * ) IWORK, integer INFO)
       SLAED1 used by SSTEDC. Computes the updated eigensystem of a diagonal matrix after
       modification by a rank-one symmetric matrix. Used when the original matrix is tridiagonal.

       Purpose:

            SLAED1 computes the updated eigensystem of a diagonal
            matrix after modification by a rank-one symmetric matrix.  This
            routine is used only for the eigenproblem which requires all
            eigenvalues and eigenvectors of a tridiagonal matrix.  SLAED7 handles
            the case in which eigenvalues only or eigenvalues and eigenvectors
            of a full symmetric matrix (which was reduced to tridiagonal form)
            are desired.

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

               where Z = Q**T*u, u is a vector of length N with ones in the
               CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.

               The eigenvectors of the original matrix are stored in Q, and the
               eigenvalues are in D.  The algorithm consists of three stages:

                  The first stage consists of deflating the size of the problem
                  when there are multiple eigenvalues or if there is a zero in
                  the Z vector.  For each such occurrence the dimension of the
                  secular equation problem is reduced by one.  This stage is
                  performed by the routine SLAED2.

                  The second stage consists of calculating the updated
                  eigenvalues. This is done by finding the roots of the secular
                  equation via the routine SLAED4 (as called by SLAED3).
                  This routine also calculates the eigenvectors of the current
                  problem.

                  The final stage consists of computing the updated eigenvectors
                  directly using the updated eigenvalues.  The eigenvectors for
                  the current problem are multiplied with the eigenvectors from
                  the overall problem.

       Parameters
           N

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

           D

                     D is REAL array, dimension (N)
                    On entry, the eigenvalues of the rank-1-perturbed matrix.
                    On exit, the eigenvalues of the repaired matrix.

           Q

                     Q is REAL array, dimension (LDQ,N)
                    On entry, the eigenvectors of the rank-1-perturbed matrix.
                    On exit, the eigenvectors of the repaired tridiagonal matrix.

           LDQ

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

           INDXQ

                     INDXQ is INTEGER array, dimension (N)
                    On entry, the permutation which separately sorts the two
                    subproblems in D into ascending order.
                    On exit, the permutation which will reintegrate the
                    subproblems back into sorted order,
                    i.e. D( INDXQ( I = 1, N ) ) will be in ascending order.

           RHO

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

           CUTPNT

                     CUTPNT is INTEGER
                    The location of the last eigenvalue in the leading sub-matrix.
                    min(1,N) <= CUTPNT <= N/2.

           WORK

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

           IWORK

                     IWORK is INTEGER array, dimension (4*N)

           INFO

                     INFO is INTEGER
                     = 0:  successful exit.
                     < 0:  if INFO = -i, the i-th argument had an illegal value.
                     > 0:  if INFO = 1, an eigenvalue did not converge

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine slaed2 (integer K, integer N, integer N1, real, dimension( * ) D, real, dimension(
       ldq, * ) Q, integer LDQ, integer, dimension( * ) INDXQ, real RHO, real, dimension( * ) Z,
       real, dimension( * ) DLAMDA, real, dimension( * ) W, real, dimension( * ) Q2, integer,
       dimension( * ) INDX, integer, dimension( * ) INDXC, integer, dimension( * ) INDXP,
       integer, dimension( * ) COLTYP, integer INFO)
       SLAED2 used by SSTEDC. Merges eigenvalues and deflates secular equation. Used when the
       original matrix is tridiagonal.

       Purpose:

            SLAED2 merges the two sets of eigenvalues together into a single
            sorted set.  Then it tries to deflate the size of the problem.
            There are two ways in which deflation can occur:  when two or more
            eigenvalues are close together or if there is a tiny entry in the
            Z vector.  For each such occurrence the order of the related secular
            equation problem is reduced by one.

       Parameters
           K

                     K is INTEGER
                    The number of non-deflated eigenvalues, and the order of the
                    related secular equation. 0 <= K <=N.

           N

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

           N1

                     N1 is INTEGER
                    The location of the last eigenvalue in the leading sub-matrix.
                    min(1,N) <= N1 <= N/2.

           D

                     D is REAL array, dimension (N)
                    On entry, D contains the eigenvalues of the two submatrices to
                    be combined.
                    On exit, D contains the trailing (N-K) updated eigenvalues
                    (those which were deflated) sorted into increasing order.

           Q

                     Q is REAL array, dimension (LDQ, N)
                    On entry, Q contains the eigenvectors of two submatrices in
                    the two square blocks with corners at (1,1), (N1,N1)
                    and (N1+1, N1+1), (N,N).
                    On exit, Q contains the trailing (N-K) updated eigenvectors
                    (those which were deflated) in its last N-K columns.

           LDQ

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

           INDXQ

                     INDXQ is INTEGER array, dimension (N)
                    The permutation which separately sorts the two sub-problems
                    in D into ascending order.  Note that elements in the second
                    half of this permutation must first have N1 added to their
                    values. Destroyed on exit.

           RHO

                     RHO is REAL
                    On entry, the off-diagonal element associated with the rank-1
                    cut which originally split the two submatrices which are now
                    being recombined.
                    On exit, RHO has been modified to the value required by
                    SLAED3.

           Z

                     Z is REAL array, dimension (N)
                    On entry, Z contains the updating vector (the last
                    row of the first sub-eigenvector matrix and the first row of
                    the second sub-eigenvector matrix).
                    On exit, the contents of Z have been destroyed by the updating
                    process.

           DLAMDA

                     DLAMDA is REAL array, dimension (N)
                    A copy of the first K eigenvalues which will be used by
                    SLAED3 to form the secular equation.

           W

                     W is REAL array, dimension (N)
                    The first k values of the final deflation-altered z-vector
                    which will be passed to SLAED3.

           Q2

                     Q2 is REAL array, dimension (N1**2+(N-N1)**2)
                    A copy of the first K eigenvectors which will be used by
                    SLAED3 in a matrix multiply (SGEMM) to solve for the new
                    eigenvectors.

           INDX

                     INDX is INTEGER array, dimension (N)
                    The permutation used to sort the contents of DLAMDA into
                    ascending order.

           INDXC

                     INDXC is INTEGER array, dimension (N)
                    The permutation used to arrange the columns of the deflated
                    Q matrix into three groups:  the first group contains non-zero
                    elements only at and above N1, the second contains
                    non-zero elements only below N1, and the third is dense.

           INDXP

                     INDXP is INTEGER array, dimension (N)
                    The permutation used to place deflated values of D at the end
                    of the array.  INDXP(1:K) points to the nondeflated D-values
                    and INDXP(K+1:N) points to the deflated eigenvalues.

           COLTYP

                     COLTYP is INTEGER array, dimension (N)
                    During execution, a label which will indicate which of the
                    following types a column in the Q2 matrix is:
                    1 : non-zero in the upper half only;
                    2 : dense;
                    3 : non-zero in the lower half only;
                    4 : deflated.
                    On exit, COLTYP(i) is the number of columns of type i,
                    for i=1 to 4 only.

           INFO

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine slaed3 (integer K, integer N, integer N1, real, dimension( * ) D, real, dimension(
       ldq, * ) Q, integer LDQ, real RHO, real, dimension( * ) DLAMDA, real, dimension( * ) Q2,
       integer, dimension( * ) INDX, integer, dimension( * ) CTOT, real, dimension( * ) W, real,
       dimension( * ) S, integer INFO)
       SLAED3 used by SSTEDC. Finds the roots of the secular equation and updates the
       eigenvectors. Used when the original matrix is tridiagonal.

       Purpose:

            SLAED3 finds the roots of the secular equation, as defined by the
            values in D, W, and RHO, between 1 and K.  It makes the
            appropriate calls to SLAED4 and then updates the eigenvectors by
            multiplying the matrix of eigenvectors of the pair of eigensystems
            being combined by the matrix of eigenvectors of the K-by-K system
            which is solved here.

            This code makes very mild assumptions about floating point
            arithmetic. It will work on machines with a guard digit in
            add/subtract, or on those binary machines without guard digits
            which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
            It could conceivably fail on hexadecimal or decimal machines
            without guard digits, but we know of none.

       Parameters
           K

                     K is INTEGER
                     The number of terms in the rational function to be solved by
                     SLAED4.  K >= 0.

           N

                     N is INTEGER
                     The number of rows and columns in the Q matrix.
                     N >= K (deflation may result in N>K).

           N1

                     N1 is INTEGER
                     The location of the last eigenvalue in the leading submatrix.
                     min(1,N) <= N1 <= N/2.

           D

                     D is REAL array, dimension (N)
                     D(I) contains the updated eigenvalues for
                     1 <= I <= K.

           Q

                     Q is REAL array, dimension (LDQ,N)
                     Initially the first K columns are used as workspace.
                     On output the columns 1 to K contain
                     the updated eigenvectors.

           LDQ

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

           RHO

                     RHO is REAL
                     The value of the parameter in the rank one update equation.
                     RHO >= 0 required.

           DLAMDA

                     DLAMDA is REAL array, dimension (K)
                     The first K elements of this array contain the old roots
                     of the deflated updating problem.  These are the poles
                     of the secular equation. May be changed on output by
                     having lowest order bit set to zero on Cray X-MP, Cray Y-MP,
                     Cray-2, or Cray C-90, as described above.

           Q2

                     Q2 is REAL array, dimension (LDQ2*N)
                     The first K columns of this matrix contain the non-deflated
                     eigenvectors for the split problem.

           INDX

                     INDX is INTEGER array, dimension (N)
                     The permutation used to arrange the columns of the deflated
                     Q matrix into three groups (see SLAED2).
                     The rows of the eigenvectors found by SLAED4 must be likewise
                     permuted before the matrix multiply can take place.

           CTOT

                     CTOT is INTEGER array, dimension (4)
                     A count of the total number of the various types of columns
                     in Q, as described in INDX.  The fourth column type is any
                     column which has been deflated.

           W

                     W is REAL array, dimension (K)
                     The first K elements of this array contain the components
                     of the deflation-adjusted updating vector. Destroyed on
                     output.

           S

                     S is REAL array, dimension (N1 + 1)*K
                     Will contain the eigenvectors of the repaired matrix which
                     will be multiplied by the previously accumulated eigenvectors
                     to update the system.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit.
                     < 0:  if INFO = -i, the i-th argument had an illegal value.
                     > 0:  if INFO = 1, an eigenvalue did not converge

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine slaed4 (integer N, integer I, real, dimension( * ) D, real, dimension( * ) Z, real,
       dimension( * ) DELTA, real RHO, real DLAM, integer INFO)
       SLAED4 used by SSTEDC. Finds a single root of the secular equation.

       Purpose:

            This subroutine computes the I-th updated eigenvalue of a symmetric
            rank-one modification to a diagonal matrix whose elements are
            given in the array d, and that

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

            and that RHO > 0.  This is arranged by the calling routine, and is
            no loss in generality.  The rank-one modified system is thus

                       diag( D )  +  RHO * Z * Z_transpose.

            where we assume the Euclidean norm of Z is 1.

            The method consists of approximating the rational functions in the
            secular equation by simpler interpolating rational functions.

       Parameters
           N

                     N is INTEGER
                    The length of all arrays.

           I

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

           D

                     D is REAL array, dimension (N)
                    The original eigenvalues.  It is assumed that they are in
                    order, D(I) < D(J)  for I < J.

           Z

                     Z is REAL array, dimension (N)
                    The components of the updating vector.

           DELTA

                     DELTA is REAL array, dimension (N)
                    If N > 2, DELTA contains (D(j) - lambda_I) in its  j-th
                    component.  If N = 1, then DELTA(1) = 1. If N = 2, see SLAED5
                    for detail. The vector DELTA contains the information necessary
                    to construct the eigenvectors by SLAED3 and SLAED9.

           RHO

                     RHO is REAL
                    The scalar in the symmetric updating formula.

           DLAM

                     DLAM is REAL
                    The computed lambda_I, the I-th updated eigenvalue.

           INFO

                     INFO is INTEGER
                    = 0:  successful exit
                    > 0:  if INFO = 1, the updating process failed.

       Internal Parameters:

             Logical variable ORGATI (origin-at-i?) is used for distinguishing
             whether D(i) or D(i+1) is treated as the origin.

                       ORGATI = .true.    origin at i
                       ORGATI = .false.   origin at i+1

              Logical variable SWTCH3 (switch-for-3-poles?) is for noting
              if we are working with THREE poles!

              MAXIT is the maximum number of iterations allowed for each
              eigenvalue.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine slaed5 (integer I, real, dimension( 2 ) D, real, dimension( 2 ) Z, real, dimension(
       2 ) DELTA, real RHO, real DLAM)
       SLAED5 used by SSTEDC. Solves the 2-by-2 secular equation.

       Purpose:

            This subroutine computes the I-th eigenvalue of a symmetric rank-one
            modification of a 2-by-2 diagonal matrix

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

            The diagonal elements in the array D are assumed to satisfy

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

            We also assume RHO > 0 and that the Euclidean norm of the vector
            Z is one.

       Parameters
           I

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

           D

                     D is REAL array, dimension (2)
                    The original eigenvalues.  We assume D(1) < D(2).

           Z

                     Z is REAL array, dimension (2)
                    The components of the updating vector.

           DELTA

                     DELTA is REAL array, dimension (2)
                    The vector DELTA contains the information necessary
                    to construct the eigenvectors.

           RHO

                     RHO is REAL
                    The scalar in the symmetric updating formula.

           DLAM

                     DLAM is REAL
                    The computed lambda_I, the I-th updated eigenvalue.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine slaed6 (integer KNITER, logical ORGATI, real RHO, real, dimension( 3 ) D, real,
       dimension( 3 ) Z, real FINIT, real TAU, integer INFO)
       SLAED6 used by SSTEDC. Computes one Newton step in solution of the secular equation.

       Purpose:

            SLAED6 computes the positive or negative root (closest to the origin)
            of
                             z(1)        z(2)        z(3)
            f(x) =   rho + --------- + ---------- + ---------
                            d(1)-x      d(2)-x      d(3)-x

            It is assumed that

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

            This routine will be called by SLAED4 when necessary. In most cases,
            the root sought is the smallest in magnitude, though it might not be
            in some extremely rare situations.

       Parameters
           KNITER

                     KNITER is INTEGER
                          Refer to SLAED4 for its significance.

           ORGATI

                     ORGATI is LOGICAL
                          If ORGATI is true, the needed root is between d(2) and
                          d(3); otherwise it is between d(1) and d(2).  See
                          SLAED4 for further details.

           RHO

                     RHO is REAL
                          Refer to the equation f(x) above.

           D

                     D is REAL array, dimension (3)
                          D satisfies d(1) < d(2) < d(3).

           Z

                     Z is REAL array, dimension (3)
                          Each of the elements in z must be positive.

           FINIT

                     FINIT is REAL
                          The value of f at 0. It is more accurate than the one
                          evaluated inside this routine (if someone wants to do
                          so).

           TAU

                     TAU is REAL
                          The root of the equation f(x).

           INFO

                     INFO is INTEGER
                          = 0: successful exit
                          > 0: if INFO = 1, failure to converge

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Further Details:

             10/02/03: This version has a few statements commented out for thread
             safety (machine parameters are computed on each entry). SJH.

             05/10/06: Modified from a new version of Ren-Cang Li, use
                Gragg-Thornton-Warner cubic convergent scheme for better stability.

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

   subroutine slaed7 (integer ICOMPQ, integer N, integer QSIZ, integer TLVLS, integer CURLVL,
       integer CURPBM, real, dimension( * ) D, real, dimension( ldq, * ) Q, integer LDQ, integer,
       dimension( * ) INDXQ, real RHO, integer CUTPNT, real, dimension( * ) QSTORE, integer,
       dimension( * ) QPTR, integer, dimension( * ) PRMPTR, integer, dimension( * ) PERM,
       integer, dimension( * ) GIVPTR, integer, dimension( 2, * ) GIVCOL, real, dimension( 2, * )
       GIVNUM, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)
       SLAED7 used by SSTEDC. Computes the updated eigensystem of a diagonal matrix after
       modification by a rank-one symmetric matrix. Used when the original matrix is dense.

       Purpose:

            SLAED7 computes the updated eigensystem of a diagonal
            matrix after modification by a rank-one symmetric matrix. This
            routine is used only for the eigenproblem which requires all
            eigenvalues and optionally eigenvectors of a dense symmetric matrix
            that has been reduced to tridiagonal form.  SLAED1 handles
            the case in which all eigenvalues and eigenvectors of a symmetric
            tridiagonal matrix are desired.

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

               where Z = Q**Tu, u is a vector of length N with ones in the
               CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.

               The eigenvectors of the original matrix are stored in Q, and the
               eigenvalues are in D.  The algorithm consists of three stages:

                  The first stage consists of deflating the size of the problem
                  when there are multiple eigenvalues or if there is a zero in
                  the Z vector.  For each such occurrence the dimension of the
                  secular equation problem is reduced by one.  This stage is
                  performed by the routine SLAED8.

                  The second stage consists of calculating the updated
                  eigenvalues. This is done by finding the roots of the secular
                  equation via the routine SLAED4 (as called by SLAED9).
                  This routine also calculates the eigenvectors of the current
                  problem.

                  The final stage consists of computing the updated eigenvectors
                  directly using the updated eigenvalues.  The eigenvectors for
                  the current problem are multiplied with the eigenvectors from
                  the overall problem.

       Parameters
           ICOMPQ

                     ICOMPQ is INTEGER
                     = 0:  Compute eigenvalues only.
                     = 1:  Compute eigenvectors of original dense symmetric matrix
                           also.  On entry, Q contains the orthogonal matrix used
                           to reduce the original matrix to tridiagonal form.

           N

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

           QSIZ

                     QSIZ is INTEGER
                    The dimension of the orthogonal matrix used to reduce
                    the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.

           TLVLS

                     TLVLS is INTEGER
                    The total number of merging levels in the overall divide and
                    conquer tree.

           CURLVL

                     CURLVL is INTEGER
                    The current level in the overall merge routine,
                    0 <= CURLVL <= TLVLS.

           CURPBM

                     CURPBM is INTEGER
                    The current problem in the current level in the overall
                    merge routine (counting from upper left to lower right).

           D

                     D is REAL array, dimension (N)
                    On entry, the eigenvalues of the rank-1-perturbed matrix.
                    On exit, the eigenvalues of the repaired matrix.

           Q

                     Q is REAL array, dimension (LDQ, N)
                    On entry, the eigenvectors of the rank-1-perturbed matrix.
                    On exit, the eigenvectors of the repaired tridiagonal matrix.

           LDQ

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

           INDXQ

                     INDXQ is INTEGER array, dimension (N)
                    The permutation which will reintegrate the subproblem just
                    solved back into sorted order, i.e., D( INDXQ( I = 1, N ) )
                    will be in ascending order.

           RHO

                     RHO is REAL
                    The subdiagonal element used to create the rank-1
                    modification.

           CUTPNT

                     CUTPNT is INTEGER
                    Contains the location of the last eigenvalue in the leading
                    sub-matrix.  min(1,N) <= CUTPNT <= N.

           QSTORE

                     QSTORE is REAL array, dimension (N**2+1)
                    Stores eigenvectors of submatrices encountered during
                    divide and conquer, packed together. QPTR points to
                    beginning of the submatrices.

           QPTR

                     QPTR is INTEGER array, dimension (N+2)
                    List of indices pointing to beginning of submatrices stored
                    in QSTORE. The submatrices are numbered starting at the
                    bottom left of the divide and conquer tree, from left to
                    right and bottom to top.

           PRMPTR

                     PRMPTR is INTEGER array, dimension (N lg N)
                    Contains a list of pointers which indicate where in PERM a
                    level's permutation is stored.  PRMPTR(i+1) - PRMPTR(i)
                    indicates the size of the permutation and also the size of
                    the full, non-deflated problem.

           PERM

                     PERM is INTEGER array, dimension (N lg N)
                    Contains the permutations (from deflation and sorting) to be
                    applied to each eigenblock.

           GIVPTR

                     GIVPTR is INTEGER array, dimension (N lg N)
                    Contains a list of pointers which indicate where in GIVCOL a
                    level's Givens rotations are stored.  GIVPTR(i+1) - GIVPTR(i)
                    indicates the number of Givens rotations.

           GIVCOL

                     GIVCOL is INTEGER array, dimension (2, N lg N)
                    Each pair of numbers indicates a pair of columns to take place
                    in a Givens rotation.

           GIVNUM

                     GIVNUM is REAL array, dimension (2, N lg N)
                    Each number indicates the S value to be used in the
                    corresponding Givens rotation.

           WORK

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

           IWORK

                     IWORK is INTEGER array, dimension (4*N)

           INFO

                     INFO is INTEGER
                     = 0:  successful exit.
                     < 0:  if INFO = -i, the i-th argument had an illegal value.
                     > 0:  if INFO = 1, an eigenvalue did not converge

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine slaed8 (integer ICOMPQ, integer K, integer N, integer QSIZ, real, dimension( * ) D,
       real, dimension( ldq, * ) Q, integer LDQ, integer, dimension( * ) INDXQ, real RHO, integer
       CUTPNT, real, dimension( * ) Z, real, dimension( * ) DLAMDA, real, dimension( ldq2, * )
       Q2, integer LDQ2, real, dimension( * ) W, integer, dimension( * ) PERM, integer GIVPTR,
       integer, dimension( 2, * ) GIVCOL, real, dimension( 2, * ) GIVNUM, integer, dimension( * )
       INDXP, integer, dimension( * ) INDX, integer INFO)
       SLAED8 used by SSTEDC. Merges eigenvalues and deflates secular equation. Used when the
       original matrix is dense.

       Purpose:

            SLAED8 merges the two sets of eigenvalues together into a single
            sorted set.  Then it tries to deflate the size of the problem.
            There are two ways in which deflation can occur:  when two or more
            eigenvalues are close together or if there is a tiny element in the
            Z vector.  For each such occurrence the order of the related secular
            equation problem is reduced by one.

       Parameters
           ICOMPQ

                     ICOMPQ is INTEGER
                     = 0:  Compute eigenvalues only.
                     = 1:  Compute eigenvectors of original dense symmetric matrix
                           also.  On entry, Q contains the orthogonal matrix used
                           to reduce the original matrix to tridiagonal form.

           K

                     K is INTEGER
                    The number of non-deflated eigenvalues, and the order of the
                    related secular equation.

           N

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

           QSIZ

                     QSIZ is INTEGER
                    The dimension of the orthogonal matrix used to reduce
                    the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.

           D

                     D is REAL array, dimension (N)
                    On entry, the eigenvalues of the two submatrices to be
                    combined.  On exit, the trailing (N-K) updated eigenvalues
                    (those which were deflated) sorted into increasing order.

           Q

                     Q is REAL array, dimension (LDQ,N)
                    If ICOMPQ = 0, Q is not referenced.  Otherwise,
                    on entry, Q contains the eigenvectors of the partially solved
                    system which has been previously updated in matrix
                    multiplies with other partially solved eigensystems.
                    On exit, Q contains the trailing (N-K) updated eigenvectors
                    (those which were deflated) in its last N-K columns.

           LDQ

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

           INDXQ

                     INDXQ is INTEGER array, dimension (N)
                    The permutation which separately sorts the two sub-problems
                    in D into ascending order.  Note that elements in the second
                    half of this permutation must first have CUTPNT added to
                    their values in order to be accurate.

           RHO

                     RHO is REAL
                    On entry, the off-diagonal element associated with the rank-1
                    cut which originally split the two submatrices which are now
                    being recombined.
                    On exit, RHO has been modified to the value required by
                    SLAED3.

           CUTPNT

                     CUTPNT is INTEGER
                    The location of the last eigenvalue in the leading
                    sub-matrix.  min(1,N) <= CUTPNT <= N.

           Z

                     Z is REAL array, dimension (N)
                    On entry, Z contains the updating vector (the last row of
                    the first sub-eigenvector matrix and the first row of the
                    second sub-eigenvector matrix).
                    On exit, the contents of Z are destroyed by the updating
                    process.

           DLAMDA

                     DLAMDA is REAL array, dimension (N)
                    A copy of the first K eigenvalues which will be used by
                    SLAED3 to form the secular equation.

           Q2

                     Q2 is REAL array, dimension (LDQ2,N)
                    If ICOMPQ = 0, Q2 is not referenced.  Otherwise,
                    a copy of the first K eigenvectors which will be used by
                    SLAED7 in a matrix multiply (SGEMM) to update the new
                    eigenvectors.

           LDQ2

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

           W

                     W is REAL array, dimension (N)
                    The first k values of the final deflation-altered z-vector and
                    will be passed to SLAED3.

           PERM

                     PERM is INTEGER array, dimension (N)
                    The permutations (from deflation and sorting) to be applied
                    to each eigenblock.

           GIVPTR

                     GIVPTR is INTEGER
                    The number of Givens rotations which took place in this
                    subproblem.

           GIVCOL

                     GIVCOL is INTEGER array, dimension (2, N)
                    Each pair of numbers indicates a pair of columns to take place
                    in a Givens rotation.

           GIVNUM

                     GIVNUM is REAL array, dimension (2, N)
                    Each number indicates the S value to be used in the
                    corresponding Givens rotation.

           INDXP

                     INDXP is INTEGER array, dimension (N)
                    The permutation used to place deflated values of D at the end
                    of the array.  INDXP(1:K) points to the nondeflated D-values
                    and INDXP(K+1:N) points to the deflated eigenvalues.

           INDX

                     INDX is INTEGER array, dimension (N)
                    The permutation used to sort the contents of D into ascending
                    order.

           INFO

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine slaed9 (integer K, integer KSTART, integer KSTOP, integer N, real, dimension( * )
       D, real, dimension( ldq, * ) Q, integer LDQ, real RHO, real, dimension( * ) DLAMDA, real,
       dimension( * ) W, real, dimension( lds, * ) S, integer LDS, integer INFO)
       SLAED9 used by SSTEDC. Finds the roots of the secular equation and updates the
       eigenvectors. Used when the original matrix is dense.

       Purpose:

            SLAED9 finds the roots of the secular equation, as defined by the
            values in D, Z, and RHO, between KSTART and KSTOP.  It makes the
            appropriate calls to SLAED4 and then stores the new matrix of
            eigenvectors for use in calculating the next level of Z vectors.

       Parameters
           K

                     K is INTEGER
                     The number of terms in the rational function to be solved by
                     SLAED4.  K >= 0.

           KSTART

                     KSTART is INTEGER

           KSTOP

                     KSTOP is INTEGER
                     The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP
                     are to be computed.  1 <= KSTART <= KSTOP <= K.

           N

                     N is INTEGER
                     The number of rows and columns in the Q matrix.
                     N >= K (delation may result in N > K).

           D

                     D is REAL array, dimension (N)
                     D(I) contains the updated eigenvalues
                     for KSTART <= I <= KSTOP.

           Q

                     Q is REAL array, dimension (LDQ,N)

           LDQ

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

           RHO

                     RHO is REAL
                     The value of the parameter in the rank one update equation.
                     RHO >= 0 required.

           DLAMDA

                     DLAMDA is REAL array, dimension (K)
                     The first K elements of this array contain the old roots
                     of the deflated updating problem.  These are the poles
                     of the secular equation.

           W

                     W is REAL array, dimension (K)
                     The first K elements of this array contain the components
                     of the deflation-adjusted updating vector.

           S

                     S is REAL array, dimension (LDS, K)
                     Will contain the eigenvectors of the repaired matrix which
                     will be stored for subsequent Z vector calculation and
                     multiplied by the previously accumulated eigenvectors
                     to update the system.

           LDS

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

           INFO

                     INFO is INTEGER
                     = 0:  successful exit.
                     < 0:  if INFO = -i, the i-th argument had an illegal value.
                     > 0:  if INFO = 1, an eigenvalue did not converge

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine slaeda (integer N, integer TLVLS, integer CURLVL, integer CURPBM, integer,
       dimension( * ) PRMPTR, integer, dimension( * ) PERM, integer, dimension( * ) GIVPTR,
       integer, dimension( 2, * ) GIVCOL, real, dimension( 2, * ) GIVNUM, real, dimension( * ) Q,
       integer, dimension( * ) QPTR, real, dimension( * ) Z, real, dimension( * ) ZTEMP, integer
       INFO)
       SLAEDA used by SSTEDC. Computes the Z vector determining the rank-one modification of the
       diagonal matrix. Used when the original matrix is dense.

       Purpose:

            SLAEDA computes the Z vector corresponding to the merge step in the
            CURLVLth step of the merge process with TLVLS steps for the CURPBMth
            problem.

       Parameters
           N

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

           TLVLS

                     TLVLS is INTEGER
                    The total number of merging levels in the overall divide and
                    conquer tree.

           CURLVL

                     CURLVL is INTEGER
                    The current level in the overall merge routine,
                    0 <= curlvl <= tlvls.

           CURPBM

                     CURPBM is INTEGER
                    The current problem in the current level in the overall
                    merge routine (counting from upper left to lower right).

           PRMPTR

                     PRMPTR is INTEGER array, dimension (N lg N)
                    Contains a list of pointers which indicate where in PERM a
                    level's permutation is stored.  PRMPTR(i+1) - PRMPTR(i)
                    indicates the size of the permutation and incidentally the
                    size of the full, non-deflated problem.

           PERM

                     PERM is INTEGER array, dimension (N lg N)
                    Contains the permutations (from deflation and sorting) to be
                    applied to each eigenblock.

           GIVPTR

                     GIVPTR is INTEGER array, dimension (N lg N)
                    Contains a list of pointers which indicate where in GIVCOL a
                    level's Givens rotations are stored.  GIVPTR(i+1) - GIVPTR(i)
                    indicates the number of Givens rotations.

           GIVCOL

                     GIVCOL is INTEGER array, dimension (2, N lg N)
                    Each pair of numbers indicates a pair of columns to take place
                    in a Givens rotation.

           GIVNUM

                     GIVNUM is REAL array, dimension (2, N lg N)
                    Each number indicates the S value to be used in the
                    corresponding Givens rotation.

           Q

                     Q is REAL array, dimension (N**2)
                    Contains the square eigenblocks from previous levels, the
                    starting positions for blocks are given by QPTR.

           QPTR

                     QPTR is INTEGER array, dimension (N+2)
                    Contains a list of pointers which indicate where in Q an
                    eigenblock is stored.  SQRT( QPTR(i+1) - QPTR(i) ) indicates
                    the size of the block.

           Z

                     Z is REAL array, dimension (N)
                    On output this vector contains the updating vector (the last
                    row of the first sub-eigenvector matrix and the first row of
                    the second sub-eigenvector matrix).

           ZTEMP

                     ZTEMP is REAL array, dimension (N)

           INFO

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine slagtf (integer N, real, dimension( * ) A, real LAMBDA, real, dimension( * ) B,
       real, dimension( * ) C, real TOL, real, dimension( * ) D, integer, dimension( * ) IN,
       integer INFO)
       SLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal
       matrix, and λ a scalar, using partial pivoting with row interchanges.

       Purpose:

            SLAGTF factorizes the matrix (T - lambda*I), where T is an n by n
            tridiagonal matrix and lambda is a scalar, as

               T - lambda*I = PLU,

            where P is a permutation matrix, L is a unit lower tridiagonal matrix
            with at most one non-zero sub-diagonal elements per column and U is
            an upper triangular matrix with at most two non-zero super-diagonal
            elements per column.

            The factorization is obtained by Gaussian elimination with partial
            pivoting and implicit row scaling.

            The parameter LAMBDA is included in the routine so that SLAGTF may
            be used, in conjunction with SLAGTS, to obtain eigenvectors of T by
            inverse iteration.

       Parameters
           N

                     N is INTEGER
                     The order of the matrix T.

           A

                     A is REAL array, dimension (N)
                     On entry, A must contain the diagonal elements of T.

                     On exit, A is overwritten by the n diagonal elements of the
                     upper triangular matrix U of the factorization of T.

           LAMBDA

                     LAMBDA is REAL
                     On entry, the scalar lambda.

           B

                     B is REAL array, dimension (N-1)
                     On entry, B must contain the (n-1) super-diagonal elements of
                     T.

                     On exit, B is overwritten by the (n-1) super-diagonal
                     elements of the matrix U of the factorization of T.

           C

                     C is REAL array, dimension (N-1)
                     On entry, C must contain the (n-1) sub-diagonal elements of
                     T.

                     On exit, C is overwritten by the (n-1) sub-diagonal elements
                     of the matrix L of the factorization of T.

           TOL

                     TOL is REAL
                     On entry, a relative tolerance used to indicate whether or
                     not the matrix (T - lambda*I) is nearly singular. TOL should
                     normally be chose as approximately the largest relative error
                     in the elements of T. For example, if the elements of T are
                     correct to about 4 significant figures, then TOL should be
                     set to about 5*10**(-4). If TOL is supplied as less than eps,
                     where eps is the relative machine precision, then the value
                     eps is used in place of TOL.

           D

                     D is REAL array, dimension (N-2)
                     On exit, D is overwritten by the (n-2) second super-diagonal
                     elements of the matrix U of the factorization of T.

           IN

                     IN is INTEGER array, dimension (N)
                     On exit, IN contains details of the permutation matrix P. If
                     an interchange occurred at the kth step of the elimination,
                     then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)
                     returns the smallest positive integer j such that

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

                     where norm( A(j) ) denotes the sum of the absolute values of
                     the jth row of the matrix A. If no such j exists then IN(n)
                     is returned as zero. If IN(n) is returned as positive, then a
                     diagonal element of U is small, indicating that
                     (T - lambda*I) is singular or nearly singular,

           INFO

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine slamrg (integer N1, integer N2, real, dimension( * ) A, integer STRD1, integer
       STRD2, integer, dimension( * ) INDEX)
       SLAMRG creates a permutation list to merge the entries of two independently sorted sets
       into a single set sorted in ascending order.

       Purpose:

            SLAMRG will create a permutation list which will merge the elements
            of A (which is composed of two independently sorted sets) into a
            single set which is sorted in ascending order.

       Parameters
           N1

                     N1 is INTEGER

           N2

                     N2 is INTEGER
                    These arguments contain the respective lengths of the two
                    sorted lists to be merged.

           A

                     A is REAL array, dimension (N1+N2)
                    The first N1 elements of A contain a list of numbers which
                    are sorted in either ascending or descending order.  Likewise
                    for the final N2 elements.

           STRD1

                     STRD1 is INTEGER

           STRD2

                     STRD2 is INTEGER
                    These are the strides to be taken through the array A.
                    Allowable strides are 1 and -1.  They indicate whether a
                    subset of A is sorted in ascending (STRDx = 1) or descending
                    (STRDx = -1) order.

           INDEX

                     INDEX is INTEGER array, dimension (N1+N2)
                    On exit this array will contain a permutation such that
                    if B( I ) = A( INDEX( I ) ) for I=1,N1+N2, then B will be
                    sorted in ascending order.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine slartgs (real X, real Y, real SIGMA, real CS, real SN)
       SLARTGS generates a plane rotation designed to introduce a bulge in implicit QR iteration
       for the bidiagonal SVD problem.

       Purpose:

            SLARTGS generates a plane rotation designed to introduce a bulge in
            Golub-Reinsch-style implicit QR iteration for the bidiagonal SVD
            problem. X and Y are the top-row entries, and SIGMA is the shift.
            The computed CS and SN define a plane rotation satisfying

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

            with R nonnegative.  If X^2 - SIGMA and X * Y are 0, then the
            rotation is by PI/2.

       Parameters
           X

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

           Y

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

           SIGMA

                     SIGMA is REAL
                     The shift.

           CS

                     CS is REAL
                     The cosine of the rotation.

           SN

                     SN is REAL
                     The sine of the rotation.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine slasq1 (integer N, real, dimension( * ) D, real, dimension( * ) E, real, dimension(
       * ) WORK, integer INFO)
       SLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.

       Purpose:

            SLASQ1 computes the singular values of a real N-by-N bidiagonal
            matrix with diagonal D and off-diagonal E. The singular values
            are computed to high relative accuracy, in the absence of
            denormalization, underflow and overflow. The algorithm was first
            presented in

            'Accurate singular values and differential qd algorithms' by K. V.
            Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
            1994,

            and the present implementation is described in 'An implementation of
            the dqds Algorithm (Positive Case)', LAPACK Working Note.

       Parameters
           N

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

           D

                     D is REAL array, dimension (N)
                   On entry, D contains the diagonal elements of the
                   bidiagonal matrix whose SVD is desired. On normal exit,
                   D contains the singular values in decreasing order.

           E

                     E is REAL array, dimension (N)
                   On entry, elements E(1:N-1) contain the off-diagonal elements
                   of the bidiagonal matrix whose SVD is desired.
                   On exit, E is overwritten.

           WORK

                     WORK is REAL array, dimension (4*N)

           INFO

                     INFO is INTEGER
                   = 0: successful exit
                   < 0: if INFO = -i, the i-th argument had an illegal value
                   > 0: the algorithm failed
                        = 1, a split was marked by a positive value in E
                        = 2, current block of Z not diagonalized after 100*N
                             iterations (in inner while loop)  On exit D and E
                             represent a matrix with the same singular values
                             which the calling subroutine could use to finish the
                             computation, or even feed back into SLASQ1
                        = 3, termination criterion of outer while loop not met
                             (program created more than N unreduced blocks)

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine slasq2 (integer N, real, dimension( * ) Z, integer INFO)
       SLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix
       associated with the qd Array Z to high relative accuracy. Used by sbdsqr and sstegr.

       Purpose:

            SLASQ2 computes all the eigenvalues of the symmetric positive
            definite tridiagonal matrix associated with the qd array Z to high
            relative accuracy are computed to high relative accuracy, in the
            absence of denormalization, underflow and overflow.

            To see the relation of Z to the tridiagonal matrix, let L be a
            unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
            let U be an upper bidiagonal matrix with 1's above and diagonal
            Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
            symmetric tridiagonal to which it is similar.

            Note : SLASQ2 defines a logical variable, IEEE, which is true
            on machines which follow ieee-754 floating-point standard in their
            handling of infinities and NaNs, and false otherwise. This variable
            is passed to SLASQ3.

       Parameters
           N

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

           Z

                     Z is REAL array, dimension ( 4*N )
                   On entry Z holds the qd array. On exit, entries 1 to N hold
                   the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
                   trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
                   N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
                   holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
                   shifts that failed.

           INFO

                     INFO is INTEGER
                   = 0: successful exit
                   < 0: if the i-th argument is a scalar and had an illegal
                        value, then INFO = -i, if the i-th argument is an
                        array and the j-entry had an illegal value, then
                        INFO = -(i*100+j)
                   > 0: the algorithm failed
                         = 1, a split was marked by a positive value in E
                         = 2, current block of Z not diagonalized after 100*N
                              iterations (in inner while loop).  On exit Z holds
                              a qd array with the same eigenvalues as the given Z.
                         = 3, termination criterion of outer while loop not met
                              (program created more than N unreduced blocks)

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Further Details:

             Local Variables: I0:N0 defines a current unreduced segment of Z.
             The shifts are accumulated in SIGMA. Iteration count is in ITER.
             Ping-pong is controlled by PP (alternates between 0 and 1).

   subroutine slasq3 (integer I0, integer N0, real, dimension( * ) Z, integer PP, real DMIN, real
       SIGMA, real DESIG, real QMAX, integer NFAIL, integer ITER, integer NDIV, logical IEEE,
       integer TTYPE, real DMIN1, real DMIN2, real DN, real DN1, real DN2, real G, real TAU)
       SLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr.

       Purpose:

            SLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.
            In case of failure it changes shifts, and tries again until output
            is positive.

       Parameters
           I0

                     I0 is INTEGER
                    First index.

           N0

                     N0 is INTEGER
                    Last index.

           Z

                     Z is REAL array, dimension ( 4*N0 )
                    Z holds the qd array.

           PP

                     PP is INTEGER
                    PP=0 for ping, PP=1 for pong.
                    PP=2 indicates that flipping was applied to the Z array
                    and that the initial tests for deflation should not be
                    performed.

           DMIN

                     DMIN is REAL
                    Minimum value of d.

           SIGMA

                     SIGMA is REAL
                    Sum of shifts used in current segment.

           DESIG

                     DESIG is REAL
                    Lower order part of SIGMA

           QMAX

                     QMAX is REAL
                    Maximum value of q.

           NFAIL

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

           ITER

                     ITER is INTEGER
                    Increment ITER by 1 for each iteration.

           NDIV

                     NDIV is INTEGER
                    Increment NDIV by 1 for each division.

           IEEE

                     IEEE is LOGICAL
                    Flag for IEEE or non IEEE arithmetic (passed to SLASQ5).

           TTYPE

                     TTYPE is INTEGER
                    Shift type.

           DMIN1

                     DMIN1 is REAL

           DMIN2

                     DMIN2 is REAL

           DN

                     DN is REAL

           DN1

                     DN1 is REAL

           DN2

                     DN2 is REAL

           G

                     G is REAL

           TAU

                     TAU is REAL

                    These are passed as arguments in order to save their values
                    between calls to SLASQ3.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine slasq4 (integer I0, integer N0, real, dimension( * ) Z, integer PP, integer N0IN,
       real DMIN, real DMIN1, real DMIN2, real DN, real DN1, real DN2, real TAU, integer TTYPE,
       real G)
       SLASQ4 computes an approximation to the smallest eigenvalue using values of d from the
       previous transform. Used by sbdsqr.

       Purpose:

            SLASQ4 computes an approximation TAU to the smallest eigenvalue
            using values of d from the previous transform.

       Parameters
           I0

                     I0 is INTEGER
                   First index.

           N0

                     N0 is INTEGER
                   Last index.

           Z

                     Z is REAL array, dimension ( 4*N0 )
                   Z holds the qd array.

           PP

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

           N0IN

                     N0IN is INTEGER
                   The value of N0 at start of EIGTEST.

           DMIN

                     DMIN is REAL
                   Minimum value of d.

           DMIN1

                     DMIN1 is REAL
                   Minimum value of d, excluding D( N0 ).

           DMIN2

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

           DN

                     DN is REAL
                   d(N)

           DN1

                     DN1 is REAL
                   d(N-1)

           DN2

                     DN2 is REAL
                   d(N-2)

           TAU

                     TAU is REAL
                   This is the shift.

           TTYPE

                     TTYPE is INTEGER
                   Shift type.

           G

                     G is REAL
                   G is passed as an argument in order to save its value between
                   calls to SLASQ4.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Further Details:

             CNST1 = 9/16

   subroutine slasq5 (integer I0, integer N0, real, dimension( * ) Z, integer PP, real TAU, real
       SIGMA, real DMIN, real DMIN1, real DMIN2, real DN, real DNM1, real DNM2, logical IEEE,
       real EPS)
        SLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.

       Purpose:

            SLASQ5 computes one dqds transform in ping-pong form, one
            version for IEEE machines another for non IEEE machines.

       Parameters
           I0

                     I0 is INTEGER
                   First index.

           N0

                     N0 is INTEGER
                   Last index.

           Z

                     Z is REAL array, dimension ( 4*N )
                   Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
                   an extra argument.

           PP

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

           TAU

                     TAU is REAL
                   This is the shift.

           SIGMA

                     SIGMA is REAL
                   This is the accumulated shift up to this step.

           DMIN

                     DMIN is REAL
                   Minimum value of d.

           DMIN1

                     DMIN1 is REAL
                   Minimum value of d, excluding D( N0 ).

           DMIN2

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

           DN

                     DN is REAL
                   d(N0), the last value of d.

           DNM1

                     DNM1 is REAL
                   d(N0-1).

           DNM2

                     DNM2 is REAL
                   d(N0-2).

           IEEE

                     IEEE is LOGICAL
                   Flag for IEEE or non IEEE arithmetic.

           EPS

                    EPS is REAL
                   This is the value of epsilon used.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

       Purpose:

            SLASQ6 computes one dqd (shift equal to zero) transform in
            ping-pong form, with protection against underflow and overflow.

       Parameters
           I0

                     I0 is INTEGER
                   First index.

           N0

                     N0 is INTEGER
                   Last index.

           Z

                     Z is REAL array, dimension ( 4*N )
                   Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
                   an extra argument.

           PP

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

           DMIN

                     DMIN is REAL
                   Minimum value of d.

           DMIN1

                     DMIN1 is REAL
                   Minimum value of d, excluding D( N0 ).

           DMIN2

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

           DN

                     DN is REAL
                   d(N0), the last value of d.

           DNM1

                     DNM1 is REAL
                   d(N0-1).

           DNM2

                     DNM2 is REAL
                   d(N0-2).

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine slasrt (character ID, integer N, real, dimension( * ) D, integer INFO)
       SLASRT sorts numbers in increasing or decreasing order.

       Purpose:

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

            Use Quick Sort, reverting to Insertion sort on arrays of
            size <= 20. Dimension of STACK limits N to about 2**32.

       Parameters
           ID

                     ID is CHARACTER*1
                     = 'I': sort D in increasing order;
                     = 'D': sort D in decreasing order.

           N

                     N is INTEGER
                     The length of the array D.

           D

                     D is REAL array, dimension (N)
                     On entry, the array to be sorted.
                     On exit, D has been sorted into increasing order
                     (D(1) <= ... <= D(N) ) or into decreasing order
                     (D(1) >= ... >= D(N) ), depending on ID.

           INFO

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine spttrf (integer N, real, dimension( * ) D, real, dimension( * ) E, integer INFO)
       SPTTRF

       Purpose:

            SPTTRF computes the L*D*L**T factorization of a real symmetric
            positive definite tridiagonal matrix A.  The factorization may also
            be regarded as having the form A = U**T*D*U.

       Parameters
           N

                     N is INTEGER
                     The order of the matrix A.  N >= 0.

           D

                     D is REAL array, dimension (N)
                     On entry, the n diagonal elements of the tridiagonal matrix
                     A.  On exit, the n diagonal elements of the diagonal matrix
                     D from the L*D*L**T factorization of A.

           E

                     E is REAL array, dimension (N-1)
                     On entry, the (n-1) subdiagonal elements of the tridiagonal
                     matrix A.  On exit, the (n-1) subdiagonal elements of the
                     unit bidiagonal factor L from the L*D*L**T factorization of A.
                     E can also be regarded as the superdiagonal of the unit
                     bidiagonal factor U from the U**T*D*U factorization of A.

           INFO

                     INFO is INTEGER
                     = 0: successful exit
                     < 0: if INFO = -k, the k-th argument had an illegal value
                     > 0: if INFO = k, the leading minor of order k is not
                          positive definite; if k < N, the factorization could not
                          be completed, while if k = N, the factorization was
                          completed, but D(N) <= 0.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

       Purpose:

            SSTEBZ computes the eigenvalues of a symmetric tridiagonal
            matrix T.  The user may ask for all eigenvalues, all eigenvalues
            in the half-open interval (VL, VU], or the IL-th through IU-th
            eigenvalues.

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

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

       Parameters
           RANGE

                     RANGE is CHARACTER*1
                     = 'A': ('All')   all eigenvalues will be found.
                     = 'V': ('Value') all eigenvalues in the half-open interval
                                      (VL, VU] will be found.
                     = 'I': ('Index') the IL-th through IU-th eigenvalues (of the
                                      entire matrix) will be found.

           ORDER

                     ORDER is CHARACTER*1
                     = 'B': ('By Block') the eigenvalues will be grouped by
                                         split-off block (see IBLOCK, ISPLIT) and
                                         ordered from smallest to largest within
                                         the block.
                     = 'E': ('Entire matrix')
                                         the eigenvalues for the entire matrix
                                         will be ordered from smallest to
                                         largest.

           N

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

           VL

                     VL is REAL

                     If RANGE='V', the lower bound of the interval to
                     be searched for eigenvalues.  Eigenvalues less than or equal
                     to VL, or greater than VU, will not be returned.  VL < VU.
                     Not referenced if RANGE = 'A' or 'I'.

           VU

                     VU is REAL

                     If RANGE='V', the upper bound of the interval to
                     be searched for eigenvalues.  Eigenvalues less than or equal
                     to VL, or greater than VU, will not be returned.  VL < VU.
                     Not referenced if RANGE = 'A' or 'I'.

           IL

                     IL is INTEGER

                     If RANGE='I', the index of the
                     smallest eigenvalue to be returned.
                     1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
                     Not referenced if RANGE = 'A' or 'V'.

           IU

                     IU is INTEGER

                     If RANGE='I', the index of the
                     largest eigenvalue to be returned.
                     1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
                     Not referenced if RANGE = 'A' or 'V'.

           ABSTOL

                     ABSTOL is REAL
                     The absolute tolerance for the eigenvalues.  An eigenvalue
                     (or cluster) is considered to be located if it has been
                     determined to lie in an interval whose width is ABSTOL or
                     less.  If ABSTOL is less than or equal to zero, then ULP*|T|
                     will be used, where |T| means the 1-norm of T.

                     Eigenvalues will be computed most accurately when ABSTOL is
                     set to twice the underflow threshold 2*SLAMCH('S'), not zero.

           D

                     D is REAL array, dimension (N)
                     The n diagonal elements of the tridiagonal matrix T.

           E

                     E is REAL array, dimension (N-1)
                     The (n-1) off-diagonal elements of the tridiagonal matrix T.

           M

                     M is INTEGER
                     The actual number of eigenvalues found. 0 <= M <= N.
                     (See also the description of INFO=2,3.)

           NSPLIT

                     NSPLIT is INTEGER
                     The number of diagonal blocks in the matrix T.
                     1 <= NSPLIT <= N.

           W

                     W is REAL array, dimension (N)
                     On exit, the first M elements of W will contain the
                     eigenvalues.  (SSTEBZ may use the remaining N-M elements as
                     workspace.)

           IBLOCK

                     IBLOCK is INTEGER array, dimension (N)
                     At each row/column j where E(j) is zero or small, the
                     matrix T is considered to split into a block diagonal
                     matrix.  On exit, if INFO = 0, IBLOCK(i) specifies to which
                     block (from 1 to the number of blocks) the eigenvalue W(i)
                     belongs.  (SSTEBZ may use the remaining N-M elements as
                     workspace.)

           ISPLIT

                     ISPLIT is INTEGER array, dimension (N)
                     The splitting points, at which T breaks up into submatrices.
                     The first submatrix consists of rows/columns 1 to ISPLIT(1),
                     the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
                     etc., and the NSPLIT-th consists of rows/columns
                     ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
                     (Only the first NSPLIT elements will actually be used, but
                     since the user cannot know a priori what value NSPLIT will
                     have, N words must be reserved for ISPLIT.)

           WORK

                     WORK is REAL array, dimension (4*N)

           IWORK

                     IWORK is INTEGER array, dimension (3*N)

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  if INFO = -i, the i-th argument had an illegal value
                     > 0:  some or all of the eigenvalues failed to converge or
                           were not computed:
                           =1 or 3: Bisection failed to converge for some
                                   eigenvalues; these eigenvalues are flagged by a
                                   negative block number.  The effect is that the
                                   eigenvalues may not be as accurate as the
                                   absolute and relative tolerances.  This is
                                   generally caused by unexpectedly inaccurate
                                   arithmetic.
                           =2 or 3: RANGE='I' only: Not all of the eigenvalues
                                   IL:IU were found.
                                   Effect: M < IU+1-IL
                                   Cause:  non-monotonic arithmetic, causing the
                                           Sturm sequence to be non-monotonic.
                                   Cure:   recalculate, using RANGE='A', and pick
                                           out eigenvalues IL:IU.  In some cases,
                                           increasing the PARAMETER 'FUDGE' may
                                           make things work.
                           = 4:    RANGE='I', and the Gershgorin interval
                                   initially used was too small.  No eigenvalues
                                   were computed.
                                   Probable cause: your machine has sloppy
                                                   floating-point arithmetic.
                                   Cure: Increase the PARAMETER 'FUDGE',
                                         recompile, and try again.

       Internal Parameters:

             RELFAC  REAL, default = 2.0e0
                     The relative tolerance.  An interval (a,b] lies within
                     'relative tolerance' if  b-a < RELFAC*ulp*max(|a|,|b|),
                     where 'ulp' is the machine precision (distance from 1 to
                     the next larger floating point number.)

             FUDGE   REAL, default = 2
                     A 'fudge factor' to widen the Gershgorin intervals.  Ideally,
                     a value of 1 should work, but on machines with sloppy
                     arithmetic, this needs to be larger.  The default for
                     publicly released versions should be large enough to handle
                     the worst machine around.  Note that this has no effect
                     on accuracy of the solution.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine sstedc (character COMPZ, integer N, real, dimension( * ) D, real, dimension( * ) E,
       real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer LWORK,
       integer, dimension( * ) IWORK, integer LIWORK, integer INFO)
       SSTEDC

       Purpose:

            SSTEDC computes all eigenvalues and, optionally, eigenvectors of a
            symmetric tridiagonal matrix using the divide and conquer method.
            The eigenvectors of a full or band real symmetric matrix can also be
            found if SSYTRD or SSPTRD or SSBTRD has been used to reduce this
            matrix to tridiagonal form.

            This code makes very mild assumptions about floating point
            arithmetic. It will work on machines with a guard digit in
            add/subtract, or on those binary machines without guard digits
            which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
            It could conceivably fail on hexadecimal or decimal machines
            without guard digits, but we know of none.  See SLAED3 for details.

       Parameters
           COMPZ

                     COMPZ is CHARACTER*1
                     = 'N':  Compute eigenvalues only.
                     = 'I':  Compute eigenvectors of tridiagonal matrix also.
                     = 'V':  Compute eigenvectors of original dense symmetric
                             matrix also.  On entry, Z contains the orthogonal
                             matrix used to reduce the original matrix to
                             tridiagonal form.

           N

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

           D

                     D is REAL array, dimension (N)
                     On entry, the diagonal elements of the tridiagonal matrix.
                     On exit, if INFO = 0, the eigenvalues in ascending order.

           E

                     E is REAL array, dimension (N-1)
                     On entry, the subdiagonal elements of the tridiagonal matrix.
                     On exit, E has been destroyed.

           Z

                     Z is REAL array, dimension (LDZ,N)
                     On entry, if COMPZ = 'V', then Z contains the orthogonal
                     matrix used in the reduction to tridiagonal form.
                     On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
                     orthonormal eigenvectors of the original symmetric matrix,
                     and if COMPZ = 'I', Z contains the orthonormal eigenvectors
                     of the symmetric tridiagonal matrix.
                     If  COMPZ = 'N', then Z is not referenced.

           LDZ

                     LDZ is INTEGER
                     The leading dimension of the array Z.  LDZ >= 1.
                     If eigenvectors are desired, then LDZ >= max(1,N).

           WORK

                     WORK is REAL array, dimension (MAX(1,LWORK))
                     On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

           LWORK

                     LWORK is INTEGER
                     The dimension of the array WORK.
                     If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
                     If COMPZ = 'V' and N > 1 then LWORK must be at least
                                    ( 1 + 3*N + 2*N*lg N + 4*N**2 ),
                                    where lg( N ) = smallest integer k such
                                    that 2**k >= N.
                     If COMPZ = 'I' and N > 1 then LWORK must be at least
                                    ( 1 + 4*N + N**2 ).
                     Note that for COMPZ = 'I' or 'V', then if N is less than or
                     equal to the minimum divide size, usually 25, then LWORK need
                     only be max(1,2*(N-1)).

                     If LWORK = -1, then a workspace query is assumed; the routine
                     only calculates the optimal size of the WORK array, returns
                     this value as the first entry of the WORK array, and no error
                     message related to LWORK is issued by XERBLA.

           IWORK

                     IWORK is INTEGER array, dimension (MAX(1,LIWORK))
                     On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.

           LIWORK

                     LIWORK is INTEGER
                     The dimension of the array IWORK.
                     If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
                     If COMPZ = 'V' and N > 1 then LIWORK must be at least
                                    ( 6 + 6*N + 5*N*lg N ).
                     If COMPZ = 'I' and N > 1 then LIWORK must be at least
                                    ( 3 + 5*N ).
                     Note that for COMPZ = 'I' or 'V', then if N is less than or
                     equal to the minimum divide size, usually 25, then LIWORK
                     need only be 1.

                     If LIWORK = -1, then a workspace query is assumed; the
                     routine only calculates the optimal size of the IWORK array,
                     returns this value as the first entry of the IWORK array, and
                     no error message related to LIWORK is issued by XERBLA.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit.
                     < 0:  if INFO = -i, the i-th argument had an illegal value.
                     > 0:  The algorithm failed to compute an eigenvalue while
                           working on the submatrix lying in rows and columns
                           INFO/(N+1) through mod(INFO,N+1).

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

   subroutine ssteqr (character COMPZ, integer N, real, dimension( * ) D, real, dimension( * ) E,
       real, dimension( ldz, * ) Z, integer LDZ, real, dimension( * ) WORK, integer INFO)
       SSTEQR

       Purpose:

            SSTEQR computes all eigenvalues and, optionally, eigenvectors of a
            symmetric tridiagonal matrix using the implicit QL or QR method.
            The eigenvectors of a full or band symmetric matrix can also be found
            if SSYTRD or SSPTRD or SSBTRD has been used to reduce this matrix to
            tridiagonal form.

       Parameters
           COMPZ

                     COMPZ is CHARACTER*1
                     = 'N':  Compute eigenvalues only.
                     = 'V':  Compute eigenvalues and eigenvectors of the original
                             symmetric matrix.  On entry, Z must contain the
                             orthogonal matrix used to reduce the original matrix
                             to tridiagonal form.
                     = 'I':  Compute eigenvalues and eigenvectors of the
                             tridiagonal matrix.  Z is initialized to the identity
                             matrix.

           N

                     N is INTEGER
                     The order of the matrix.  N >= 0.

           D

                     D is REAL array, dimension (N)
                     On entry, the diagonal elements of the tridiagonal matrix.
                     On exit, if INFO = 0, the eigenvalues in ascending order.

           E

                     E is REAL array, dimension (N-1)
                     On entry, the (n-1) subdiagonal elements of the tridiagonal
                     matrix.
                     On exit, E has been destroyed.

           Z

                     Z is REAL array, dimension (LDZ, N)
                     On entry, if  COMPZ = 'V', then Z contains the orthogonal
                     matrix used in the reduction to tridiagonal form.
                     On exit, if INFO = 0, then if  COMPZ = 'V', Z contains the
                     orthonormal eigenvectors of the original symmetric matrix,
                     and if COMPZ = 'I', Z contains the orthonormal eigenvectors
                     of the symmetric tridiagonal matrix.
                     If COMPZ = 'N', then Z is not referenced.

           LDZ

                     LDZ is INTEGER
                     The leading dimension of the array Z.  LDZ >= 1, and if
                     eigenvectors are desired, then  LDZ >= max(1,N).

           WORK

                     WORK is REAL array, dimension (max(1,2*N-2))
                     If COMPZ = 'N', then WORK is not referenced.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  if INFO = -i, the i-th argument had an illegal value
                     > 0:  the algorithm has failed to find all the eigenvalues in
                           a total of 30*N iterations; if INFO = i, then i
                           elements of E have not converged to zero; on exit, D
                           and E contain the elements of a symmetric tridiagonal
                           matrix which is orthogonally similar to the original
                           matrix.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine ssterf (integer N, real, dimension( * ) D, real, dimension( * ) E, integer INFO)
       SSTERF

       Purpose:

            SSTERF computes all eigenvalues of a symmetric tridiagonal matrix
            using the Pal-Walker-Kahan variant of the QL or QR algorithm.

       Parameters
           N

                     N is INTEGER
                     The order of the matrix.  N >= 0.

           D

                     D is REAL array, dimension (N)
                     On entry, the n diagonal elements of the tridiagonal matrix.
                     On exit, if INFO = 0, the eigenvalues in ascending order.

           E

                     E is REAL array, dimension (N-1)
                     On entry, the (n-1) subdiagonal elements of the tridiagonal
                     matrix.
                     On exit, E has been destroyed.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  if INFO = -i, the i-th argument had an illegal value
                     > 0:  the algorithm failed to find all of the eigenvalues in
                           a total of 30*N iterations; if INFO = i, then i
                           elements of E have not converged to zero.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

Author

       Generated automatically by Doxygen for LAPACK from the source code.