Provided by: liblapack-doc_3.8.0-2_all bug

NAME

       doubleOTHERauxiliary

SYNOPSIS

   Functions
       subroutine dlabrd (M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y, LDY)
           DLABRD reduces the first nb rows and columns of a general matrix to a bidiagonal form.
       subroutine dlacn2 (N, V, X, ISGN, EST, KASE, ISAVE)
           DLACN2 estimates the 1-norm of a square matrix, using reverse communication for
           evaluating matrix-vector products.
       subroutine dlacon (N, V, X, ISGN, EST, KASE)
           DLACON estimates the 1-norm of a square matrix, using reverse communication for
           evaluating matrix-vector products.
       subroutine dladiv (A, B, C, D, P, Q)
           DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
       subroutine dladiv1 (A, B, C, D, P, Q)
       double precision function dladiv2 (A, B, C, D, R, T)
       subroutine dlaein (RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B, LDB, WORK, EPS3, SMLNUM,
           BIGNUM, INFO)
           DLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by
           inverse iteration.
       subroutine dlaexc (WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK, INFO)
           DLAEXC swaps adjacent diagonal blocks of a real upper quasi-triangular matrix in Schur
           canonical form, by an orthogonal similarity transformation.
       subroutine dlag2 (A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, WI)
           DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with
           scaling as necessary to avoid over-/underflow.
       subroutine dlag2s (M, N, A, LDA, SA, LDSA, INFO)
           DLAG2S converts a double precision matrix to a single precision matrix.
       subroutine dlags2 (UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV, SNV, CSQ, SNQ)
           DLAGS2 computes 2-by-2 orthogonal matrices U, V, and Q, and applies them to matrices A
           and B such that the rows of the transformed A and B are parallel.
       subroutine dlagtm (TRANS, N, NRHS, ALPHA, DL, D, DU, X, LDX, BETA, B, LDB)
           DLAGTM performs a matrix-matrix product of the form C = αAB+βC, where A is a
           tridiagonal matrix, B and C are rectangular matrices, and α and β are scalars, which
           may be 0, 1, or -1.
       subroutine dlagv2 (A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL, CSR, SNR)
           DLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil
           (A,B) where B is upper triangular.
       subroutine dlahqr (WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, INFO)
           DLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,
           using the double-shift/single-shift QR algorithm.
       subroutine dlahr2 (N, K, NB, A, LDA, TAU, T, LDT, Y, LDY)
           DLAHR2 reduces the specified number of first columns of a general rectangular matrix A
           so that elements below the specified subdiagonal are zero, and returns auxiliary
           matrices which are needed to apply the transformation to the unreduced part of A.
       subroutine dlaic1 (JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
           DLAIC1 applies one step of incremental condition estimation.
       subroutine dlaln2 (LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, LDB, WR, WI, X, LDX,
           SCALE, XNORM, INFO)
           DLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
       double precision function dlangt (NORM, N, DL, D, DU)
           DLANGT returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest
           absolute value of any element of a general tridiagonal matrix.
       double precision function dlanhs (NORM, N, A, LDA, WORK)
           DLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest
           absolute value of any element of an upper Hessenberg matrix.
       double precision function dlansb (NORM, UPLO, N, K, AB, LDAB, WORK)
           DLANSB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,
           or the element of largest absolute value of a symmetric band matrix.
       double precision function dlansp (NORM, UPLO, N, AP, WORK)
           DLANSP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,
           or the element of largest absolute value of a symmetric matrix supplied in packed
           form.
       double precision function dlantb (NORM, UPLO, DIAG, N, K, AB, LDAB, WORK)
           DLANTB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,
           or the element of largest absolute value of a triangular band matrix.
       double precision function dlantp (NORM, UPLO, DIAG, N, AP, WORK)
           DLANTP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,
           or the element of largest absolute value of a triangular matrix supplied in packed
           form.
       double precision function dlantr (NORM, UPLO, DIAG, M, N, A, LDA, WORK)
           DLANTR returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,
           or the element of largest absolute value of a trapezoidal or triangular matrix.
       subroutine dlanv2 (A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN)
           DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in
           standard form.
       subroutine dlapll (N, X, INCX, Y, INCY, SSMIN)
           DLAPLL measures the linear dependence of two vectors.
       subroutine dlapmr (FORWRD, M, N, X, LDX, K)
           DLAPMR rearranges rows of a matrix as specified by a permutation vector.
       subroutine dlapmt (FORWRD, M, N, X, LDX, K)
           DLAPMT performs a forward or backward permutation of the columns of a matrix.
       subroutine dlaqp2 (M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, WORK)
           DLAQP2 computes a QR factorization with column pivoting of the matrix block.
       subroutine dlaqps (M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, VN2, AUXV, F, LDF)
           DLAQPS computes a step of QR factorization with column pivoting of a real m-by-n
           matrix A by using BLAS level 3.
       subroutine dlaqr0 (WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, WORK,
           LWORK, INFO)
           DLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices
           from the Schur decomposition.
       subroutine dlaqr1 (N, H, LDH, SR1, SI1, SR2, SI2, V)
           DLAQR1 sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3
           matrix H and specified shifts.
       subroutine dlaqr2 (WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND,
           SR, SI, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)
           DLAQR2 performs the orthogonal similarity transformation of a Hessenberg matrix to
           detect and deflate fully converged eigenvalues from a trailing principal submatrix
           (aggressive early deflation).
       subroutine dlaqr3 (WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND,
           SR, SI, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)
           DLAQR3 performs the orthogonal similarity transformation of a Hessenberg matrix to
           detect and deflate fully converged eigenvalues from a trailing principal submatrix
           (aggressive early deflation).
       subroutine dlaqr4 (WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, WORK,
           LWORK, INFO)
           DLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices
           from the Schur decomposition.
       subroutine dlaqr5 (WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, SR, SI, H, LDH, ILOZ,
           IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH)
           DLAQR5 performs a single small-bulge multi-shift QR sweep.
       subroutine dlaqsb (UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED)
           DLAQSB scales a symmetric/Hermitian band matrix, using scaling factors computed by
           spbequ.
       subroutine dlaqsp (UPLO, N, AP, S, SCOND, AMAX, EQUED)
           DLAQSP scales a symmetric/Hermitian matrix in packed storage, using scaling factors
           computed by sppequ.
       subroutine dlaqtr (LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK, INFO)
           DLAQTR solves a real quasi-triangular system of equations, or a complex quasi-
           triangular system of special form, in real arithmetic.
       subroutine dlar1v (N, B1, BN, LAMBDA, D, L, LD, LLD, PIVMIN, GAPTOL, Z, WANTNC, NEGCNT,
           ZTZ, MINGMA, R, ISUPPZ, NRMINV, RESID, RQCORR, WORK)
           DLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1
           through bn of the tridiagonal matrix LDLT - λI.
       subroutine dlar2v (N, X, Y, Z, INCX, C, S, INCC)
           DLAR2V applies a vector of plane rotations with real cosines and real sines from both
           sides to a sequence of 2-by-2 symmetric/Hermitian matrices.
       subroutine dlarf (SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
           DLARF applies an elementary reflector to a general rectangular matrix.
       subroutine dlarfb (SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK,
           LDWORK)
           DLARFB applies a block reflector or its transpose to a general rectangular matrix.
       subroutine dlarfg (N, ALPHA, X, INCX, TAU)
           DLARFG generates an elementary reflector (Householder matrix).
       subroutine dlarfgp (N, ALPHA, X, INCX, TAU)
           DLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
       subroutine dlarft (DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
           DLARFT forms the triangular factor T of a block reflector H = I - vtvH
       subroutine dlarfx (SIDE, M, N, V, TAU, C, LDC, WORK)
           DLARFX applies an elementary reflector to a general rectangular matrix, with loop
           unrolling when the reflector has order ≤ 10.
       subroutine dlargv (N, X, INCX, Y, INCY, C, INCC)
           DLARGV generates a vector of plane rotations with real cosines and real sines.
       subroutine dlarrv (N, VL, VU, D, L, PIVMIN, ISPLIT, M, DOL, DOU, MINRGP, RTOL1, RTOL2, W,
           WERR, WGAP, IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ, WORK, IWORK, INFO)
           DLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and
           the eigenvalues of L D LT.
       subroutine dlartv (N, X, INCX, Y, INCY, C, S, INCC)
           DLARTV applies a vector of plane rotations with real cosines and real sines to the
           elements of a pair of vectors.
       subroutine dlaswp (N, A, LDA, K1, K2, IPIV, INCX)
           DLASWP performs a series of row interchanges on a general rectangular matrix.
       subroutine dlat2s (UPLO, N, A, LDA, SA, LDSA, INFO)
           DLAT2S converts a double-precision triangular matrix to a single-precision triangular
           matrix.
       subroutine dlatbs (UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE, CNORM, INFO)
           DLATBS solves a triangular banded system of equations.
       subroutine dlatdf (IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, JPIV)
           DLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes
           a contribution to the reciprocal Dif-estimate.
       subroutine dlatps (UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE, CNORM, INFO)
           DLATPS solves a triangular system of equations with the matrix held in packed storage.
       subroutine dlatrd (UPLO, N, NB, A, LDA, E, TAU, W, LDW)
           DLATRD reduces the first nb rows and columns of a symmetric/Hermitian matrix A to real
           tridiagonal form by an orthogonal similarity transformation.
       subroutine dlatrs (UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
           DLATRS solves a triangular system of equations with the scale factor set to prevent
           overflow.
       subroutine dlauu2 (UPLO, N, A, LDA, INFO)
           DLAUU2 computes the product UUH or LHL, where U and L are upper or lower triangular
           matrices (unblocked algorithm).
       subroutine dlauum (UPLO, N, A, LDA, INFO)
           DLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular
           matrices (blocked algorithm).
       subroutine drscl (N, SA, SX, INCX)
           DRSCL multiplies a vector by the reciprocal of a real scalar.
       subroutine dtprfb (SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V, LDV, T, LDT, A, LDA, B,
           LDB, WORK, LDWORK)
           DTPRFB applies a real or complex 'triangular-pentagonal' blocked reflector to a real
           or complex matrix, which is composed of two blocks.
       subroutine slatrd (UPLO, N, NB, A, LDA, E, TAU, W, LDW)
           SLATRD reduces the first nb rows and columns of a symmetric/Hermitian matrix A to real
           tridiagonal form by an orthogonal similarity transformation.

Detailed Description

       This is the group of double other auxiliary routines

Function Documentation

   subroutine dlabrd (integer M, integer N, integer NB, double precision, dimension( lda, * ) A,
       integer LDA, double precision, dimension( * ) D, double precision, dimension( * ) E,
       double precision, dimension( * ) TAUQ, double precision, dimension( * ) TAUP, double
       precision, dimension( ldx, * ) X, integer LDX, double precision, dimension( ldy, * ) Y,
       integer LDY)
       DLABRD reduces the first nb rows and columns of a general matrix to a bidiagonal form.

       Purpose:

            DLABRD reduces the first NB rows and columns of a real general
            m by n matrix A to upper or lower bidiagonal form by an orthogonal
            transformation Q**T * A * P, and returns the matrices X and Y which
            are needed to apply the transformation to the unreduced part of A.

            If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower
            bidiagonal form.

            This is an auxiliary routine called by DGEBRD

       Parameters:
           M

                     M is INTEGER
                     The number of rows in the matrix A.

           N

                     N is INTEGER
                     The number of columns in the matrix A.

           NB

                     NB is INTEGER
                     The number of leading rows and columns of A to be reduced.

           A

                     A is DOUBLE PRECISION array, dimension (LDA,N)
                     On entry, the m by n general matrix to be reduced.
                     On exit, the first NB rows and columns of the matrix are
                     overwritten; the rest of the array is unchanged.
                     If m >= n, elements on and below the diagonal in the first NB
                       columns, with the array TAUQ, represent the orthogonal
                       matrix Q as a product of elementary reflectors; and
                       elements above the diagonal in the first NB rows, with the
                       array TAUP, represent the orthogonal matrix P as a product
                       of elementary reflectors.
                     If m < n, elements below the diagonal in the first NB
                       columns, with the array TAUQ, represent the orthogonal
                       matrix Q as a product of elementary reflectors, and
                       elements on and above the diagonal in the first NB rows,
                       with the array TAUP, represent the orthogonal matrix P as
                       a product of elementary reflectors.
                     See Further Details.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,M).

           D

                     D is DOUBLE PRECISION array, dimension (NB)
                     The diagonal elements of the first NB rows and columns of
                     the reduced matrix.  D(i) = A(i,i).

           E

                     E is DOUBLE PRECISION array, dimension (NB)
                     The off-diagonal elements of the first NB rows and columns of
                     the reduced matrix.

           TAUQ

                     TAUQ is DOUBLE PRECISION array, dimension (NB)
                     The scalar factors of the elementary reflectors which
                     represent the orthogonal matrix Q. See Further Details.

           TAUP

                     TAUP is DOUBLE PRECISION array, dimension (NB)
                     The scalar factors of the elementary reflectors which
                     represent the orthogonal matrix P. See Further Details.

           X

                     X is DOUBLE PRECISION array, dimension (LDX,NB)
                     The m-by-nb matrix X required to update the unreduced part
                     of A.

           LDX

                     LDX is INTEGER
                     The leading dimension of the array X. LDX >= max(1,M).

           Y

                     Y is DOUBLE PRECISION array, dimension (LDY,NB)
                     The n-by-nb matrix Y required to update the unreduced part
                     of A.

           LDY

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

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           June 2017

       Further Details:

             The matrices Q and P are represented as products of elementary
             reflectors:

                Q = H(1) H(2) . . . H(nb)  and  P = G(1) G(2) . . . G(nb)

             Each H(i) and G(i) has the form:

                H(i) = I - tauq * v * v**T  and G(i) = I - taup * u * u**T

             where tauq and taup are real scalars, and v and u are real vectors.

             If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in
             A(i:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is stored on exit in
             A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).

             If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in
             A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in
             A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).

             The elements of the vectors v and u together form the m-by-nb matrix
             V and the nb-by-n matrix U**T which are needed, with X and Y, to apply
             the transformation to the unreduced part of the matrix, using a block
             update of the form:  A := A - V*Y**T - X*U**T.

             The contents of A on exit are illustrated by the following examples
             with nb = 2:

             m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):

               (  1   1   u1  u1  u1 )           (  1   u1  u1  u1  u1  u1 )
               (  v1  1   1   u2  u2 )           (  1   1   u2  u2  u2  u2 )
               (  v1  v2  a   a   a  )           (  v1  1   a   a   a   a  )
               (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )
               (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )
               (  v1  v2  a   a   a  )

             where a denotes an element of the original matrix which is unchanged,
             vi denotes an element of the vector defining H(i), and ui an element
             of the vector defining G(i).

   subroutine dlacn2 (integer N, double precision, dimension( * ) V, double precision, dimension(
       * ) X, integer, dimension( * ) ISGN, double precision EST, integer KASE, integer,
       dimension( 3 ) ISAVE)
       DLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating
       matrix-vector products.

       Purpose:

            DLACN2 estimates the 1-norm of a square, real matrix A.
            Reverse communication is used for evaluating matrix-vector products.

       Parameters:
           N

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

           V

                     V is DOUBLE PRECISION array, dimension (N)
                    On the final return, V = A*W,  where  EST = norm(V)/norm(W)
                    (W is not returned).

           X

                     X is DOUBLE PRECISION array, dimension (N)
                    On an intermediate return, X should be overwritten by
                          A * X,   if KASE=1,
                          A**T * X,  if KASE=2,
                    and DLACN2 must be re-called with all the other parameters
                    unchanged.

           ISGN

                     ISGN is INTEGER array, dimension (N)

           EST

                     EST is DOUBLE PRECISION
                    On entry with KASE = 1 or 2 and ISAVE(1) = 3, EST should be
                    unchanged from the previous call to DLACN2.
                    On exit, EST is an estimate (a lower bound) for norm(A).

           KASE

                     KASE is INTEGER
                    On the initial call to DLACN2, KASE should be 0.
                    On an intermediate return, KASE will be 1 or 2, indicating
                    whether X should be overwritten by A * X  or A**T * X.
                    On the final return from DLACN2, KASE will again be 0.

           ISAVE

                     ISAVE is INTEGER array, dimension (3)
                    ISAVE is used to save variables between calls to DLACN2

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

       Further Details:

             Originally named SONEST, dated March 16, 1988.

             This is a thread safe version of DLACON, which uses the array ISAVE
             in place of a SAVE statement, as follows:

                DLACON     DLACN2
                 JUMP     ISAVE(1)
                 J        ISAVE(2)
                 ITER     ISAVE(3)

       Contributors:
           Nick Higham, University of Manchester

       References:
           N.J. Higham, 'FORTRAN codes for estimating the one-norm of
             a real or complex matrix, with applications to condition estimation', ACM Trans.
           Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.

   subroutine dlacon (integer N, double precision, dimension( * ) V, double precision, dimension(
       * ) X, integer, dimension( * ) ISGN, double precision EST, integer KASE)
       DLACON estimates the 1-norm of a square matrix, using reverse communication for evaluating
       matrix-vector products.

       Purpose:

            DLACON estimates the 1-norm of a square, real matrix A.
            Reverse communication is used for evaluating matrix-vector products.

       Parameters:
           N

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

           V

                     V is DOUBLE PRECISION array, dimension (N)
                    On the final return, V = A*W,  where  EST = norm(V)/norm(W)
                    (W is not returned).

           X

                     X is DOUBLE PRECISION array, dimension (N)
                    On an intermediate return, X should be overwritten by
                          A * X,   if KASE=1,
                          A**T * X,  if KASE=2,
                    and DLACON must be re-called with all the other parameters
                    unchanged.

           ISGN

                     ISGN is INTEGER array, dimension (N)

           EST

                     EST is DOUBLE PRECISION
                    On entry with KASE = 1 or 2 and JUMP = 3, EST should be
                    unchanged from the previous call to DLACON.
                    On exit, EST is an estimate (a lower bound) for norm(A).

           KASE

                     KASE is INTEGER
                    On the initial call to DLACON, KASE should be 0.
                    On an intermediate return, KASE will be 1 or 2, indicating
                    whether X should be overwritten by A * X  or A**T * X.
                    On the final return from DLACON, KASE will again be 0.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

       Contributors:
           Nick Higham, University of Manchester.
            Originally named SONEST, dated March 16, 1988.

       References:
           N.J. Higham, 'FORTRAN codes for estimating the one-norm of
             a real or complex matrix, with applications to condition estimation', ACM Trans.
           Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.

   subroutine dladiv (double precision A, double precision B, double precision C, double
       precision D, double precision P, double precision Q)
       DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.

       Purpose:

            DLADIV performs complex division in  real arithmetic

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

            The algorithm is due to Michael Baudin and Robert L. Smith
            and can be found in the paper
            "A Robust Complex Division in Scilab"

       Parameters:
           A

                     A is DOUBLE PRECISION

           B

                     B is DOUBLE PRECISION

           C

                     C is DOUBLE PRECISION

           D

                     D is DOUBLE PRECISION
                     The scalars a, b, c, and d in the above expression.

           P

                     P is DOUBLE PRECISION

           Q

                     Q is DOUBLE PRECISION
                     The scalars p and q in the above expression.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           January 2013

   subroutine dlaein (logical RIGHTV, logical NOINIT, integer N, double precision, dimension(
       ldh, * ) H, integer LDH, double precision WR, double precision WI, double precision,
       dimension( * ) VR, double precision, dimension( * ) VI, double precision, dimension( ldb,
       * ) B, integer LDB, double precision, dimension( * ) WORK, double precision EPS3, double
       precision SMLNUM, double precision BIGNUM, integer INFO)
       DLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by
       inverse iteration.

       Purpose:

            DLAEIN uses inverse iteration to find a right or left eigenvector
            corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg
            matrix H.

       Parameters:
           RIGHTV

                     RIGHTV is LOGICAL
                     = .TRUE. : compute right eigenvector;
                     = .FALSE.: compute left eigenvector.

           NOINIT

                     NOINIT is LOGICAL
                     = .TRUE. : no initial vector supplied in (VR,VI).
                     = .FALSE.: initial vector supplied in (VR,VI).

           N

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

           H

                     H is DOUBLE PRECISION array, dimension (LDH,N)
                     The upper Hessenberg matrix H.

           LDH

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

           WR

                     WR is DOUBLE PRECISION

           WI

                     WI is DOUBLE PRECISION
                     The real and imaginary parts of the eigenvalue of H whose
                     corresponding right or left eigenvector is to be computed.

           VR

                     VR is DOUBLE PRECISION array, dimension (N)

           VI

                     VI is DOUBLE PRECISION array, dimension (N)
                     On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain
                     a real starting vector for inverse iteration using the real
                     eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI
                     must contain the real and imaginary parts of a complex
                     starting vector for inverse iteration using the complex
                     eigenvalue (WR,WI); otherwise VR and VI need not be set.
                     On exit, if WI = 0.0 (real eigenvalue), VR contains the
                     computed real eigenvector; if WI.ne.0.0 (complex eigenvalue),
                     VR and VI contain the real and imaginary parts of the
                     computed complex eigenvector. The eigenvector is normalized
                     so that the component of largest magnitude has magnitude 1;
                     here the magnitude of a complex number (x,y) is taken to be
                     |x| + |y|.
                     VI is not referenced if WI = 0.0.

           B

                     B is DOUBLE PRECISION array, dimension (LDB,N)

           LDB

                     LDB is INTEGER
                     The leading dimension of the array B.  LDB >= N+1.

           WORK

                     WORK is DOUBLE PRECISION array, dimension (N)

           EPS3

                     EPS3 is DOUBLE PRECISION
                     A small machine-dependent value which is used to perturb
                     close eigenvalues, and to replace zero pivots.

           SMLNUM

                     SMLNUM is DOUBLE PRECISION
                     A machine-dependent value close to the underflow threshold.

           BIGNUM

                     BIGNUM is DOUBLE PRECISION
                     A machine-dependent value close to the overflow threshold.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     = 1:  inverse iteration did not converge; VR is set to the
                           last iterate, and so is VI if WI.ne.0.0.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlaexc (logical WANTQ, integer N, double precision, dimension( ldt, * ) T, integer
       LDT, double precision, dimension( ldq, * ) Q, integer LDQ, integer J1, integer N1, integer
       N2, double precision, dimension( * ) WORK, integer INFO)
       DLAEXC swaps adjacent diagonal blocks of a real upper quasi-triangular matrix in Schur
       canonical form, by an orthogonal similarity transformation.

       Purpose:

            DLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
            an upper quasi-triangular matrix T by an orthogonal similarity
            transformation.

            T must be in Schur canonical form, that is, block upper triangular
            with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block
            has its diagonal elemnts equal and its off-diagonal elements of
            opposite sign.

       Parameters:
           WANTQ

                     WANTQ is LOGICAL
                     = .TRUE. : accumulate the transformation in the matrix Q;
                     = .FALSE.: do not accumulate the transformation.

           N

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

           T

                     T is DOUBLE PRECISION array, dimension (LDT,N)
                     On entry, the upper quasi-triangular matrix T, in Schur
                     canonical form.
                     On exit, the updated matrix T, again in Schur canonical form.

           LDT

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

           Q

                     Q is DOUBLE PRECISION array, dimension (LDQ,N)
                     On entry, if WANTQ is .TRUE., the orthogonal matrix Q.
                     On exit, if WANTQ is .TRUE., the updated matrix Q.
                     If WANTQ is .FALSE., Q is not referenced.

           LDQ

                     LDQ is INTEGER
                     The leading dimension of the array Q.
                     LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.

           J1

                     J1 is INTEGER
                     The index of the first row of the first block T11.

           N1

                     N1 is INTEGER
                     The order of the first block T11. N1 = 0, 1 or 2.

           N2

                     N2 is INTEGER
                     The order of the second block T22. N2 = 0, 1 or 2.

           WORK

                     WORK is DOUBLE PRECISION array, dimension (N)

           INFO

                     INFO is INTEGER
                     = 0: successful exit
                     = 1: the transformed matrix T would be too far from Schur
                          form; the blocks are not swapped and T and Q are
                          unchanged.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlag2 (double precision, dimension( lda, * ) A, integer LDA, double precision,
       dimension( ldb, * ) B, integer LDB, double precision SAFMIN, double precision SCALE1,
       double precision SCALE2, double precision WR1, double precision WR2, double precision WI)
       DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as
       necessary to avoid over-/underflow.

       Purpose:

            DLAG2 computes the eigenvalues of a 2 x 2 generalized eigenvalue
            problem  A - w B, with scaling as necessary to avoid over-/underflow.

            The scaling factor "s" results in a modified eigenvalue equation

                s A - w B

            where  s  is a non-negative scaling factor chosen so that  w,  w B,
            and  s A  do not overflow and, if possible, do not underflow, either.

       Parameters:
           A

                     A is DOUBLE PRECISION array, dimension (LDA, 2)
                     On entry, the 2 x 2 matrix A.  It is assumed that its 1-norm
                     is less than 1/SAFMIN.  Entries less than
                     sqrt(SAFMIN)*norm(A) are subject to being treated as zero.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= 2.

           B

                     B is DOUBLE PRECISION array, dimension (LDB, 2)
                     On entry, the 2 x 2 upper triangular matrix B.  It is
                     assumed that the one-norm of B is less than 1/SAFMIN.  The
                     diagonals should be at least sqrt(SAFMIN) times the largest
                     element of B (in absolute value); if a diagonal is smaller
                     than that, then  +/- sqrt(SAFMIN) will be used instead of
                     that diagonal.

           LDB

                     LDB is INTEGER
                     The leading dimension of the array B.  LDB >= 2.

           SAFMIN

                     SAFMIN is DOUBLE PRECISION
                     The smallest positive number s.t. 1/SAFMIN does not
                     overflow.  (This should always be DLAMCH('S') -- it is an
                     argument in order to avoid having to call DLAMCH frequently.)

           SCALE1

                     SCALE1 is DOUBLE PRECISION
                     A scaling factor used to avoid over-/underflow in the
                     eigenvalue equation which defines the first eigenvalue.  If
                     the eigenvalues are complex, then the eigenvalues are
                     ( WR1  +/-  WI i ) / SCALE1  (which may lie outside the
                     exponent range of the machine), SCALE1=SCALE2, and SCALE1
                     will always be positive.  If the eigenvalues are real, then
                     the first (real) eigenvalue is  WR1 / SCALE1 , but this may
                     overflow or underflow, and in fact, SCALE1 may be zero or
                     less than the underflow threshold if the exact eigenvalue
                     is sufficiently large.

           SCALE2

                     SCALE2 is DOUBLE PRECISION
                     A scaling factor used to avoid over-/underflow in the
                     eigenvalue equation which defines the second eigenvalue.  If
                     the eigenvalues are complex, then SCALE2=SCALE1.  If the
                     eigenvalues are real, then the second (real) eigenvalue is
                     WR2 / SCALE2 , but this may overflow or underflow, and in
                     fact, SCALE2 may be zero or less than the underflow
                     threshold if the exact eigenvalue is sufficiently large.

           WR1

                     WR1 is DOUBLE PRECISION
                     If the eigenvalue is real, then WR1 is SCALE1 times the
                     eigenvalue closest to the (2,2) element of A B**(-1).  If the
                     eigenvalue is complex, then WR1=WR2 is SCALE1 times the real
                     part of the eigenvalues.

           WR2

                     WR2 is DOUBLE PRECISION
                     If the eigenvalue is real, then WR2 is SCALE2 times the
                     other eigenvalue.  If the eigenvalue is complex, then
                     WR1=WR2 is SCALE1 times the real part of the eigenvalues.

           WI

                     WI is DOUBLE PRECISION
                     If the eigenvalue is real, then WI is zero.  If the
                     eigenvalue is complex, then WI is SCALE1 times the imaginary
                     part of the eigenvalues.  WI will always be non-negative.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           June 2016

   subroutine dlag2s (integer M, integer N, double precision, dimension( lda, * ) A, integer LDA,
       real, dimension( ldsa, * ) SA, integer LDSA, integer INFO)
       DLAG2S converts a double precision matrix to a single precision matrix.

       Purpose:

            DLAG2S converts a DOUBLE PRECISION matrix, SA, to a SINGLE
            PRECISION matrix, A.

            RMAX is the overflow for the SINGLE PRECISION arithmetic
            DLAG2S checks that all the entries of A are between -RMAX and
            RMAX. If not the conversion is aborted and a flag is raised.

            This is an auxiliary routine so there is no argument checking.

       Parameters:
           M

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

           N

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

           A

                     A is DOUBLE PRECISION array, dimension (LDA,N)
                     On entry, the M-by-N coefficient matrix A.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,M).

           SA

                     SA is REAL array, dimension (LDSA,N)
                     On exit, if INFO=0, the M-by-N coefficient matrix SA; if
                     INFO>0, the content of SA is unspecified.

           LDSA

                     LDSA is INTEGER
                     The leading dimension of the array SA.  LDSA >= max(1,M).

           INFO

                     INFO is INTEGER
                     = 0:  successful exit.
                     = 1:  an entry of the matrix A is greater than the SINGLE
                           PRECISION overflow threshold, in this case, the content
                           of SA in exit is unspecified.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlags2 (logical UPPER, double precision A1, double precision A2, double precision
       A3, double precision B1, double precision B2, double precision B3, double precision CSU,
       double precision SNU, double precision CSV, double precision SNV, double precision CSQ,
       double precision SNQ)
       DLAGS2 computes 2-by-2 orthogonal matrices U, V, and Q, and applies them to matrices A and
       B such that the rows of the transformed A and B are parallel.

       Purpose:

            DLAGS2 computes 2-by-2 orthogonal matrices U, V and Q, such
            that if ( UPPER ) then

                      U**T *A*Q = U**T *( A1 A2 )*Q = ( x  0  )
                                        ( 0  A3 )     ( x  x  )
            and
                      V**T*B*Q = V**T *( B1 B2 )*Q = ( x  0  )
                                       ( 0  B3 )     ( x  x  )

            or if ( .NOT.UPPER ) then

                      U**T *A*Q = U**T *( A1 0  )*Q = ( x  x  )
                                        ( A2 A3 )     ( 0  x  )
            and
                      V**T*B*Q = V**T*( B1 0  )*Q = ( x  x  )
                                      ( B2 B3 )     ( 0  x  )

            The rows of the transformed A and B are parallel, where

              U = (  CSU  SNU ), V = (  CSV SNV ), Q = (  CSQ   SNQ )
                  ( -SNU  CSU )      ( -SNV CSV )      ( -SNQ   CSQ )

            Z**T denotes the transpose of Z.

       Parameters:
           UPPER

                     UPPER is LOGICAL
                     = .TRUE.: the input matrices A and B are upper triangular.
                     = .FALSE.: the input matrices A and B are lower triangular.

           A1

                     A1 is DOUBLE PRECISION

           A2

                     A2 is DOUBLE PRECISION

           A3

                     A3 is DOUBLE PRECISION
                     On entry, A1, A2 and A3 are elements of the input 2-by-2
                     upper (lower) triangular matrix A.

           B1

                     B1 is DOUBLE PRECISION

           B2

                     B2 is DOUBLE PRECISION

           B3

                     B3 is DOUBLE PRECISION
                     On entry, B1, B2 and B3 are elements of the input 2-by-2
                     upper (lower) triangular matrix B.

           CSU

                     CSU is DOUBLE PRECISION

           SNU

                     SNU is DOUBLE PRECISION
                     The desired orthogonal matrix U.

           CSV

                     CSV is DOUBLE PRECISION

           SNV

                     SNV is DOUBLE PRECISION
                     The desired orthogonal matrix V.

           CSQ

                     CSQ is DOUBLE PRECISION

           SNQ

                     SNQ is DOUBLE PRECISION
                     The desired orthogonal matrix Q.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlagtm (character TRANS, integer N, integer NRHS, double precision ALPHA, double
       precision, dimension( * ) DL, double precision, dimension( * ) D, double precision,
       dimension( * ) DU, double precision, dimension( ldx, * ) X, integer LDX, double precision
       BETA, double precision, dimension( ldb, * ) B, integer LDB)
       DLAGTM performs a matrix-matrix product of the form C = αAB+βC, where A is a tridiagonal
       matrix, B and C are rectangular matrices, and α and β are scalars, which may be 0, 1, or
       -1.

       Purpose:

            DLAGTM performs a matrix-vector product of the form

               B := alpha * A * X + beta * B

            where A is a tridiagonal matrix of order N, B and X are N by NRHS
            matrices, and alpha and beta are real scalars, each of which may be
            0., 1., or -1.

       Parameters:
           TRANS

                     TRANS is CHARACTER*1
                     Specifies the operation applied to A.
                     = 'N':  No transpose, B := alpha * A * X + beta * B
                     = 'T':  Transpose,    B := alpha * A'* X + beta * B
                     = 'C':  Conjugate transpose = Transpose

           N

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

           NRHS

                     NRHS is INTEGER
                     The number of right hand sides, i.e., the number of columns
                     of the matrices X and B.

           ALPHA

                     ALPHA is DOUBLE PRECISION
                     The scalar alpha.  ALPHA must be 0., 1., or -1.; otherwise,
                     it is assumed to be 0.

           DL

                     DL is DOUBLE PRECISION array, dimension (N-1)
                     The (n-1) sub-diagonal elements of T.

           D

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

           DU

                     DU is DOUBLE PRECISION array, dimension (N-1)
                     The (n-1) super-diagonal elements of T.

           X

                     X is DOUBLE PRECISION array, dimension (LDX,NRHS)
                     The N by NRHS matrix X.

           LDX

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

           BETA

                     BETA is DOUBLE PRECISION
                     The scalar beta.  BETA must be 0., 1., or -1.; otherwise,
                     it is assumed to be 1.

           B

                     B is DOUBLE PRECISION array, dimension (LDB,NRHS)
                     On entry, the N by NRHS matrix B.
                     On exit, B is overwritten by the matrix expression
                     B := alpha * A * X + beta * B.

           LDB

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

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlagv2 (double precision, dimension( lda, * ) A, integer LDA, double precision,
       dimension( ldb, * ) B, integer LDB, double precision, dimension( 2 ) ALPHAR, double
       precision, dimension( 2 ) ALPHAI, double precision, dimension( 2 ) BETA, double precision
       CSL, double precision SNL, double precision CSR, double precision SNR)
       DLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A,B)
       where B is upper triangular.

       Purpose:

            DLAGV2 computes the Generalized Schur factorization of a real 2-by-2
            matrix pencil (A,B) where B is upper triangular. This routine
            computes orthogonal (rotation) matrices given by CSL, SNL and CSR,
            SNR such that

            1) if the pencil (A,B) has two real eigenvalues (include 0/0 or 1/0
               types), then

               [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
               [  0  a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]

               [ b11 b12 ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
               [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ],

            2) if the pencil (A,B) has a pair of complex conjugate eigenvalues,
               then

               [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
               [ a21 a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]

               [ b11  0  ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
               [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ]

               where b11 >= b22 > 0.

       Parameters:
           A

                     A is DOUBLE PRECISION array, dimension (LDA, 2)
                     On entry, the 2 x 2 matrix A.
                     On exit, A is overwritten by the ``A-part'' of the
                     generalized Schur form.

           LDA

                     LDA is INTEGER
                     THe leading dimension of the array A.  LDA >= 2.

           B

                     B is DOUBLE PRECISION array, dimension (LDB, 2)
                     On entry, the upper triangular 2 x 2 matrix B.
                     On exit, B is overwritten by the ``B-part'' of the
                     generalized Schur form.

           LDB

                     LDB is INTEGER
                     THe leading dimension of the array B.  LDB >= 2.

           ALPHAR

                     ALPHAR is DOUBLE PRECISION array, dimension (2)

           ALPHAI

                     ALPHAI is DOUBLE PRECISION array, dimension (2)

           BETA

                     BETA is DOUBLE PRECISION array, dimension (2)
                     (ALPHAR(k)+i*ALPHAI(k))/BETA(k) are the eigenvalues of the
                     pencil (A,B), k=1,2, i = sqrt(-1).  Note that BETA(k) may
                     be zero.

           CSL

                     CSL is DOUBLE PRECISION
                     The cosine of the left rotation matrix.

           SNL

                     SNL is DOUBLE PRECISION
                     The sine of the left rotation matrix.

           CSR

                     CSR is DOUBLE PRECISION
                     The cosine of the right rotation matrix.

           SNR

                     SNR is DOUBLE PRECISION
                     The sine of the right rotation matrix.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

       Contributors:
           Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA

   subroutine dlahqr (logical WANTT, logical WANTZ, integer N, integer ILO, integer IHI, double
       precision, dimension( ldh, * ) H, integer LDH, double precision, dimension( * ) WR, double
       precision, dimension( * ) WI, integer ILOZ, integer IHIZ, double precision, dimension(
       ldz, * ) Z, integer LDZ, integer INFO)
       DLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,
       using the double-shift/single-shift QR algorithm.

       Purpose:

               DLAHQR is an auxiliary routine called by DHSEQR to update the
               eigenvalues and Schur decomposition already computed by DHSEQR, by
               dealing with the Hessenberg submatrix in rows and columns ILO to
               IHI.

       Parameters:
           WANTT

                     WANTT is LOGICAL
                     = .TRUE. : the full Schur form T is required;
                     = .FALSE.: only eigenvalues are required.

           WANTZ

                     WANTZ is LOGICAL
                     = .TRUE. : the matrix of Schur vectors Z is required;
                     = .FALSE.: Schur vectors are not required.

           N

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

           ILO

                     ILO is INTEGER

           IHI

                     IHI is INTEGER
                     It is assumed that H is already upper quasi-triangular in
                     rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
                     ILO = 1). DLAHQR works primarily with the Hessenberg
                     submatrix in rows and columns ILO to IHI, but applies
                     transformations to all of H if WANTT is .TRUE..
                     1 <= ILO <= max(1,IHI); IHI <= N.

           H

                     H is DOUBLE PRECISION array, dimension (LDH,N)
                     On entry, the upper Hessenberg matrix H.
                     On exit, if INFO is zero and if WANTT is .TRUE., H is upper
                     quasi-triangular in rows and columns ILO:IHI, with any
                     2-by-2 diagonal blocks in standard form. If INFO is zero
                     and WANTT is .FALSE., the contents of H are unspecified on
                     exit.  The output state of H if INFO is nonzero is given
                     below under the description of INFO.

           LDH

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

           WR

                     WR is DOUBLE PRECISION array, dimension (N)

           WI

                     WI is DOUBLE PRECISION array, dimension (N)
                     The real and imaginary parts, respectively, of the computed
                     eigenvalues ILO to IHI are stored in the corresponding
                     elements of WR and WI. If two eigenvalues are computed as a
                     complex conjugate pair, they are stored in consecutive
                     elements of WR and WI, say the i-th and (i+1)th, with
                     WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
                     eigenvalues are stored in the same order as on the diagonal
                     of the Schur form returned in H, with WR(i) = H(i,i), and, if
                     H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
                     WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).

           ILOZ

                     ILOZ is INTEGER

           IHIZ

                     IHIZ is INTEGER
                     Specify the rows of Z to which transformations must be
                     applied if WANTZ is .TRUE..
                     1 <= ILOZ <= ILO; IHI <= IHIZ <= N.

           Z

                     Z is DOUBLE PRECISION array, dimension (LDZ,N)
                     If WANTZ is .TRUE., on entry Z must contain the current
                     matrix Z of transformations accumulated by DHSEQR, and on
                     exit Z has been updated; transformations are applied only to
                     the submatrix Z(ILOZ:IHIZ,ILO:IHI).
                     If WANTZ is .FALSE., Z is not referenced.

           LDZ

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

           INFO

                     INFO is INTEGER
                      =   0: successful exit
                     .GT. 0: If INFO = i, DLAHQR failed to compute all the
                             eigenvalues ILO to IHI in a total of 30 iterations
                             per eigenvalue; elements i+1:ihi of WR and WI
                             contain those eigenvalues which have been
                             successfully computed.

                             If INFO .GT. 0 and WANTT is .FALSE., then on exit,
                             the remaining unconverged eigenvalues are the
                             eigenvalues of the upper Hessenberg matrix rows
                             and columns ILO thorugh INFO of the final, output
                             value of H.

                             If INFO .GT. 0 and WANTT is .TRUE., then on exit
                     (*)       (initial value of H)*U  = U*(final value of H)
                             where U is an orthognal matrix.    The final
                             value of H is upper Hessenberg and triangular in
                             rows and columns INFO+1 through IHI.

                             If INFO .GT. 0 and WANTZ is .TRUE., then on exit
                                 (final value of Z)  = (initial value of Z)*U
                             where U is the orthogonal matrix in (*)
                             (regardless of the value of WANTT.)

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

       Further Details:

                02-96 Based on modifications by
                David Day, Sandia National Laboratory, USA

                12-04 Further modifications by
                Ralph Byers, University of Kansas, USA
                This is a modified version of DLAHQR from LAPACK version 3.0.
                It is (1) more robust against overflow and underflow and
                (2) adopts the more conservative Ahues & Tisseur stopping
                criterion (LAWN 122, 1997).

   subroutine dlahr2 (integer N, integer K, integer NB, double precision, dimension( lda, * ) A,
       integer LDA, double precision, dimension( nb ) TAU, double precision, dimension( ldt, nb )
       T, integer LDT, double precision, dimension( ldy, nb ) Y, integer LDY)
       DLAHR2 reduces the specified number of first columns of a general rectangular matrix A so
       that elements below the specified subdiagonal are zero, and returns auxiliary matrices
       which are needed to apply the transformation to the unreduced part of A.

       Purpose:

            DLAHR2 reduces the first NB columns of A real general n-BY-(n-k+1)
            matrix A so that elements below the k-th subdiagonal are zero. The
            reduction is performed by an orthogonal similarity transformation
            Q**T * A * Q. The routine returns the matrices V and T which determine
            Q as a block reflector I - V*T*V**T, and also the matrix Y = A * V * T.

            This is an auxiliary routine called by DGEHRD.

       Parameters:
           N

                     N is INTEGER
                     The order of the matrix A.

           K

                     K is INTEGER
                     The offset for the reduction. Elements below the k-th
                     subdiagonal in the first NB columns are reduced to zero.
                     K < N.

           NB

                     NB is INTEGER
                     The number of columns to be reduced.

           A

                     A is DOUBLE PRECISION array, dimension (LDA,N-K+1)
                     On entry, the n-by-(n-k+1) general matrix A.
                     On exit, the elements on and above the k-th subdiagonal in
                     the first NB columns are overwritten with the corresponding
                     elements of the reduced matrix; the elements below the k-th
                     subdiagonal, with the array TAU, represent the matrix Q as a
                     product of elementary reflectors. The other columns of A are
                     unchanged. See Further Details.

           LDA

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

           TAU

                     TAU is DOUBLE PRECISION array, dimension (NB)
                     The scalar factors of the elementary reflectors. See Further
                     Details.

           T

                     T is DOUBLE PRECISION array, dimension (LDT,NB)
                     The upper triangular matrix T.

           LDT

                     LDT is INTEGER
                     The leading dimension of the array T.  LDT >= NB.

           Y

                     Y is DOUBLE PRECISION array, dimension (LDY,NB)
                     The n-by-nb matrix Y.

           LDY

                     LDY is INTEGER
                     The leading dimension of the array Y. LDY >= N.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

       Further Details:

             The matrix Q is represented as a product of nb elementary reflectors

                Q = H(1) H(2) . . . H(nb).

             Each H(i) has the form

                H(i) = I - tau * v * v**T

             where tau is a real scalar, and v is a real vector with
             v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
             A(i+k+1:n,i), and tau in TAU(i).

             The elements of the vectors v together form the (n-k+1)-by-nb matrix
             V which is needed, with T and Y, to apply the transformation to the
             unreduced part of the matrix, using an update of the form:
             A := (I - V*T*V**T) * (A - Y*V**T).

             The contents of A on exit are illustrated by the following example
             with n = 7, k = 3 and nb = 2:

                ( a   a   a   a   a )
                ( a   a   a   a   a )
                ( a   a   a   a   a )
                ( h   h   a   a   a )
                ( v1  h   a   a   a )
                ( v1  v2  a   a   a )
                ( v1  v2  a   a   a )

             where a denotes an element of the original matrix A, h denotes a
             modified element of the upper Hessenberg matrix H, and vi denotes an
             element of the vector defining H(i).

             This subroutine is a slight modification of LAPACK-3.0's DLAHRD
             incorporating improvements proposed by Quintana-Orti and Van de
             Gejin. Note that the entries of A(1:K,2:NB) differ from those
             returned by the original LAPACK-3.0's DLAHRD routine. (This
             subroutine is not backward compatible with LAPACK-3.0's DLAHRD.)

       References:
           Gregorio Quintana-Orti and Robert van de Geijn, 'Improving the
             performance of reduction to Hessenberg form,' ACM Transactions on Mathematical
           Software, 32(2):180-194, June 2006.

   subroutine dlaic1 (integer JOB, integer J, double precision, dimension( j ) X, double
       precision SEST, double precision, dimension( j ) W, double precision GAMMA, double
       precision SESTPR, double precision S, double precision C)
       DLAIC1 applies one step of incremental condition estimation.

       Purpose:

            DLAIC1 applies one step of incremental condition estimation in
            its simplest version:

            Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
            lower triangular matrix L, such that
                     twonorm(L*x) = sest
            Then DLAIC1 computes sestpr, s, c such that
            the vector
                            [ s*x ]
                     xhat = [  c  ]
            is an approximate singular vector of
                            [ L       0  ]
                     Lhat = [ w**T gamma ]
            in the sense that
                     twonorm(Lhat*xhat) = sestpr.

            Depending on JOB, an estimate for the largest or smallest singular
            value is computed.

            Note that [s c]**T and sestpr**2 is an eigenpair of the system

                diag(sest*sest, 0) + [alpha  gamma] * [ alpha ]
                                                      [ gamma ]

            where  alpha =  x**T*w.

       Parameters:
           JOB

                     JOB is INTEGER
                     = 1: an estimate for the largest singular value is computed.
                     = 2: an estimate for the smallest singular value is computed.

           J

                     J is INTEGER
                     Length of X and W

           X

                     X is DOUBLE PRECISION array, dimension (J)
                     The j-vector x.

           SEST

                     SEST is DOUBLE PRECISION
                     Estimated singular value of j by j matrix L

           W

                     W is DOUBLE PRECISION array, dimension (J)
                     The j-vector w.

           GAMMA

                     GAMMA is DOUBLE PRECISION
                     The diagonal element gamma.

           SESTPR

                     SESTPR is DOUBLE PRECISION
                     Estimated singular value of (j+1) by (j+1) matrix Lhat.

           S

                     S is DOUBLE PRECISION
                     Sine needed in forming xhat.

           C

                     C is DOUBLE PRECISION
                     Cosine needed in forming xhat.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlaln2 (logical LTRANS, integer NA, integer NW, double precision SMIN, double
       precision CA, double precision, dimension( lda, * ) A, integer LDA, double precision D1,
       double precision D2, double precision, dimension( ldb, * ) B, integer LDB, double
       precision WR, double precision WI, double precision, dimension( ldx, * ) X, integer LDX,
       double precision SCALE, double precision XNORM, integer INFO)
       DLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.

       Purpose:

            DLALN2 solves a system of the form  (ca A - w D ) X = s B
            or (ca A**T - w D) X = s B   with possible scaling ("s") and
            perturbation of A.  (A**T means A-transpose.)

            A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA
            real diagonal matrix, w is a real or complex value, and X and B are
            NA x 1 matrices -- real if w is real, complex if w is complex.  NA
            may be 1 or 2.

            If w is complex, X and B are represented as NA x 2 matrices,
            the first column of each being the real part and the second
            being the imaginary part.

            "s" is a scaling factor (.LE. 1), computed by DLALN2, which is
            so chosen that X can be computed without overflow.  X is further
            scaled if necessary to assure that norm(ca A - w D)*norm(X) is less
            than overflow.

            If both singular values of (ca A - w D) are less than SMIN,
            SMIN*identity will be used instead of (ca A - w D).  If only one
            singular value is less than SMIN, one element of (ca A - w D) will be
            perturbed enough to make the smallest singular value roughly SMIN.
            If both singular values are at least SMIN, (ca A - w D) will not be
            perturbed.  In any case, the perturbation will be at most some small
            multiple of max( SMIN, ulp*norm(ca A - w D) ).  The singular values
            are computed by infinity-norm approximations, and thus will only be
            correct to a factor of 2 or so.

            Note: all input quantities are assumed to be smaller than overflow
            by a reasonable factor.  (See BIGNUM.)

       Parameters:
           LTRANS

                     LTRANS is LOGICAL
                     =.TRUE.:  A-transpose will be used.
                     =.FALSE.: A will be used (not transposed.)

           NA

                     NA is INTEGER
                     The size of the matrix A.  It may (only) be 1 or 2.

           NW

                     NW is INTEGER
                     1 if "w" is real, 2 if "w" is complex.  It may only be 1
                     or 2.

           SMIN

                     SMIN is DOUBLE PRECISION
                     The desired lower bound on the singular values of A.  This
                     should be a safe distance away from underflow or overflow,
                     say, between (underflow/machine precision) and  (machine
                     precision * overflow ).  (See BIGNUM and ULP.)

           CA

                     CA is DOUBLE PRECISION
                     The coefficient c, which A is multiplied by.

           A

                     A is DOUBLE PRECISION array, dimension (LDA,NA)
                     The NA x NA matrix A.

           LDA

                     LDA is INTEGER
                     The leading dimension of A.  It must be at least NA.

           D1

                     D1 is DOUBLE PRECISION
                     The 1,1 element in the diagonal matrix D.

           D2

                     D2 is DOUBLE PRECISION
                     The 2,2 element in the diagonal matrix D.  Not used if NA=1.

           B

                     B is DOUBLE PRECISION array, dimension (LDB,NW)
                     The NA x NW matrix B (right-hand side).  If NW=2 ("w" is
                     complex), column 1 contains the real part of B and column 2
                     contains the imaginary part.

           LDB

                     LDB is INTEGER
                     The leading dimension of B.  It must be at least NA.

           WR

                     WR is DOUBLE PRECISION
                     The real part of the scalar "w".

           WI

                     WI is DOUBLE PRECISION
                     The imaginary part of the scalar "w".  Not used if NW=1.

           X

                     X is DOUBLE PRECISION array, dimension (LDX,NW)
                     The NA x NW matrix X (unknowns), as computed by DLALN2.
                     If NW=2 ("w" is complex), on exit, column 1 will contain
                     the real part of X and column 2 will contain the imaginary
                     part.

           LDX

                     LDX is INTEGER
                     The leading dimension of X.  It must be at least NA.

           SCALE

                     SCALE is DOUBLE PRECISION
                     The scale factor that B must be multiplied by to insure
                     that overflow does not occur when computing X.  Thus,
                     (ca A - w D) X  will be SCALE*B, not B (ignoring
                     perturbations of A.)  It will be at most 1.

           XNORM

                     XNORM is DOUBLE PRECISION
                     The infinity-norm of X, when X is regarded as an NA x NW
                     real matrix.

           INFO

                     INFO is INTEGER
                     An error flag.  It will be set to zero if no error occurs,
                     a negative number if an argument is in error, or a positive
                     number if  ca A - w D  had to be perturbed.
                     The possible values are:
                     = 0: No error occurred, and (ca A - w D) did not have to be
                            perturbed.
                     = 1: (ca A - w D) had to be perturbed to make its smallest
                          (or only) singular value greater than SMIN.
                     NOTE: In the interests of speed, this routine does not
                           check the inputs for errors.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   double precision function dlangt (character NORM, integer N, double precision, dimension( * )
       DL, double precision, dimension( * ) D, double precision, dimension( * ) DU)
       DLANGT returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest
       absolute value of any element of a general tridiagonal matrix.

       Purpose:

            DLANGT  returns the value of the one norm,  or the Frobenius norm, or
            the  infinity norm,  or the  element of  largest absolute value  of a
            real tridiagonal matrix A.

       Returns:
           DLANGT

               DLANGT = ( max(abs(A(i,j))), NORM = 'M' or 'm'
                        (
                        ( norm1(A),         NORM = '1', 'O' or 'o'
                        (
                        ( normI(A),         NORM = 'I' or 'i'
                        (
                        ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

            where  norm1  denotes the  one norm of a matrix (maximum column sum),
            normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
            normF  denotes the  Frobenius norm of a matrix (square root of sum of
            squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.

       Parameters:
           NORM

                     NORM is CHARACTER*1
                     Specifies the value to be returned in DLANGT as described
                     above.

           N

                     N is INTEGER
                     The order of the matrix A.  N >= 0.  When N = 0, DLANGT is
                     set to zero.

           DL

                     DL is DOUBLE PRECISION array, dimension (N-1)
                     The (n-1) sub-diagonal elements of A.

           D

                     D is DOUBLE PRECISION array, dimension (N)
                     The diagonal elements of A.

           DU

                     DU is DOUBLE PRECISION array, dimension (N-1)
                     The (n-1) super-diagonal elements of A.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   double precision function dlanhs (character NORM, integer N, double precision, dimension( lda,
       * ) A, integer LDA, double precision, dimension( * ) WORK)
       DLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest
       absolute value of any element of an upper Hessenberg matrix.

       Purpose:

            DLANHS  returns the value of the one norm,  or the Frobenius norm, or
            the  infinity norm,  or the  element of  largest absolute value  of a
            Hessenberg matrix A.

       Returns:
           DLANHS

               DLANHS = ( max(abs(A(i,j))), NORM = 'M' or 'm'
                        (
                        ( norm1(A),         NORM = '1', 'O' or 'o'
                        (
                        ( normI(A),         NORM = 'I' or 'i'
                        (
                        ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

            where  norm1  denotes the  one norm of a matrix (maximum column sum),
            normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
            normF  denotes the  Frobenius norm of a matrix (square root of sum of
            squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.

       Parameters:
           NORM

                     NORM is CHARACTER*1
                     Specifies the value to be returned in DLANHS as described
                     above.

           N

                     N is INTEGER
                     The order of the matrix A.  N >= 0.  When N = 0, DLANHS is
                     set to zero.

           A

                     A is DOUBLE PRECISION array, dimension (LDA,N)
                     The n by n upper Hessenberg matrix A; the part of A below the
                     first sub-diagonal is not referenced.

           LDA

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

           WORK

                     WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
                     where LWORK >= N when NORM = 'I'; otherwise, WORK is not
                     referenced.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   double precision function dlansb (character NORM, character UPLO, integer N, integer K, double
       precision, dimension( ldab, * ) AB, integer LDAB, double precision, dimension( * ) WORK)
       DLANSB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or
       the element of largest absolute value of a symmetric band matrix.

       Purpose:

            DLANSB  returns the value of the one norm,  or the Frobenius norm, or
            the  infinity norm,  or the element of  largest absolute value  of an
            n by n symmetric band matrix A,  with k super-diagonals.

       Returns:
           DLANSB

               DLANSB = ( max(abs(A(i,j))), NORM = 'M' or 'm'
                        (
                        ( norm1(A),         NORM = '1', 'O' or 'o'
                        (
                        ( normI(A),         NORM = 'I' or 'i'
                        (
                        ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

            where  norm1  denotes the  one norm of a matrix (maximum column sum),
            normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
            normF  denotes the  Frobenius norm of a matrix (square root of sum of
            squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.

       Parameters:
           NORM

                     NORM is CHARACTER*1
                     Specifies the value to be returned in DLANSB as described
                     above.

           UPLO

                     UPLO is CHARACTER*1
                     Specifies whether the upper or lower triangular part of the
                     band matrix A is supplied.
                     = 'U':  Upper triangular part is supplied
                     = 'L':  Lower triangular part is supplied

           N

                     N is INTEGER
                     The order of the matrix A.  N >= 0.  When N = 0, DLANSB is
                     set to zero.

           K

                     K is INTEGER
                     The number of super-diagonals or sub-diagonals of the
                     band matrix A.  K >= 0.

           AB

                     AB is DOUBLE PRECISION array, dimension (LDAB,N)
                     The upper or lower triangle of the symmetric band matrix A,
                     stored in the first K+1 rows of AB.  The j-th column of A is
                     stored in the j-th column of the array AB as follows:
                     if UPLO = 'U', AB(k+1+i-j,j) = A(i,j) for max(1,j-k)<=i<=j;
                     if UPLO = 'L', AB(1+i-j,j)   = A(i,j) for j<=i<=min(n,j+k).

           LDAB

                     LDAB is INTEGER
                     The leading dimension of the array AB.  LDAB >= K+1.

           WORK

                     WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
                     where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
                     WORK is not referenced.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   double precision function dlansp (character NORM, character UPLO, integer N, double precision,
       dimension( * ) AP, double precision, dimension( * ) WORK)
       DLANSP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or
       the element of largest absolute value of a symmetric matrix supplied in packed form.

       Purpose:

            DLANSP  returns the value of the one norm,  or the Frobenius norm, or
            the  infinity norm,  or the  element of  largest absolute value  of a
            real symmetric matrix A,  supplied in packed form.

       Returns:
           DLANSP

               DLANSP = ( max(abs(A(i,j))), NORM = 'M' or 'm'
                        (
                        ( norm1(A),         NORM = '1', 'O' or 'o'
                        (
                        ( normI(A),         NORM = 'I' or 'i'
                        (
                        ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

            where  norm1  denotes the  one norm of a matrix (maximum column sum),
            normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
            normF  denotes the  Frobenius norm of a matrix (square root of sum of
            squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.

       Parameters:
           NORM

                     NORM is CHARACTER*1
                     Specifies the value to be returned in DLANSP as described
                     above.

           UPLO

                     UPLO is CHARACTER*1
                     Specifies whether the upper or lower triangular part of the
                     symmetric matrix A is supplied.
                     = 'U':  Upper triangular part of A is supplied
                     = 'L':  Lower triangular part of A is supplied

           N

                     N is INTEGER
                     The order of the matrix A.  N >= 0.  When N = 0, DLANSP is
                     set to zero.

           AP

                     AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
                     The upper or lower triangle of the symmetric matrix A, packed
                     columnwise in a linear array.  The j-th column of A is stored
                     in the array AP as follows:
                     if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
                     if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.

           WORK

                     WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
                     where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
                     WORK is not referenced.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   double precision function dlantb (character NORM, character UPLO, character DIAG, integer N,
       integer K, double precision, dimension( ldab, * ) AB, integer LDAB, double precision,
       dimension( * ) WORK)
       DLANTB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or
       the element of largest absolute value of a triangular band matrix.

       Purpose:

            DLANTB  returns the value of the one norm,  or the Frobenius norm, or
            the  infinity norm,  or the element of  largest absolute value  of an
            n by n triangular band matrix A,  with ( k + 1 ) diagonals.

       Returns:
           DLANTB

               DLANTB = ( max(abs(A(i,j))), NORM = 'M' or 'm'
                        (
                        ( norm1(A),         NORM = '1', 'O' or 'o'
                        (
                        ( normI(A),         NORM = 'I' or 'i'
                        (
                        ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

            where  norm1  denotes the  one norm of a matrix (maximum column sum),
            normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
            normF  denotes the  Frobenius norm of a matrix (square root of sum of
            squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.

       Parameters:
           NORM

                     NORM is CHARACTER*1
                     Specifies the value to be returned in DLANTB as described
                     above.

           UPLO

                     UPLO is CHARACTER*1
                     Specifies whether the matrix A is upper or lower triangular.
                     = 'U':  Upper triangular
                     = 'L':  Lower triangular

           DIAG

                     DIAG is CHARACTER*1
                     Specifies whether or not the matrix A is unit triangular.
                     = 'N':  Non-unit triangular
                     = 'U':  Unit triangular

           N

                     N is INTEGER
                     The order of the matrix A.  N >= 0.  When N = 0, DLANTB is
                     set to zero.

           K

                     K is INTEGER
                     The number of super-diagonals of the matrix A if UPLO = 'U',
                     or the number of sub-diagonals of the matrix A if UPLO = 'L'.
                     K >= 0.

           AB

                     AB is DOUBLE PRECISION array, dimension (LDAB,N)
                     The upper or lower triangular band matrix A, stored in the
                     first k+1 rows of AB.  The j-th column of A is stored
                     in the j-th column of the array AB as follows:
                     if UPLO = 'U', AB(k+1+i-j,j) = A(i,j) for max(1,j-k)<=i<=j;
                     if UPLO = 'L', AB(1+i-j,j)   = A(i,j) for j<=i<=min(n,j+k).
                     Note that when DIAG = 'U', the elements of the array AB
                     corresponding to the diagonal elements of the matrix A are
                     not referenced, but are assumed to be one.

           LDAB

                     LDAB is INTEGER
                     The leading dimension of the array AB.  LDAB >= K+1.

           WORK

                     WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
                     where LWORK >= N when NORM = 'I'; otherwise, WORK is not
                     referenced.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   double precision function dlantp (character NORM, character UPLO, character DIAG, integer N,
       double precision, dimension( * ) AP, double precision, dimension( * ) WORK)
       DLANTP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or
       the element of largest absolute value of a triangular matrix supplied in packed form.

       Purpose:

            DLANTP  returns the value of the one norm,  or the Frobenius norm, or
            the  infinity norm,  or the  element of  largest absolute value  of a
            triangular matrix A, supplied in packed form.

       Returns:
           DLANTP

               DLANTP = ( max(abs(A(i,j))), NORM = 'M' or 'm'
                        (
                        ( norm1(A),         NORM = '1', 'O' or 'o'
                        (
                        ( normI(A),         NORM = 'I' or 'i'
                        (
                        ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

            where  norm1  denotes the  one norm of a matrix (maximum column sum),
            normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
            normF  denotes the  Frobenius norm of a matrix (square root of sum of
            squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.

       Parameters:
           NORM

                     NORM is CHARACTER*1
                     Specifies the value to be returned in DLANTP as described
                     above.

           UPLO

                     UPLO is CHARACTER*1
                     Specifies whether the matrix A is upper or lower triangular.
                     = 'U':  Upper triangular
                     = 'L':  Lower triangular

           DIAG

                     DIAG is CHARACTER*1
                     Specifies whether or not the matrix A is unit triangular.
                     = 'N':  Non-unit triangular
                     = 'U':  Unit triangular

           N

                     N is INTEGER
                     The order of the matrix A.  N >= 0.  When N = 0, DLANTP is
                     set to zero.

           AP

                     AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
                     The upper or lower triangular matrix A, packed columnwise in
                     a linear array.  The j-th column of A is stored in the array
                     AP as follows:
                     if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
                     if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
                     Note that when DIAG = 'U', the elements of the array AP
                     corresponding to the diagonal elements of the matrix A are
                     not referenced, but are assumed to be one.

           WORK

                     WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
                     where LWORK >= N when NORM = 'I'; otherwise, WORK is not
                     referenced.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   double precision function dlantr (character NORM, character UPLO, character DIAG, integer M,
       integer N, double precision, dimension( lda, * ) A, integer LDA, double precision,
       dimension( * ) WORK)
       DLANTR returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or
       the element of largest absolute value of a trapezoidal or triangular matrix.

       Purpose:

            DLANTR  returns the value of the one norm,  or the Frobenius norm, or
            the  infinity norm,  or the  element of  largest absolute value  of a
            trapezoidal or triangular matrix A.

       Returns:
           DLANTR

               DLANTR = ( max(abs(A(i,j))), NORM = 'M' or 'm'
                        (
                        ( norm1(A),         NORM = '1', 'O' or 'o'
                        (
                        ( normI(A),         NORM = 'I' or 'i'
                        (
                        ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

            where  norm1  denotes the  one norm of a matrix (maximum column sum),
            normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
            normF  denotes the  Frobenius norm of a matrix (square root of sum of
            squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.

       Parameters:
           NORM

                     NORM is CHARACTER*1
                     Specifies the value to be returned in DLANTR as described
                     above.

           UPLO

                     UPLO is CHARACTER*1
                     Specifies whether the matrix A is upper or lower trapezoidal.
                     = 'U':  Upper trapezoidal
                     = 'L':  Lower trapezoidal
                     Note that A is triangular instead of trapezoidal if M = N.

           DIAG

                     DIAG is CHARACTER*1
                     Specifies whether or not the matrix A has unit diagonal.
                     = 'N':  Non-unit diagonal
                     = 'U':  Unit diagonal

           M

                     M is INTEGER
                     The number of rows of the matrix A.  M >= 0, and if
                     UPLO = 'U', M <= N.  When M = 0, DLANTR is set to zero.

           N

                     N is INTEGER
                     The number of columns of the matrix A.  N >= 0, and if
                     UPLO = 'L', N <= M.  When N = 0, DLANTR is set to zero.

           A

                     A is DOUBLE PRECISION array, dimension (LDA,N)
                     The trapezoidal matrix A (A is triangular if M = N).
                     If UPLO = 'U', the leading m by n upper trapezoidal part of
                     the array A contains the upper trapezoidal matrix, and the
                     strictly lower triangular part of A is not referenced.
                     If UPLO = 'L', the leading m by n lower trapezoidal part of
                     the array A contains the lower trapezoidal matrix, and the
                     strictly upper triangular part of A is not referenced.  Note
                     that when DIAG = 'U', the diagonal elements of A are not
                     referenced and are assumed to be one.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(M,1).

           WORK

                     WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
                     where LWORK >= M when NORM = 'I'; otherwise, WORK is not
                     referenced.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlanv2 (double precision A, double precision B, double precision C, double
       precision D, double precision RT1R, double precision RT1I, double precision RT2R, double
       precision RT2I, double precision CS, double precision SN)
       DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard
       form.

       Purpose:

            DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric
            matrix in standard form:

                 [ A  B ] = [ CS -SN ] [ AA  BB ] [ CS  SN ]
                 [ C  D ]   [ SN  CS ] [ CC  DD ] [-SN  CS ]

            where either
            1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or
            2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex
            conjugate eigenvalues.

       Parameters:
           A

                     A is DOUBLE PRECISION

           B

                     B is DOUBLE PRECISION

           C

                     C is DOUBLE PRECISION

           D

                     D is DOUBLE PRECISION
                     On entry, the elements of the input matrix.
                     On exit, they are overwritten by the elements of the
                     standardised Schur form.

           RT1R

                     RT1R is DOUBLE PRECISION

           RT1I

                     RT1I is DOUBLE PRECISION

           RT2R

                     RT2R is DOUBLE PRECISION

           RT2I

                     RT2I is DOUBLE PRECISION
                     The real and imaginary parts of the eigenvalues. If the
                     eigenvalues are a complex conjugate pair, RT1I > 0.

           CS

                     CS is DOUBLE PRECISION

           SN

                     SN is DOUBLE PRECISION
                     Parameters of the rotation matrix.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

       Further Details:

             Modified by V. Sima, Research Institute for Informatics, Bucharest,
             Romania, to reduce the risk of cancellation errors,
             when computing real eigenvalues, and to ensure, if possible, that
             abs(RT1R) >= abs(RT2R).

   subroutine dlapll (integer N, double precision, dimension( * ) X, integer INCX, double
       precision, dimension( * ) Y, integer INCY, double precision SSMIN)
       DLAPLL measures the linear dependence of two vectors.

       Purpose:

            Given two column vectors X and Y, let

                                 A = ( X Y ).

            The subroutine first computes the QR factorization of A = Q*R,
            and then computes the SVD of the 2-by-2 upper triangular matrix R.
            The smaller singular value of R is returned in SSMIN, which is used
            as the measurement of the linear dependency of the vectors X and Y.

       Parameters:
           N

                     N is INTEGER
                     The length of the vectors X and Y.

           X

                     X is DOUBLE PRECISION array,
                                    dimension (1+(N-1)*INCX)
                     On entry, X contains the N-vector X.
                     On exit, X is overwritten.

           INCX

                     INCX is INTEGER
                     The increment between successive elements of X. INCX > 0.

           Y

                     Y is DOUBLE PRECISION array,
                                    dimension (1+(N-1)*INCY)
                     On entry, Y contains the N-vector Y.
                     On exit, Y is overwritten.

           INCY

                     INCY is INTEGER
                     The increment between successive elements of Y. INCY > 0.

           SSMIN

                     SSMIN is DOUBLE PRECISION
                     The smallest singular value of the N-by-2 matrix A = ( X Y ).

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlapmr (logical FORWRD, integer M, integer N, double precision, dimension( ldx, * )
       X, integer LDX, integer, dimension( * ) K)
       DLAPMR rearranges rows of a matrix as specified by a permutation vector.

       Purpose:

            DLAPMR rearranges the rows of the M by N matrix X as specified
            by the permutation K(1),K(2),...,K(M) of the integers 1,...,M.
            If FORWRD = .TRUE.,  forward permutation:

                 X(K(I),*) is moved X(I,*) for I = 1,2,...,M.

            If FORWRD = .FALSE., backward permutation:

                 X(I,*) is moved to X(K(I),*) for I = 1,2,...,M.

       Parameters:
           FORWRD

                     FORWRD is LOGICAL
                     = .TRUE., forward permutation
                     = .FALSE., backward permutation

           M

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

           N

                     N is INTEGER
                     The number of columns of the matrix X. N >= 0.

           X

                     X is DOUBLE PRECISION array, dimension (LDX,N)
                     On entry, the M by N matrix X.
                     On exit, X contains the permuted matrix X.

           LDX

                     LDX is INTEGER
                     The leading dimension of the array X, LDX >= MAX(1,M).

           K

                     K is INTEGER array, dimension (M)
                     On entry, K contains the permutation vector. K is used as
                     internal workspace, but reset to its original value on
                     output.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlapmt (logical FORWRD, integer M, integer N, double precision, dimension( ldx, * )
       X, integer LDX, integer, dimension( * ) K)
       DLAPMT performs a forward or backward permutation of the columns of a matrix.

       Purpose:

            DLAPMT rearranges the columns of the M by N matrix X as specified
            by the permutation K(1),K(2),...,K(N) of the integers 1,...,N.
            If FORWRD = .TRUE.,  forward permutation:

                 X(*,K(J)) is moved X(*,J) for J = 1,2,...,N.

            If FORWRD = .FALSE., backward permutation:

                 X(*,J) is moved to X(*,K(J)) for J = 1,2,...,N.

       Parameters:
           FORWRD

                     FORWRD is LOGICAL
                     = .TRUE., forward permutation
                     = .FALSE., backward permutation

           M

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

           N

                     N is INTEGER
                     The number of columns of the matrix X. N >= 0.

           X

                     X is DOUBLE PRECISION array, dimension (LDX,N)
                     On entry, the M by N matrix X.
                     On exit, X contains the permuted matrix X.

           LDX

                     LDX is INTEGER
                     The leading dimension of the array X, LDX >= MAX(1,M).

           K

                     K is INTEGER array, dimension (N)
                     On entry, K contains the permutation vector. K is used as
                     internal workspace, but reset to its original value on
                     output.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlaqp2 (integer M, integer N, integer OFFSET, double precision, dimension( lda, * )
       A, integer LDA, integer, dimension( * ) JPVT, double precision, dimension( * ) TAU, double
       precision, dimension( * ) VN1, double precision, dimension( * ) VN2, double precision,
       dimension( * ) WORK)
       DLAQP2 computes a QR factorization with column pivoting of the matrix block.

       Purpose:

            DLAQP2 computes a QR factorization with column pivoting of
            the block A(OFFSET+1:M,1:N).
            The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.

       Parameters:
           M

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

           N

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

           OFFSET

                     OFFSET is INTEGER
                     The number of rows of the matrix A that must be pivoted
                     but no factorized. OFFSET >= 0.

           A

                     A is DOUBLE PRECISION array, dimension (LDA,N)
                     On entry, the M-by-N matrix A.
                     On exit, the upper triangle of block A(OFFSET+1:M,1:N) is
                     the triangular factor obtained; the elements in block
                     A(OFFSET+1:M,1:N) below the diagonal, together with the
                     array TAU, represent the orthogonal matrix Q as a product of
                     elementary reflectors. Block A(1:OFFSET,1:N) has been
                     accordingly pivoted, but no factorized.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A. LDA >= max(1,M).

           JPVT

                     JPVT is INTEGER array, dimension (N)
                     On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
                     to the front of A*P (a leading column); if JPVT(i) = 0,
                     the i-th column of A is a free column.
                     On exit, if JPVT(i) = k, then the i-th column of A*P
                     was the k-th column of A.

           TAU

                     TAU is DOUBLE PRECISION array, dimension (min(M,N))
                     The scalar factors of the elementary reflectors.

           VN1

                     VN1 is DOUBLE PRECISION array, dimension (N)
                     The vector with the partial column norms.

           VN2

                     VN2 is DOUBLE PRECISION array, dimension (N)
                     The vector with the exact column norms.

           WORK

                     WORK is DOUBLE PRECISION array, dimension (N)

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

       Contributors:
           G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain X. Sun, Computer
           Science Dept., Duke University, USA
            Partial column norm updating strategy modified on April 2011 Z. Drmac and Z.
           Bujanovic, Dept. of Mathematics, University of Zagreb, Croatia.

       References:
           LAPACK Working Note 176

   subroutine dlaqps (integer M, integer N, integer OFFSET, integer NB, integer KB, double
       precision, dimension( lda, * ) A, integer LDA, integer, dimension( * ) JPVT, double
       precision, dimension( * ) TAU, double precision, dimension( * ) VN1, double precision,
       dimension( * ) VN2, double precision, dimension( * ) AUXV, double precision, dimension(
       ldf, * ) F, integer LDF)
       DLAQPS computes a step of QR factorization with column pivoting of a real m-by-n matrix A
       by using BLAS level 3.

       Purpose:

            DLAQPS computes a step of QR factorization with column pivoting
            of a real M-by-N matrix A by using Blas-3.  It tries to factorize
            NB columns from A starting from the row OFFSET+1, and updates all
            of the matrix with Blas-3 xGEMM.

            In some cases, due to catastrophic cancellations, it cannot
            factorize NB columns.  Hence, the actual number of factorized
            columns is returned in KB.

            Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.

       Parameters:
           M

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

           N

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

           OFFSET

                     OFFSET is INTEGER
                     The number of rows of A that have been factorized in
                     previous steps.

           NB

                     NB is INTEGER
                     The number of columns to factorize.

           KB

                     KB is INTEGER
                     The number of columns actually factorized.

           A

                     A is DOUBLE PRECISION array, dimension (LDA,N)
                     On entry, the M-by-N matrix A.
                     On exit, block A(OFFSET+1:M,1:KB) is the triangular
                     factor obtained and block A(1:OFFSET,1:N) has been
                     accordingly pivoted, but no factorized.
                     The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has
                     been updated.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A. LDA >= max(1,M).

           JPVT

                     JPVT is INTEGER array, dimension (N)
                     JPVT(I) = K <==> Column K of the full matrix A has been
                     permuted into position I in AP.

           TAU

                     TAU is DOUBLE PRECISION array, dimension (KB)
                     The scalar factors of the elementary reflectors.

           VN1

                     VN1 is DOUBLE PRECISION array, dimension (N)
                     The vector with the partial column norms.

           VN2

                     VN2 is DOUBLE PRECISION array, dimension (N)
                     The vector with the exact column norms.

           AUXV

                     AUXV is DOUBLE PRECISION array, dimension (NB)
                     Auxiliar vector.

           F

                     F is DOUBLE PRECISION array, dimension (LDF,NB)
                     Matrix F**T = L*Y**T*A.

           LDF

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

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

       Contributors:
           G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain X. Sun, Computer
           Science Dept., Duke University, USA
            Partial column norm updating strategy modified on April 2011 Z. Drmac and Z.
           Bujanovic, Dept. of Mathematics, University of Zagreb, Croatia.

       References:
           LAPACK Working Note 176

   subroutine dlaqr0 (logical WANTT, logical WANTZ, integer N, integer ILO, integer IHI, double
       precision, dimension( ldh, * ) H, integer LDH, double precision, dimension( * ) WR, double
       precision, dimension( * ) WI, integer ILOZ, integer IHIZ, double precision, dimension(
       ldz, * ) Z, integer LDZ, double precision, dimension( * ) WORK, integer LWORK, integer
       INFO)
       DLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from
       the Schur decomposition.

       Purpose:

               DLAQR0 computes the eigenvalues of a Hessenberg matrix H
               and, optionally, the matrices T and Z from the Schur decomposition
               H = Z T Z**T, where T is an upper quasi-triangular matrix (the
               Schur form), and Z is the orthogonal matrix of Schur vectors.

               Optionally Z may be postmultiplied into an input orthogonal
               matrix Q so that this routine can give the Schur factorization
               of a matrix A which has been reduced to the Hessenberg form H
               by the orthogonal matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.

       Parameters:
           WANTT

                     WANTT is LOGICAL
                     = .TRUE. : the full Schur form T is required;
                     = .FALSE.: only eigenvalues are required.

           WANTZ

                     WANTZ is LOGICAL
                     = .TRUE. : the matrix of Schur vectors Z is required;
                     = .FALSE.: Schur vectors are not required.

           N

                     N is INTEGER
                      The order of the matrix H.  N .GE. 0.

           ILO

                     ILO is INTEGER

           IHI

                     IHI is INTEGER
                      It is assumed that H is already upper triangular in rows
                      and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
                      H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
                      previous call to DGEBAL, and then passed to DGEHRD when the
                      matrix output by DGEBAL is reduced to Hessenberg form.
                      Otherwise, ILO and IHI should be set to 1 and N,
                      respectively.  If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
                      If N = 0, then ILO = 1 and IHI = 0.

           H

                     H is DOUBLE PRECISION array, dimension (LDH,N)
                      On entry, the upper Hessenberg matrix H.
                      On exit, if INFO = 0 and WANTT is .TRUE., then H contains
                      the upper quasi-triangular matrix T from the Schur
                      decomposition (the Schur form); 2-by-2 diagonal blocks
                      (corresponding to complex conjugate pairs of eigenvalues)
                      are returned in standard form, with H(i,i) = H(i+1,i+1)
                      and H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and WANTT is
                      .FALSE., then the contents of H are unspecified on exit.
                      (The output value of H when INFO.GT.0 is given under the
                      description of INFO below.)

                      This subroutine may explicitly set H(i,j) = 0 for i.GT.j and
                      j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.

           LDH

                     LDH is INTEGER
                      The leading dimension of the array H. LDH .GE. max(1,N).

           WR

                     WR is DOUBLE PRECISION array, dimension (IHI)

           WI

                     WI is DOUBLE PRECISION array, dimension (IHI)
                      The real and imaginary parts, respectively, of the computed
                      eigenvalues of H(ILO:IHI,ILO:IHI) are stored in WR(ILO:IHI)
                      and WI(ILO:IHI). If two eigenvalues are computed as a
                      complex conjugate pair, they are stored in consecutive
                      elements of WR and WI, say the i-th and (i+1)th, with
                      WI(i) .GT. 0 and WI(i+1) .LT. 0. If WANTT is .TRUE., then
                      the eigenvalues are stored in the same order as on the
                      diagonal of the Schur form returned in H, with
                      WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal
                      block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
                      WI(i+1) = -WI(i).

           ILOZ

                     ILOZ is INTEGER

           IHIZ

                     IHIZ is INTEGER
                      Specify the rows of Z to which transformations must be
                      applied if WANTZ is .TRUE..
                      1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N.

           Z

                     Z is DOUBLE PRECISION array, dimension (LDZ,IHI)
                      If WANTZ is .FALSE., then Z is not referenced.
                      If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
                      replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
                      orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
                      (The output value of Z when INFO.GT.0 is given under
                      the description of INFO below.)

           LDZ

                     LDZ is INTEGER
                      The leading dimension of the array Z.  if WANTZ is .TRUE.
                      then LDZ.GE.MAX(1,IHIZ).  Otherwize, LDZ.GE.1.

           WORK

                     WORK is DOUBLE PRECISION array, dimension LWORK
                      On exit, if LWORK = -1, WORK(1) returns an estimate of
                      the optimal value for LWORK.

           LWORK

                     LWORK is INTEGER
                      The dimension of the array WORK.  LWORK .GE. max(1,N)
                      is sufficient, but LWORK typically as large as 6*N may
                      be required for optimal performance.  A workspace query
                      to determine the optimal workspace size is recommended.

                      If LWORK = -1, then DLAQR0 does a workspace query.
                      In this case, DLAQR0 checks the input parameters and
                      estimates the optimal workspace size for the given
                      values of N, ILO and IHI.  The estimate is returned
                      in WORK(1).  No error message related to LWORK is
                      issued by XERBLA.  Neither H nor Z are accessed.

           INFO

                     INFO is INTEGER
                        =  0:  successful exit
                      .GT. 0:  if INFO = i, DLAQR0 failed to compute all of
                           the eigenvalues.  Elements 1:ilo-1 and i+1:n of WR
                           and WI contain those eigenvalues which have been
                           successfully computed.  (Failures are rare.)

                           If INFO .GT. 0 and WANT is .FALSE., then on exit,
                           the remaining unconverged eigenvalues are the eigen-
                           values of the upper Hessenberg matrix rows and
                           columns ILO through INFO of the final, output
                           value of H.

                           If INFO .GT. 0 and WANTT is .TRUE., then on exit

                      (*)  (initial value of H)*U  = U*(final value of H)

                           where U is an orthogonal matrix.  The final
                           value of H is upper Hessenberg and quasi-triangular
                           in rows and columns INFO+1 through IHI.

                           If INFO .GT. 0 and WANTZ is .TRUE., then on exit

                             (final value of Z(ILO:IHI,ILOZ:IHIZ)
                              =  (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U

                           where U is the orthogonal matrix in (*) (regard-
                           less of the value of WANTT.)

                           If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
                           accessed.

       Contributors:
           Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

       References:
           K. Braman, R. Byers and R. Mathias, The Multi-Shift QR Algorithm Part I: Maintaining
           Well Focused Shifts, and Level 3 Performance, SIAM Journal of Matrix Analysis, volume
           23, pages 929--947, 2002.
            K. Braman, R. Byers and R. Mathias, The Multi-Shift QR Algorithm Part II: Aggressive
           Early Deflation, SIAM Journal of Matrix Analysis, volume 23, pages 948--973, 2002.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlaqr1 (integer N, double precision, dimension( ldh, * ) H, integer LDH, double
       precision SR1, double precision SI1, double precision SR2, double precision SI2, double
       precision, dimension( * ) V)
       DLAQR1 sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3
       matrix H and specified shifts.

       Purpose:

                 Given a 2-by-2 or 3-by-3 matrix H, DLAQR1 sets v to a
                 scalar multiple of the first column of the product

                 (*)  K = (H - (sr1 + i*si1)*I)*(H - (sr2 + i*si2)*I)

                 scaling to avoid overflows and most underflows. It
                 is assumed that either

                         1) sr1 = sr2 and si1 = -si2
                     or
                         2) si1 = si2 = 0.

                 This is useful for starting double implicit shift bulges
                 in the QR algorithm.

       Parameters:
           N

                     N is INTEGER
                         Order of the matrix H. N must be either 2 or 3.

           H

                     H is DOUBLE PRECISION array, dimension (LDH,N)
                         The 2-by-2 or 3-by-3 matrix H in (*).

           LDH

                     LDH is INTEGER
                         The leading dimension of H as declared in
                         the calling procedure.  LDH.GE.N

           SR1

                     SR1 is DOUBLE PRECISION

           SI1

                     SI1 is DOUBLE PRECISION

           SR2

                     SR2 is DOUBLE PRECISION

           SI2

                     SI2 is DOUBLE PRECISION
                         The shifts in (*).

           V

                     V is DOUBLE PRECISION array, dimension (N)
                         A scalar multiple of the first column of the
                         matrix K in (*).

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           June 2017

       Contributors:
           Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

   subroutine dlaqr2 (logical WANTT, logical WANTZ, integer N, integer KTOP, integer KBOT,
       integer NW, double precision, dimension( ldh, * ) H, integer LDH, integer ILOZ, integer
       IHIZ, double precision, dimension( ldz, * ) Z, integer LDZ, integer NS, integer ND, double
       precision, dimension( * ) SR, double precision, dimension( * ) SI, double precision,
       dimension( ldv, * ) V, integer LDV, integer NH, double precision, dimension( ldt, * ) T,
       integer LDT, integer NV, double precision, dimension( ldwv, * ) WV, integer LDWV, double
       precision, dimension( * ) WORK, integer LWORK)
       DLAQR2 performs the orthogonal similarity transformation of a Hessenberg matrix to detect
       and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive
       early deflation).

       Purpose:

               DLAQR2 is identical to DLAQR3 except that it avoids
               recursion by calling DLAHQR instead of DLAQR4.

               Aggressive early deflation:

               This subroutine accepts as input an upper Hessenberg matrix
               H and performs an orthogonal similarity transformation
               designed to detect and deflate fully converged eigenvalues from
               a trailing principal submatrix.  On output H has been over-
               written by a new Hessenberg matrix that is a perturbation of
               an orthogonal similarity transformation of H.  It is to be
               hoped that the final version of H has many zero subdiagonal
               entries.

       Parameters:
           WANTT

                     WANTT is LOGICAL
                     If .TRUE., then the Hessenberg matrix H is fully updated
                     so that the quasi-triangular Schur factor may be
                     computed (in cooperation with the calling subroutine).
                     If .FALSE., then only enough of H is updated to preserve
                     the eigenvalues.

           WANTZ

                     WANTZ is LOGICAL
                     If .TRUE., then the orthogonal matrix Z is updated so
                     so that the orthogonal Schur factor may be computed
                     (in cooperation with the calling subroutine).
                     If .FALSE., then Z is not referenced.

           N

                     N is INTEGER
                     The order of the matrix H and (if WANTZ is .TRUE.) the
                     order of the orthogonal matrix Z.

           KTOP

                     KTOP is INTEGER
                     It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0.
                     KBOT and KTOP together determine an isolated block
                     along the diagonal of the Hessenberg matrix.

           KBOT

                     KBOT is INTEGER
                     It is assumed without a check that either
                     KBOT = N or H(KBOT+1,KBOT)=0.  KBOT and KTOP together
                     determine an isolated block along the diagonal of the
                     Hessenberg matrix.

           NW

                     NW is INTEGER
                     Deflation window size.  1 .LE. NW .LE. (KBOT-KTOP+1).

           H

                     H is DOUBLE PRECISION array, dimension (LDH,N)
                     On input the initial N-by-N section of H stores the
                     Hessenberg matrix undergoing aggressive early deflation.
                     On output H has been transformed by an orthogonal
                     similarity transformation, perturbed, and the returned
                     to Hessenberg form that (it is to be hoped) has some
                     zero subdiagonal entries.

           LDH

                     LDH is INTEGER
                     Leading dimension of H just as declared in the calling
                     subroutine.  N .LE. LDH

           ILOZ

                     ILOZ is INTEGER

           IHIZ

                     IHIZ is INTEGER
                     Specify the rows of Z to which transformations must be
                     applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.

           Z

                     Z is DOUBLE PRECISION array, dimension (LDZ,N)
                     IF WANTZ is .TRUE., then on output, the orthogonal
                     similarity transformation mentioned above has been
                     accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right.
                     If WANTZ is .FALSE., then Z is unreferenced.

           LDZ

                     LDZ is INTEGER
                     The leading dimension of Z just as declared in the
                     calling subroutine.  1 .LE. LDZ.

           NS

                     NS is INTEGER
                     The number of unconverged (ie approximate) eigenvalues
                     returned in SR and SI that may be used as shifts by the
                     calling subroutine.

           ND

                     ND is INTEGER
                     The number of converged eigenvalues uncovered by this
                     subroutine.

           SR

                     SR is DOUBLE PRECISION array, dimension (KBOT)

           SI

                     SI is DOUBLE PRECISION array, dimension (KBOT)
                     On output, the real and imaginary parts of approximate
                     eigenvalues that may be used for shifts are stored in
                     SR(KBOT-ND-NS+1) through SR(KBOT-ND) and
                     SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively.
                     The real and imaginary parts of converged eigenvalues
                     are stored in SR(KBOT-ND+1) through SR(KBOT) and
                     SI(KBOT-ND+1) through SI(KBOT), respectively.

           V

                     V is DOUBLE PRECISION array, dimension (LDV,NW)
                     An NW-by-NW work array.

           LDV

                     LDV is INTEGER
                     The leading dimension of V just as declared in the
                     calling subroutine.  NW .LE. LDV

           NH

                     NH is INTEGER
                     The number of columns of T.  NH.GE.NW.

           T

                     T is DOUBLE PRECISION array, dimension (LDT,NW)

           LDT

                     LDT is INTEGER
                     The leading dimension of T just as declared in the
                     calling subroutine.  NW .LE. LDT

           NV

                     NV is INTEGER
                     The number of rows of work array WV available for
                     workspace.  NV.GE.NW.

           WV

                     WV is DOUBLE PRECISION array, dimension (LDWV,NW)

           LDWV

                     LDWV is INTEGER
                     The leading dimension of W just as declared in the
                     calling subroutine.  NW .LE. LDV

           WORK

                     WORK is DOUBLE PRECISION array, dimension (LWORK)
                     On exit, WORK(1) is set to an estimate of the optimal value
                     of LWORK for the given values of N, NW, KTOP and KBOT.

           LWORK

                     LWORK is INTEGER
                     The dimension of the work array WORK.  LWORK = 2*NW
                     suffices, but greater efficiency may result from larger
                     values of LWORK.

                     If LWORK = -1, then a workspace query is assumed; DLAQR2
                     only estimates the optimal workspace size for the given
                     values of N, NW, KTOP and KBOT.  The estimate is returned
                     in WORK(1).  No error message related to LWORK is issued
                     by XERBLA.  Neither H nor Z are accessed.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           June 2017

       Contributors:
           Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

   subroutine dlaqr3 (logical WANTT, logical WANTZ, integer N, integer KTOP, integer KBOT,
       integer NW, double precision, dimension( ldh, * ) H, integer LDH, integer ILOZ, integer
       IHIZ, double precision, dimension( ldz, * ) Z, integer LDZ, integer NS, integer ND, double
       precision, dimension( * ) SR, double precision, dimension( * ) SI, double precision,
       dimension( ldv, * ) V, integer LDV, integer NH, double precision, dimension( ldt, * ) T,
       integer LDT, integer NV, double precision, dimension( ldwv, * ) WV, integer LDWV, double
       precision, dimension( * ) WORK, integer LWORK)
       DLAQR3 performs the orthogonal similarity transformation of a Hessenberg matrix to detect
       and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive
       early deflation).

       Purpose:

               Aggressive early deflation:

               DLAQR3 accepts as input an upper Hessenberg matrix
               H and performs an orthogonal similarity transformation
               designed to detect and deflate fully converged eigenvalues from
               a trailing principal submatrix.  On output H has been over-
               written by a new Hessenberg matrix that is a perturbation of
               an orthogonal similarity transformation of H.  It is to be
               hoped that the final version of H has many zero subdiagonal
               entries.

       Parameters:
           WANTT

                     WANTT is LOGICAL
                     If .TRUE., then the Hessenberg matrix H is fully updated
                     so that the quasi-triangular Schur factor may be
                     computed (in cooperation with the calling subroutine).
                     If .FALSE., then only enough of H is updated to preserve
                     the eigenvalues.

           WANTZ

                     WANTZ is LOGICAL
                     If .TRUE., then the orthogonal matrix Z is updated so
                     so that the orthogonal Schur factor may be computed
                     (in cooperation with the calling subroutine).
                     If .FALSE., then Z is not referenced.

           N

                     N is INTEGER
                     The order of the matrix H and (if WANTZ is .TRUE.) the
                     order of the orthogonal matrix Z.

           KTOP

                     KTOP is INTEGER
                     It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0.
                     KBOT and KTOP together determine an isolated block
                     along the diagonal of the Hessenberg matrix.

           KBOT

                     KBOT is INTEGER
                     It is assumed without a check that either
                     KBOT = N or H(KBOT+1,KBOT)=0.  KBOT and KTOP together
                     determine an isolated block along the diagonal of the
                     Hessenberg matrix.

           NW

                     NW is INTEGER
                     Deflation window size.  1 .LE. NW .LE. (KBOT-KTOP+1).

           H

                     H is DOUBLE PRECISION array, dimension (LDH,N)
                     On input the initial N-by-N section of H stores the
                     Hessenberg matrix undergoing aggressive early deflation.
                     On output H has been transformed by an orthogonal
                     similarity transformation, perturbed, and the returned
                     to Hessenberg form that (it is to be hoped) has some
                     zero subdiagonal entries.

           LDH

                     LDH is INTEGER
                     Leading dimension of H just as declared in the calling
                     subroutine.  N .LE. LDH

           ILOZ

                     ILOZ is INTEGER

           IHIZ

                     IHIZ is INTEGER
                     Specify the rows of Z to which transformations must be
                     applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.

           Z

                     Z is DOUBLE PRECISION array, dimension (LDZ,N)
                     IF WANTZ is .TRUE., then on output, the orthogonal
                     similarity transformation mentioned above has been
                     accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right.
                     If WANTZ is .FALSE., then Z is unreferenced.

           LDZ

                     LDZ is INTEGER
                     The leading dimension of Z just as declared in the
                     calling subroutine.  1 .LE. LDZ.

           NS

                     NS is INTEGER
                     The number of unconverged (ie approximate) eigenvalues
                     returned in SR and SI that may be used as shifts by the
                     calling subroutine.

           ND

                     ND is INTEGER
                     The number of converged eigenvalues uncovered by this
                     subroutine.

           SR

                     SR is DOUBLE PRECISION array, dimension (KBOT)

           SI

                     SI is DOUBLE PRECISION array, dimension (KBOT)
                     On output, the real and imaginary parts of approximate
                     eigenvalues that may be used for shifts are stored in
                     SR(KBOT-ND-NS+1) through SR(KBOT-ND) and
                     SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively.
                     The real and imaginary parts of converged eigenvalues
                     are stored in SR(KBOT-ND+1) through SR(KBOT) and
                     SI(KBOT-ND+1) through SI(KBOT), respectively.

           V

                     V is DOUBLE PRECISION array, dimension (LDV,NW)
                     An NW-by-NW work array.

           LDV

                     LDV is INTEGER
                     The leading dimension of V just as declared in the
                     calling subroutine.  NW .LE. LDV

           NH

                     NH is INTEGER
                     The number of columns of T.  NH.GE.NW.

           T

                     T is DOUBLE PRECISION array, dimension (LDT,NW)

           LDT

                     LDT is INTEGER
                     The leading dimension of T just as declared in the
                     calling subroutine.  NW .LE. LDT

           NV

                     NV is INTEGER
                     The number of rows of work array WV available for
                     workspace.  NV.GE.NW.

           WV

                     WV is DOUBLE PRECISION array, dimension (LDWV,NW)

           LDWV

                     LDWV is INTEGER
                     The leading dimension of W just as declared in the
                     calling subroutine.  NW .LE. LDV

           WORK

                     WORK is DOUBLE PRECISION array, dimension (LWORK)
                     On exit, WORK(1) is set to an estimate of the optimal value
                     of LWORK for the given values of N, NW, KTOP and KBOT.

           LWORK

                     LWORK is INTEGER
                     The dimension of the work array WORK.  LWORK = 2*NW
                     suffices, but greater efficiency may result from larger
                     values of LWORK.

                     If LWORK = -1, then a workspace query is assumed; DLAQR3
                     only estimates the optimal workspace size for the given
                     values of N, NW, KTOP and KBOT.  The estimate is returned
                     in WORK(1).  No error message related to LWORK is issued
                     by XERBLA.  Neither H nor Z are accessed.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           June 2016

       Contributors:
           Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

   subroutine dlaqr4 (logical WANTT, logical WANTZ, integer N, integer ILO, integer IHI, double
       precision, dimension( ldh, * ) H, integer LDH, double precision, dimension( * ) WR, double
       precision, dimension( * ) WI, integer ILOZ, integer IHIZ, double precision, dimension(
       ldz, * ) Z, integer LDZ, double precision, dimension( * ) WORK, integer LWORK, integer
       INFO)
       DLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from
       the Schur decomposition.

       Purpose:

               DLAQR4 implements one level of recursion for DLAQR0.
               It is a complete implementation of the small bulge multi-shift
               QR algorithm.  It may be called by DLAQR0 and, for large enough
               deflation window size, it may be called by DLAQR3.  This
               subroutine is identical to DLAQR0 except that it calls DLAQR2
               instead of DLAQR3.

               DLAQR4 computes the eigenvalues of a Hessenberg matrix H
               and, optionally, the matrices T and Z from the Schur decomposition
               H = Z T Z**T, where T is an upper quasi-triangular matrix (the
               Schur form), and Z is the orthogonal matrix of Schur vectors.

               Optionally Z may be postmultiplied into an input orthogonal
               matrix Q so that this routine can give the Schur factorization
               of a matrix A which has been reduced to the Hessenberg form H
               by the orthogonal matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.

       Parameters:
           WANTT

                     WANTT is LOGICAL
                     = .TRUE. : the full Schur form T is required;
                     = .FALSE.: only eigenvalues are required.

           WANTZ

                     WANTZ is LOGICAL
                     = .TRUE. : the matrix of Schur vectors Z is required;
                     = .FALSE.: Schur vectors are not required.

           N

                     N is INTEGER
                      The order of the matrix H.  N .GE. 0.

           ILO

                     ILO is INTEGER

           IHI

                     IHI is INTEGER
                      It is assumed that H is already upper triangular in rows
                      and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
                      H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
                      previous call to DGEBAL, and then passed to DGEHRD when the
                      matrix output by DGEBAL is reduced to Hessenberg form.
                      Otherwise, ILO and IHI should be set to 1 and N,
                      respectively.  If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
                      If N = 0, then ILO = 1 and IHI = 0.

           H

                     H is DOUBLE PRECISION array, dimension (LDH,N)
                      On entry, the upper Hessenberg matrix H.
                      On exit, if INFO = 0 and WANTT is .TRUE., then H contains
                      the upper quasi-triangular matrix T from the Schur
                      decomposition (the Schur form); 2-by-2 diagonal blocks
                      (corresponding to complex conjugate pairs of eigenvalues)
                      are returned in standard form, with H(i,i) = H(i+1,i+1)
                      and H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and WANTT is
                      .FALSE., then the contents of H are unspecified on exit.
                      (The output value of H when INFO.GT.0 is given under the
                      description of INFO below.)

                      This subroutine may explicitly set H(i,j) = 0 for i.GT.j and
                      j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.

           LDH

                     LDH is INTEGER
                      The leading dimension of the array H. LDH .GE. max(1,N).

           WR

                     WR is DOUBLE PRECISION array, dimension (IHI)

           WI

                     WI is DOUBLE PRECISION array, dimension (IHI)
                      The real and imaginary parts, respectively, of the computed
                      eigenvalues of H(ILO:IHI,ILO:IHI) are stored in WR(ILO:IHI)
                      and WI(ILO:IHI). If two eigenvalues are computed as a
                      complex conjugate pair, they are stored in consecutive
                      elements of WR and WI, say the i-th and (i+1)th, with
                      WI(i) .GT. 0 and WI(i+1) .LT. 0. If WANTT is .TRUE., then
                      the eigenvalues are stored in the same order as on the
                      diagonal of the Schur form returned in H, with
                      WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal
                      block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
                      WI(i+1) = -WI(i).

           ILOZ

                     ILOZ is INTEGER

           IHIZ

                     IHIZ is INTEGER
                      Specify the rows of Z to which transformations must be
                      applied if WANTZ is .TRUE..
                      1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N.

           Z

                     Z is DOUBLE PRECISION array, dimension (LDZ,IHI)
                      If WANTZ is .FALSE., then Z is not referenced.
                      If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
                      replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
                      orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
                      (The output value of Z when INFO.GT.0 is given under
                      the description of INFO below.)

           LDZ

                     LDZ is INTEGER
                      The leading dimension of the array Z.  if WANTZ is .TRUE.
                      then LDZ.GE.MAX(1,IHIZ).  Otherwize, LDZ.GE.1.

           WORK

                     WORK is DOUBLE PRECISION array, dimension LWORK
                      On exit, if LWORK = -1, WORK(1) returns an estimate of
                      the optimal value for LWORK.

           LWORK

                     LWORK is INTEGER
                      The dimension of the array WORK.  LWORK .GE. max(1,N)
                      is sufficient, but LWORK typically as large as 6*N may
                      be required for optimal performance.  A workspace query
                      to determine the optimal workspace size is recommended.

                      If LWORK = -1, then DLAQR4 does a workspace query.
                      In this case, DLAQR4 checks the input parameters and
                      estimates the optimal workspace size for the given
                      values of N, ILO and IHI.  The estimate is returned
                      in WORK(1).  No error message related to LWORK is
                      issued by XERBLA.  Neither H nor Z are accessed.

           INFO

                     INFO is INTEGER
                        =  0:  successful exit
                      .GT. 0:  if INFO = i, DLAQR4 failed to compute all of
                           the eigenvalues.  Elements 1:ilo-1 and i+1:n of WR
                           and WI contain those eigenvalues which have been
                           successfully computed.  (Failures are rare.)

                           If INFO .GT. 0 and WANT is .FALSE., then on exit,
                           the remaining unconverged eigenvalues are the eigen-
                           values of the upper Hessenberg matrix rows and
                           columns ILO through INFO of the final, output
                           value of H.

                           If INFO .GT. 0 and WANTT is .TRUE., then on exit

                      (*)  (initial value of H)*U  = U*(final value of H)

                           where U is a orthogonal matrix.  The final
                           value of  H is upper Hessenberg and triangular in
                           rows and columns INFO+1 through IHI.

                           If INFO .GT. 0 and WANTZ is .TRUE., then on exit

                             (final value of Z(ILO:IHI,ILOZ:IHIZ)
                              =  (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U

                           where U is the orthogonal matrix in (*) (regard-
                           less of the value of WANTT.)

                           If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
                           accessed.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

       Contributors:
           Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

       References:
           K. Braman, R. Byers and R. Mathias, The Multi-Shift QR Algorithm Part I: Maintaining
           Well Focused Shifts, and Level 3 Performance, SIAM Journal of Matrix Analysis, volume
           23, pages 929--947, 2002.
            K. Braman, R. Byers and R. Mathias, The Multi-Shift QR Algorithm Part II: Aggressive
           Early Deflation, SIAM Journal of Matrix Analysis, volume 23, pages 948--973, 2002.

   subroutine dlaqr5 (logical WANTT, logical WANTZ, integer KACC22, integer N, integer KTOP,
       integer KBOT, integer NSHFTS, double precision, dimension( * ) SR, double precision,
       dimension( * ) SI, double precision, dimension( ldh, * ) H, integer LDH, integer ILOZ,
       integer IHIZ, double precision, dimension( ldz, * ) Z, integer LDZ, double precision,
       dimension( ldv, * ) V, integer LDV, double precision, dimension( ldu, * ) U, integer LDU,
       integer NV, double precision, dimension( ldwv, * ) WV, integer LDWV, integer NH, double
       precision, dimension( ldwh, * ) WH, integer LDWH)
       DLAQR5 performs a single small-bulge multi-shift QR sweep.

       Purpose:

               DLAQR5, called by DLAQR0, performs a
               single small-bulge multi-shift QR sweep.

       Parameters:
           WANTT

                     WANTT is LOGICAL
                        WANTT = .true. if the quasi-triangular Schur factor
                        is being computed.  WANTT is set to .false. otherwise.

           WANTZ

                     WANTZ is LOGICAL
                        WANTZ = .true. if the orthogonal Schur factor is being
                        computed.  WANTZ is set to .false. otherwise.

           KACC22

                     KACC22 is INTEGER with value 0, 1, or 2.
                        Specifies the computation mode of far-from-diagonal
                        orthogonal updates.
                   = 0: DLAQR5 does not accumulate reflections and does not
                        use matrix-matrix multiply to update far-from-diagonal
                        matrix entries.
                   = 1: DLAQR5 accumulates reflections and uses matrix-matrix
                        multiply to update the far-from-diagonal matrix entries.
                   = 2: DLAQR5 accumulates reflections, uses matrix-matrix
                        multiply to update the far-from-diagonal matrix entries,
                        and takes advantage of 2-by-2 block structure during
                        matrix multiplies.

           N

                     N is INTEGER
                        N is the order of the Hessenberg matrix H upon which this
                        subroutine operates.

           KTOP

                     KTOP is INTEGER

           KBOT

                     KBOT is INTEGER
                        These are the first and last rows and columns of an
                        isolated diagonal block upon which the QR sweep is to be
                        applied. It is assumed without a check that
                                  either KTOP = 1  or   H(KTOP,KTOP-1) = 0
                        and
                                  either KBOT = N  or   H(KBOT+1,KBOT) = 0.

           NSHFTS

                     NSHFTS is INTEGER
                        NSHFTS gives the number of simultaneous shifts.  NSHFTS
                        must be positive and even.

           SR

                     SR is DOUBLE PRECISION array, dimension (NSHFTS)

           SI

                     SI is DOUBLE PRECISION array, dimension (NSHFTS)
                        SR contains the real parts and SI contains the imaginary
                        parts of the NSHFTS shifts of origin that define the
                        multi-shift QR sweep.  On output SR and SI may be
                        reordered.

           H

                     H is DOUBLE PRECISION array, dimension (LDH,N)
                        On input H contains a Hessenberg matrix.  On output a
                        multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied
                        to the isolated diagonal block in rows and columns KTOP
                        through KBOT.

           LDH

                     LDH is INTEGER
                        LDH is the leading dimension of H just as declared in the
                        calling procedure.  LDH.GE.MAX(1,N).

           ILOZ

                     ILOZ is INTEGER

           IHIZ

                     IHIZ is INTEGER
                        Specify the rows of Z to which transformations must be
                        applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N

           Z

                     Z is DOUBLE PRECISION array, dimension (LDZ,IHIZ)
                        If WANTZ = .TRUE., then the QR Sweep orthogonal
                        similarity transformation is accumulated into
                        Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right.
                        If WANTZ = .FALSE., then Z is unreferenced.

           LDZ

                     LDZ is INTEGER
                        LDA is the leading dimension of Z just as declared in
                        the calling procedure. LDZ.GE.N.

           V

                     V is DOUBLE PRECISION array, dimension (LDV,NSHFTS/2)

           LDV

                     LDV is INTEGER
                        LDV is the leading dimension of V as declared in the
                        calling procedure.  LDV.GE.3.

           U

                     U is DOUBLE PRECISION array, dimension (LDU,3*NSHFTS-3)

           LDU

                     LDU is INTEGER
                        LDU is the leading dimension of U just as declared in the
                        in the calling subroutine.  LDU.GE.3*NSHFTS-3.

           NH

                     NH is INTEGER
                        NH is the number of columns in array WH available for
                        workspace. NH.GE.1.

           WH

                     WH is DOUBLE PRECISION array, dimension (LDWH,NH)

           LDWH

                     LDWH is INTEGER
                        Leading dimension of WH just as declared in the
                        calling procedure.  LDWH.GE.3*NSHFTS-3.

           NV

                     NV is INTEGER
                        NV is the number of rows in WV agailable for workspace.
                        NV.GE.1.

           WV

                     WV is DOUBLE PRECISION array, dimension (LDWV,3*NSHFTS-3)

           LDWV

                     LDWV is INTEGER
                        LDWV is the leading dimension of WV as declared in the
                        in the calling subroutine.  LDWV.GE.NV.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           June 2016

       Contributors:
           Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

       References:
           K. Braman, R. Byers and R. Mathias, The Multi-Shift QR Algorithm Part I: Maintaining
           Well Focused Shifts, and Level 3 Performance, SIAM Journal of Matrix Analysis, volume
           23, pages 929--947, 2002.

   subroutine dlaqsb (character UPLO, integer N, integer KD, double precision, dimension( ldab, *
       ) AB, integer LDAB, double precision, dimension( * ) S, double precision SCOND, double
       precision AMAX, character EQUED)
       DLAQSB scales a symmetric/Hermitian band matrix, using scaling factors computed by spbequ.

       Purpose:

            DLAQSB equilibrates a symmetric band matrix A using the scaling
            factors in the vector S.

       Parameters:
           UPLO

                     UPLO is CHARACTER*1
                     Specifies whether the upper or lower triangular part of the
                     symmetric matrix A is stored.
                     = 'U':  Upper triangular
                     = 'L':  Lower triangular

           N

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

           KD

                     KD is INTEGER
                     The number of super-diagonals of the matrix A if UPLO = 'U',
                     or the number of sub-diagonals if UPLO = 'L'.  KD >= 0.

           AB

                     AB is DOUBLE PRECISION array, dimension (LDAB,N)
                     On entry, the upper or lower triangle of the symmetric band
                     matrix A, stored in the first KD+1 rows of the array.  The
                     j-th column of A is stored in the j-th column of the array AB
                     as follows:
                     if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
                     if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).

                     On exit, if INFO = 0, the triangular factor U or L from the
                     Cholesky factorization A = U**T*U or A = L*L**T of the band
                     matrix A, in the same storage format as A.

           LDAB

                     LDAB is INTEGER
                     The leading dimension of the array AB.  LDAB >= KD+1.

           S

                     S is DOUBLE PRECISION array, dimension (N)
                     The scale factors for A.

           SCOND

                     SCOND is DOUBLE PRECISION
                     Ratio of the smallest S(i) to the largest S(i).

           AMAX

                     AMAX is DOUBLE PRECISION
                     Absolute value of largest matrix entry.

           EQUED

                     EQUED is CHARACTER*1
                     Specifies whether or not equilibration was done.
                     = 'N':  No equilibration.
                     = 'Y':  Equilibration was done, i.e., A has been replaced by
                             diag(S) * A * diag(S).

       Internal Parameters:

             THRESH is a threshold value used to decide if scaling should be done
             based on the ratio of the scaling factors.  If SCOND < THRESH,
             scaling is done.

             LARGE and SMALL are threshold values used to decide if scaling should
             be done based on the absolute size of the largest matrix element.
             If AMAX > LARGE or AMAX < SMALL, scaling is done.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlaqsp (character UPLO, integer N, double precision, dimension( * ) AP, double
       precision, dimension( * ) S, double precision SCOND, double precision AMAX, character
       EQUED)
       DLAQSP scales a symmetric/Hermitian matrix in packed storage, using scaling factors
       computed by sppequ.

       Purpose:

            DLAQSP equilibrates a symmetric matrix A using the scaling factors
            in the vector S.

       Parameters:
           UPLO

                     UPLO is CHARACTER*1
                     Specifies whether the upper or lower triangular part of the
                     symmetric matrix A is stored.
                     = 'U':  Upper triangular
                     = 'L':  Lower triangular

           N

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

           AP

                     AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
                     On entry, the upper or lower triangle of the symmetric matrix
                     A, packed columnwise in a linear array.  The j-th column of A
                     is stored in the array AP as follows:
                     if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
                     if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.

                     On exit, the equilibrated matrix:  diag(S) * A * diag(S), in
                     the same storage format as A.

           S

                     S is DOUBLE PRECISION array, dimension (N)
                     The scale factors for A.

           SCOND

                     SCOND is DOUBLE PRECISION
                     Ratio of the smallest S(i) to the largest S(i).

           AMAX

                     AMAX is DOUBLE PRECISION
                     Absolute value of largest matrix entry.

           EQUED

                     EQUED is CHARACTER*1
                     Specifies whether or not equilibration was done.
                     = 'N':  No equilibration.
                     = 'Y':  Equilibration was done, i.e., A has been replaced by
                             diag(S) * A * diag(S).

       Internal Parameters:

             THRESH is a threshold value used to decide if scaling should be done
             based on the ratio of the scaling factors.  If SCOND < THRESH,
             scaling is done.

             LARGE and SMALL are threshold values used to decide if scaling should
             be done based on the absolute size of the largest matrix element.
             If AMAX > LARGE or AMAX < SMALL, scaling is done.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlaqtr (logical LTRAN, logical LREAL, integer N, double precision, dimension( ldt,
       * ) T, integer LDT, double precision, dimension( * ) B, double precision W, double
       precision SCALE, double precision, dimension( * ) X, double precision, dimension( * )
       WORK, integer INFO)
       DLAQTR solves a real quasi-triangular system of equations, or a complex quasi-triangular
       system of special form, in real arithmetic.

       Purpose:

            DLAQTR solves the real quasi-triangular system

                         op(T)*p = scale*c,               if LREAL = .TRUE.

            or the complex quasi-triangular systems

                       op(T + iB)*(p+iq) = scale*(c+id),  if LREAL = .FALSE.

            in real arithmetic, where T is upper quasi-triangular.
            If LREAL = .FALSE., then the first diagonal block of T must be
            1 by 1, B is the specially structured matrix

                           B = [ b(1) b(2) ... b(n) ]
                               [       w            ]
                               [           w        ]
                               [              .     ]
                               [                 w  ]

            op(A) = A or A**T, A**T denotes the transpose of
            matrix A.

            On input, X = [ c ].  On output, X = [ p ].
                          [ d ]                  [ q ]

            This subroutine is designed for the condition number estimation
            in routine DTRSNA.

       Parameters:
           LTRAN

                     LTRAN is LOGICAL
                     On entry, LTRAN specifies the option of conjugate transpose:
                        = .FALSE.,    op(T+i*B) = T+i*B,
                        = .TRUE.,     op(T+i*B) = (T+i*B)**T.

           LREAL

                     LREAL is LOGICAL
                     On entry, LREAL specifies the input matrix structure:
                        = .FALSE.,    the input is complex
                        = .TRUE.,     the input is real

           N

                     N is INTEGER
                     On entry, N specifies the order of T+i*B. N >= 0.

           T

                     T is DOUBLE PRECISION array, dimension (LDT,N)
                     On entry, T contains a matrix in Schur canonical form.
                     If LREAL = .FALSE., then the first diagonal block of T mu
                     be 1 by 1.

           LDT

                     LDT is INTEGER
                     The leading dimension of the matrix T. LDT >= max(1,N).

           B

                     B is DOUBLE PRECISION array, dimension (N)
                     On entry, B contains the elements to form the matrix
                     B as described above.
                     If LREAL = .TRUE., B is not referenced.

           W

                     W is DOUBLE PRECISION
                     On entry, W is the diagonal element of the matrix B.
                     If LREAL = .TRUE., W is not referenced.

           SCALE

                     SCALE is DOUBLE PRECISION
                     On exit, SCALE is the scale factor.

           X

                     X is DOUBLE PRECISION array, dimension (2*N)
                     On entry, X contains the right hand side of the system.
                     On exit, X is overwritten by the solution.

           WORK

                     WORK is DOUBLE PRECISION array, dimension (N)

           INFO

                     INFO is INTEGER
                     On exit, INFO is set to
                        0: successful exit.
                          1: the some diagonal 1 by 1 block has been perturbed by
                             a small number SMIN to keep nonsingularity.
                          2: the some diagonal 2 by 2 block has been perturbed by
                             a small number in DLALN2 to keep nonsingularity.
                     NOTE: In the interests of speed, this routine does not
                           check the inputs for errors.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlar1v (integer N, integer B1, integer BN, double precision LAMBDA, double
       precision, dimension( * ) D, double precision, dimension( * ) L, double precision,
       dimension( * ) LD, double precision, dimension( * ) LLD, double precision PIVMIN, double
       precision GAPTOL, double precision, dimension( * ) Z, logical WANTNC, integer NEGCNT,
       double precision ZTZ, double precision MINGMA, integer R, integer, dimension( * ) ISUPPZ,
       double precision NRMINV, double precision RESID, double precision RQCORR, double
       precision, dimension( * ) WORK)
       DLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1
       through bn of the tridiagonal matrix LDLT - λI.

       Purpose:

            DLAR1V computes the (scaled) r-th column of the inverse of
            the sumbmatrix in rows B1 through BN of the tridiagonal matrix
            L D L**T - sigma I. When sigma is close to an eigenvalue, the
            computed vector is an accurate eigenvector. Usually, r corresponds
            to the index where the eigenvector is largest in magnitude.
            The following steps accomplish this computation :
            (a) Stationary qd transform,  L D L**T - sigma I = L(+) D(+) L(+)**T,
            (b) Progressive qd transform, L D L**T - sigma I = U(-) D(-) U(-)**T,
            (c) Computation of the diagonal elements of the inverse of
                L D L**T - sigma I by combining the above transforms, and choosing
                r as the index where the diagonal of the inverse is (one of the)
                largest in magnitude.
            (d) Computation of the (scaled) r-th column of the inverse using the
                twisted factorization obtained by combining the top part of the
                the stationary and the bottom part of the progressive transform.

       Parameters:
           N

                     N is INTEGER
                      The order of the matrix L D L**T.

           B1

                     B1 is INTEGER
                      First index of the submatrix of L D L**T.

           BN

                     BN is INTEGER
                      Last index of the submatrix of L D L**T.

           LAMBDA

                     LAMBDA is DOUBLE PRECISION
                      The shift. In order to compute an accurate eigenvector,
                      LAMBDA should be a good approximation to an eigenvalue
                      of L D L**T.

           L

                     L is DOUBLE PRECISION array, dimension (N-1)
                      The (n-1) subdiagonal elements of the unit bidiagonal matrix
                      L, in elements 1 to N-1.

           D

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

           LD

                     LD is DOUBLE PRECISION array, dimension (N-1)
                      The n-1 elements L(i)*D(i).

           LLD

                     LLD is DOUBLE PRECISION array, dimension (N-1)
                      The n-1 elements L(i)*L(i)*D(i).

           PIVMIN

                     PIVMIN is DOUBLE PRECISION
                      The minimum pivot in the Sturm sequence.

           GAPTOL

                     GAPTOL is DOUBLE PRECISION
                      Tolerance that indicates when eigenvector entries are negligible
                      w.r.t. their contribution to the residual.

           Z

                     Z is DOUBLE PRECISION array, dimension (N)
                      On input, all entries of Z must be set to 0.
                      On output, Z contains the (scaled) r-th column of the
                      inverse. The scaling is such that Z(R) equals 1.

           WANTNC

                     WANTNC is LOGICAL
                      Specifies whether NEGCNT has to be computed.

           NEGCNT

                     NEGCNT is INTEGER
                      If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin
                      in the  matrix factorization L D L**T, and NEGCNT = -1 otherwise.

           ZTZ

                     ZTZ is DOUBLE PRECISION
                      The square of the 2-norm of Z.

           MINGMA

                     MINGMA is DOUBLE PRECISION
                      The reciprocal of the largest (in magnitude) diagonal
                      element of the inverse of L D L**T - sigma I.

           R

                     R is INTEGER
                      The twist index for the twisted factorization used to
                      compute Z.
                      On input, 0 <= R <= N. If R is input as 0, R is set to
                      the index where (L D L**T - sigma I)^{-1} is largest
                      in magnitude. If 1 <= R <= N, R is unchanged.
                      On output, R contains the twist index used to compute Z.
                      Ideally, R designates the position of the maximum entry in the
                      eigenvector.

           ISUPPZ

                     ISUPPZ is INTEGER array, dimension (2)
                      The support of the vector in Z, i.e., the vector Z is
                      nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ).

           NRMINV

                     NRMINV is DOUBLE PRECISION
                      NRMINV = 1/SQRT( ZTZ )

           RESID

                     RESID is DOUBLE PRECISION
                      The residual of the FP vector.
                      RESID = ABS( MINGMA )/SQRT( ZTZ )

           RQCORR

                     RQCORR is DOUBLE PRECISION
                      The Rayleigh Quotient correction to LAMBDA.
                      RQCORR = MINGMA*TMP

           WORK

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

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

       Contributors:
           Beresford Parlett, University of California, Berkeley, USA
            Jim Demmel, University of California, Berkeley, USA
            Inderjit Dhillon, University of Texas, Austin, USA
            Osni Marques, LBNL/NERSC, USA
            Christof Voemel, University of California, Berkeley, USA

   subroutine dlar2v (integer N, double precision, dimension( * ) X, double precision, dimension(
       * ) Y, double precision, dimension( * ) Z, integer INCX, double precision, dimension( * )
       C, double precision, dimension( * ) S, integer INCC)
       DLAR2V applies a vector of plane rotations with real cosines and real sines from both
       sides to a sequence of 2-by-2 symmetric/Hermitian matrices.

       Purpose:

            DLAR2V applies a vector of real plane rotations from both sides to
            a sequence of 2-by-2 real symmetric matrices, defined by the elements
            of the vectors x, y and z. For i = 1,2,...,n

               ( x(i)  z(i) ) := (  c(i)  s(i) ) ( x(i)  z(i) ) ( c(i) -s(i) )
               ( z(i)  y(i) )    ( -s(i)  c(i) ) ( z(i)  y(i) ) ( s(i)  c(i) )

       Parameters:
           N

                     N is INTEGER
                     The number of plane rotations to be applied.

           X

                     X is DOUBLE PRECISION array,
                                    dimension (1+(N-1)*INCX)
                     The vector x.

           Y

                     Y is DOUBLE PRECISION array,
                                    dimension (1+(N-1)*INCX)
                     The vector y.

           Z

                     Z is DOUBLE PRECISION array,
                                    dimension (1+(N-1)*INCX)
                     The vector z.

           INCX

                     INCX is INTEGER
                     The increment between elements of X, Y and Z. INCX > 0.

           C

                     C is DOUBLE PRECISION array, dimension (1+(N-1)*INCC)
                     The cosines of the plane rotations.

           S

                     S is DOUBLE PRECISION array, dimension (1+(N-1)*INCC)
                     The sines of the plane rotations.

           INCC

                     INCC is INTEGER
                     The increment between elements of C and S. INCC > 0.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlarf (character SIDE, integer M, integer N, double precision, dimension( * ) V,
       integer INCV, double precision TAU, double precision, dimension( ldc, * ) C, integer LDC,
       double precision, dimension( * ) WORK)
       DLARF applies an elementary reflector to a general rectangular matrix.

       Purpose:

            DLARF applies a real elementary reflector H to a real m by n matrix
            C, from either the left or the right. H is represented in the form

                  H = I - tau * v * v**T

            where tau is a real scalar and v is a real vector.

            If tau = 0, then H is taken to be the unit matrix.

       Parameters:
           SIDE

                     SIDE is CHARACTER*1
                     = 'L': form  H * C
                     = 'R': form  C * H

           M

                     M is INTEGER
                     The number of rows of the matrix C.

           N

                     N is INTEGER
                     The number of columns of the matrix C.

           V

                     V is DOUBLE PRECISION array, dimension
                                (1 + (M-1)*abs(INCV)) if SIDE = 'L'
                             or (1 + (N-1)*abs(INCV)) if SIDE = 'R'
                     The vector v in the representation of H. V is not used if
                     TAU = 0.

           INCV

                     INCV is INTEGER
                     The increment between elements of v. INCV <> 0.

           TAU

                     TAU is DOUBLE PRECISION
                     The value tau in the representation of H.

           C

                     C is DOUBLE PRECISION array, dimension (LDC,N)
                     On entry, the m by n matrix C.
                     On exit, C is overwritten by the matrix H * C if SIDE = 'L',
                     or C * H if SIDE = 'R'.

           LDC

                     LDC is INTEGER
                     The leading dimension of the array C. LDC >= max(1,M).

           WORK

                     WORK is DOUBLE PRECISION array, dimension
                                    (N) if SIDE = 'L'
                                 or (M) if SIDE = 'R'

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlarfb (character SIDE, character TRANS, character DIRECT, character STOREV,
       integer M, integer N, integer K, double precision, dimension( ldv, * ) V, integer LDV,
       double precision, dimension( ldt, * ) T, integer LDT, double precision, dimension( ldc, *
       ) C, integer LDC, double precision, dimension( ldwork, * ) WORK, integer LDWORK)
       DLARFB applies a block reflector or its transpose to a general rectangular matrix.

       Purpose:

            DLARFB applies a real block reflector H or its transpose H**T to a
            real m by n matrix C, from either the left or the right.

       Parameters:
           SIDE

                     SIDE is CHARACTER*1
                     = 'L': apply H or H**T from the Left
                     = 'R': apply H or H**T from the Right

           TRANS

                     TRANS is CHARACTER*1
                     = 'N': apply H (No transpose)
                     = 'T': apply H**T (Transpose)

           DIRECT

                     DIRECT is CHARACTER*1
                     Indicates how H is formed from a product of elementary
                     reflectors
                     = 'F': H = H(1) H(2) . . . H(k) (Forward)
                     = 'B': H = H(k) . . . H(2) H(1) (Backward)

           STOREV

                     STOREV is CHARACTER*1
                     Indicates how the vectors which define the elementary
                     reflectors are stored:
                     = 'C': Columnwise
                     = 'R': Rowwise

           M

                     M is INTEGER
                     The number of rows of the matrix C.

           N

                     N is INTEGER
                     The number of columns of the matrix C.

           K

                     K is INTEGER
                     The order of the matrix T (= the number of elementary
                     reflectors whose product defines the block reflector).

           V

                     V is DOUBLE PRECISION array, dimension
                                           (LDV,K) if STOREV = 'C'
                                           (LDV,M) if STOREV = 'R' and SIDE = 'L'
                                           (LDV,N) if STOREV = 'R' and SIDE = 'R'
                     The matrix V. See Further Details.

           LDV

                     LDV is INTEGER
                     The leading dimension of the array V.
                     If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
                     if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
                     if STOREV = 'R', LDV >= K.

           T

                     T is DOUBLE PRECISION array, dimension (LDT,K)
                     The triangular k by k matrix T in the representation of the
                     block reflector.

           LDT

                     LDT is INTEGER
                     The leading dimension of the array T. LDT >= K.

           C

                     C is DOUBLE PRECISION array, dimension (LDC,N)
                     On entry, the m by n matrix C.
                     On exit, C is overwritten by H*C or H**T*C or C*H or C*H**T.

           LDC

                     LDC is INTEGER
                     The leading dimension of the array C. LDC >= max(1,M).

           WORK

                     WORK is DOUBLE PRECISION array, dimension (LDWORK,K)

           LDWORK

                     LDWORK is INTEGER
                     The leading dimension of the array WORK.
                     If SIDE = 'L', LDWORK >= max(1,N);
                     if SIDE = 'R', LDWORK >= max(1,M).

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           June 2013

       Further Details:

             The shape of the matrix V and the storage of the vectors which define
             the H(i) is best illustrated by the following example with n = 5 and
             k = 3. The elements equal to 1 are not stored; the corresponding
             array elements are modified but restored on exit. The rest of the
             array is not used.

             DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':

                          V = (  1       )                 V = (  1 v1 v1 v1 v1 )
                              ( v1  1    )                     (     1 v2 v2 v2 )
                              ( v1 v2  1 )                     (        1 v3 v3 )
                              ( v1 v2 v3 )
                              ( v1 v2 v3 )

             DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':

                          V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
                              ( v1 v2 v3 )                     ( v2 v2 v2  1    )
                              (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
                              (     1 v3 )
                              (        1 )

   subroutine dlarfg (integer N, double precision ALPHA, double precision, dimension( * ) X,
       integer INCX, double precision TAU)
       DLARFG generates an elementary reflector (Householder matrix).

       Purpose:

            DLARFG generates a real elementary reflector H of order n, such
            that

                  H * ( alpha ) = ( beta ),   H**T * H = I.
                      (   x   )   (   0  )

            where alpha and beta are scalars, and x is an (n-1)-element real
            vector. H is represented in the form

                  H = I - tau * ( 1 ) * ( 1 v**T ) ,
                                ( v )

            where tau is a real scalar and v is a real (n-1)-element
            vector.

            If the elements of x are all zero, then tau = 0 and H is taken to be
            the unit matrix.

            Otherwise  1 <= tau <= 2.

       Parameters:
           N

                     N is INTEGER
                     The order of the elementary reflector.

           ALPHA

                     ALPHA is DOUBLE PRECISION
                     On entry, the value alpha.
                     On exit, it is overwritten with the value beta.

           X

                     X is DOUBLE PRECISION array, dimension
                                    (1+(N-2)*abs(INCX))
                     On entry, the vector x.
                     On exit, it is overwritten with the vector v.

           INCX

                     INCX is INTEGER
                     The increment between elements of X. INCX > 0.

           TAU

                     TAU is DOUBLE PRECISION
                     The value tau.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           November 2017

   subroutine dlarfgp (integer N, double precision ALPHA, double precision, dimension( * ) X,
       integer INCX, double precision TAU)
       DLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.

       Purpose:

            DLARFGP generates a real elementary reflector H of order n, such
            that

                  H * ( alpha ) = ( beta ),   H**T * H = I.
                      (   x   )   (   0  )

            where alpha and beta are scalars, beta is non-negative, and x is
            an (n-1)-element real vector.  H is represented in the form

                  H = I - tau * ( 1 ) * ( 1 v**T ) ,
                                ( v )

            where tau is a real scalar and v is a real (n-1)-element
            vector.

            If the elements of x are all zero, then tau = 0 and H is taken to be
            the unit matrix.

       Parameters:
           N

                     N is INTEGER
                     The order of the elementary reflector.

           ALPHA

                     ALPHA is DOUBLE PRECISION
                     On entry, the value alpha.
                     On exit, it is overwritten with the value beta.

           X

                     X is DOUBLE PRECISION array, dimension
                                    (1+(N-2)*abs(INCX))
                     On entry, the vector x.
                     On exit, it is overwritten with the vector v.

           INCX

                     INCX is INTEGER
                     The increment between elements of X. INCX > 0.

           TAU

                     TAU is DOUBLE PRECISION
                     The value tau.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           November 2017

   subroutine dlarft (character DIRECT, character STOREV, integer N, integer K, double precision,
       dimension( ldv, * ) V, integer LDV, double precision, dimension( * ) TAU, double
       precision, dimension( ldt, * ) T, integer LDT)
       DLARFT forms the triangular factor T of a block reflector H = I - vtvH

       Purpose:

            DLARFT forms the triangular factor T of a real block reflector H
            of order n, which is defined as a product of k elementary reflectors.

            If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;

            If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.

            If STOREV = 'C', the vector which defines the elementary reflector
            H(i) is stored in the i-th column of the array V, and

               H  =  I - V * T * V**T

            If STOREV = 'R', the vector which defines the elementary reflector
            H(i) is stored in the i-th row of the array V, and

               H  =  I - V**T * T * V

       Parameters:
           DIRECT

                     DIRECT is CHARACTER*1
                     Specifies the order in which the elementary reflectors are
                     multiplied to form the block reflector:
                     = 'F': H = H(1) H(2) . . . H(k) (Forward)
                     = 'B': H = H(k) . . . H(2) H(1) (Backward)

           STOREV

                     STOREV is CHARACTER*1
                     Specifies how the vectors which define the elementary
                     reflectors are stored (see also Further Details):
                     = 'C': columnwise
                     = 'R': rowwise

           N

                     N is INTEGER
                     The order of the block reflector H. N >= 0.

           K

                     K is INTEGER
                     The order of the triangular factor T (= the number of
                     elementary reflectors). K >= 1.

           V

                     V is DOUBLE PRECISION array, dimension
                                          (LDV,K) if STOREV = 'C'
                                          (LDV,N) if STOREV = 'R'
                     The matrix V. See further details.

           LDV

                     LDV is INTEGER
                     The leading dimension of the array V.
                     If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.

           TAU

                     TAU is DOUBLE PRECISION array, dimension (K)
                     TAU(i) must contain the scalar factor of the elementary
                     reflector H(i).

           T

                     T is DOUBLE PRECISION array, dimension (LDT,K)
                     The k by k triangular factor T of the block reflector.
                     If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
                     lower triangular. The rest of the array is not used.

           LDT

                     LDT is INTEGER
                     The leading dimension of the array T. LDT >= K.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

       Further Details:

             The shape of the matrix V and the storage of the vectors which define
             the H(i) is best illustrated by the following example with n = 5 and
             k = 3. The elements equal to 1 are not stored.

             DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':

                          V = (  1       )                 V = (  1 v1 v1 v1 v1 )
                              ( v1  1    )                     (     1 v2 v2 v2 )
                              ( v1 v2  1 )                     (        1 v3 v3 )
                              ( v1 v2 v3 )
                              ( v1 v2 v3 )

             DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':

                          V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
                              ( v1 v2 v3 )                     ( v2 v2 v2  1    )
                              (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
                              (     1 v3 )
                              (        1 )

   subroutine dlarfx (character SIDE, integer M, integer N, double precision, dimension( * ) V,
       double precision TAU, double precision, dimension( ldc, * ) C, integer LDC, double
       precision, dimension( * ) WORK)
       DLARFX applies an elementary reflector to a general rectangular matrix, with loop
       unrolling when the reflector has order ≤ 10.

       Purpose:

            DLARFX applies a real elementary reflector H to a real m by n
            matrix C, from either the left or the right. H is represented in the
            form

                  H = I - tau * v * v**T

            where tau is a real scalar and v is a real vector.

            If tau = 0, then H is taken to be the unit matrix

            This version uses inline code if H has order < 11.

       Parameters:
           SIDE

                     SIDE is CHARACTER*1
                     = 'L': form  H * C
                     = 'R': form  C * H

           M

                     M is INTEGER
                     The number of rows of the matrix C.

           N

                     N is INTEGER
                     The number of columns of the matrix C.

           V

                     V is DOUBLE PRECISION array, dimension (M) if SIDE = 'L'
                                                or (N) if SIDE = 'R'
                     The vector v in the representation of H.

           TAU

                     TAU is DOUBLE PRECISION
                     The value tau in the representation of H.

           C

                     C is DOUBLE PRECISION array, dimension (LDC,N)
                     On entry, the m by n matrix C.
                     On exit, C is overwritten by the matrix H * C if SIDE = 'L',
                     or C * H if SIDE = 'R'.

           LDC

                     LDC is INTEGER
                     The leading dimension of the array C. LDA >= (1,M).

           WORK

                     WORK is DOUBLE PRECISION array, dimension
                                 (N) if SIDE = 'L'
                                 or (M) if SIDE = 'R'
                     WORK is not referenced if H has order < 11.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlargv (integer N, double precision, dimension( * ) X, integer INCX, double
       precision, dimension( * ) Y, integer INCY, double precision, dimension( * ) C, integer
       INCC)
       DLARGV generates a vector of plane rotations with real cosines and real sines.

       Purpose:

            DLARGV generates a vector of real plane rotations, determined by
            elements of the real vectors x and y. For i = 1,2,...,n

               (  c(i)  s(i) ) ( x(i) ) = ( a(i) )
               ( -s(i)  c(i) ) ( y(i) ) = (   0  )

       Parameters:
           N

                     N is INTEGER
                     The number of plane rotations to be generated.

           X

                     X is DOUBLE PRECISION array,
                                    dimension (1+(N-1)*INCX)
                     On entry, the vector x.
                     On exit, x(i) is overwritten by a(i), for i = 1,...,n.

           INCX

                     INCX is INTEGER
                     The increment between elements of X. INCX > 0.

           Y

                     Y is DOUBLE PRECISION array,
                                    dimension (1+(N-1)*INCY)
                     On entry, the vector y.
                     On exit, the sines of the plane rotations.

           INCY

                     INCY is INTEGER
                     The increment between elements of Y. INCY > 0.

           C

                     C is DOUBLE PRECISION array, dimension (1+(N-1)*INCC)
                     The cosines of the plane rotations.

           INCC

                     INCC is INTEGER
                     The increment between elements of C. INCC > 0.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlarrv (integer N, double precision VL, double precision VU, double precision,
       dimension( * ) D, double precision, dimension( * ) L, double precision PIVMIN, integer,
       dimension( * ) ISPLIT, integer M, integer DOL, integer DOU, double precision MINRGP,
       double precision RTOL1, double precision RTOL2, double precision, dimension( * ) W, double
       precision, dimension( * ) WERR, double precision, dimension( * ) WGAP, integer, dimension(
       * ) IBLOCK, integer, dimension( * ) INDEXW, double precision, dimension( * ) GERS, double
       precision, dimension( ldz, * ) Z, integer LDZ, integer, dimension( * ) ISUPPZ, double
       precision, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)
       DLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the
       eigenvalues of L D LT.

       Purpose:

            DLARRV computes the eigenvectors of the tridiagonal matrix
            T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T.
            The input eigenvalues should have been computed by DLARRE.

       Parameters:
           N

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

           VL

                     VL is DOUBLE PRECISION
                     Lower bound of the interval that contains the desired
                     eigenvalues. VL < VU. Needed to compute gaps on the left or right
                     end of the extremal eigenvalues in the desired RANGE.

           VU

                     VU is DOUBLE PRECISION
                     Upper bound of the interval that contains the desired
                     eigenvalues. VL < VU.
                     Note: VU is currently not used by this implementation of DLARRV, VU is
                     passed to DLARRV because it could be used compute gaps on the right end
                     of the extremal eigenvalues. However, with not much initial accuracy in
                     LAMBDA and VU, the formula can lead to an overestimation of the right gap
                     and thus to inadequately early RQI 'convergence'. This is currently
                     prevented this by forcing a small right gap. And so it turns out that VU
                     is currently not used by this implementation of DLARRV.

           D

                     D is DOUBLE PRECISION array, dimension (N)
                     On entry, the N diagonal elements of the diagonal matrix D.
                     On exit, D may be overwritten.

           L

                     L is DOUBLE PRECISION array, dimension (N)
                     On entry, the (N-1) subdiagonal elements of the unit
                     bidiagonal matrix L are in elements 1 to N-1 of L
                     (if the matrix is not split.) At the end of each block
                     is stored the corresponding shift as given by DLARRE.
                     On exit, L is overwritten.

           PIVMIN

                     PIVMIN is DOUBLE PRECISION
                     The minimum pivot allowed in the Sturm sequence.

           ISPLIT

                     ISPLIT is INTEGER array, dimension (N)
                     The splitting points, at which T breaks up into blocks.
                     The first block consists of rows/columns 1 to
                     ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
                     through ISPLIT( 2 ), etc.

           M

                     M is INTEGER
                     The total number of input eigenvalues.  0 <= M <= N.

           DOL

                     DOL is INTEGER

           DOU

                     DOU is INTEGER
                     If the user wants to compute only selected eigenvectors from all
                     the eigenvalues supplied, he can specify an index range DOL:DOU.
                     Or else the setting DOL=1, DOU=M should be applied.
                     Note that DOL and DOU refer to the order in which the eigenvalues
                     are stored in W.
                     If the user wants to compute only selected eigenpairs, then
                     the columns DOL-1 to DOU+1 of the eigenvector space Z contain the
                     computed eigenvectors. All other columns of Z are set to zero.

           MINRGP

                     MINRGP is DOUBLE PRECISION

           RTOL1

                     RTOL1 is DOUBLE PRECISION

           RTOL2

                     RTOL2 is DOUBLE PRECISION
                      Parameters for bisection.
                      An interval [LEFT,RIGHT] has converged if
                      RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )

           W

                     W is DOUBLE PRECISION array, dimension (N)
                     The first M elements of W contain the APPROXIMATE eigenvalues for
                     which eigenvectors are to be computed.  The eigenvalues
                     should be grouped by split-off block and ordered from
                     smallest to largest within the block ( The output array
                     W from DLARRE is expected here ). Furthermore, they are with
                     respect to the shift of the corresponding root representation
                     for their block. On exit, W holds the eigenvalues of the
                     UNshifted matrix.

           WERR

                     WERR is DOUBLE PRECISION array, dimension (N)
                     The first M elements contain the semiwidth of the uncertainty
                     interval of the corresponding eigenvalue in W

           WGAP

                     WGAP is DOUBLE PRECISION array, dimension (N)
                     The separation from the right neighbor eigenvalue in W.

           IBLOCK

                     IBLOCK is INTEGER array, dimension (N)
                     The indices of the blocks (submatrices) associated with the
                     corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
                     W(i) belongs to the first block from the top, =2 if W(i)
                     belongs to the second block, etc.

           INDEXW

                     INDEXW is INTEGER array, dimension (N)
                     The indices of the eigenvalues within each block (submatrix);
                     for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
                     i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.

           GERS

                     GERS is DOUBLE PRECISION array, dimension (2*N)
                     The N Gerschgorin intervals (the i-th Gerschgorin interval
                     is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should
                     be computed from the original UNshifted matrix.

           Z

                     Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
                     If INFO = 0, the first M columns of Z contain the
                     orthonormal eigenvectors of the matrix T
                     corresponding to the input eigenvalues, with the i-th
                     column of Z holding the eigenvector associated with W(i).
                     Note: the user must ensure that at least max(1,M) columns are
                     supplied in the array Z.

           LDZ

                     LDZ is INTEGER
                     The leading dimension of the array Z.  LDZ >= 1, and if
                     JOBZ = 'V', LDZ >= max(1,N).

           ISUPPZ

                     ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
                     The support of the eigenvectors in Z, i.e., the indices
                     indicating the nonzero elements in Z. The I-th eigenvector
                     is nonzero only in elements ISUPPZ( 2*I-1 ) through
                     ISUPPZ( 2*I ).

           WORK

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

           IWORK

                     IWORK is INTEGER array, dimension (7*N)

           INFO

                     INFO is INTEGER
                     = 0:  successful exit

                     > 0:  A problem occurred in DLARRV.
                     < 0:  One of the called subroutines signaled an internal problem.
                           Needs inspection of the corresponding parameter IINFO
                           for further information.

                     =-1:  Problem in DLARRB when refining a child's eigenvalues.
                     =-2:  Problem in DLARRF when computing the RRR of a child.
                           When a child is inside a tight cluster, it can be difficult
                           to find an RRR. A partial remedy from the user's point of
                           view is to make the parameter MINRGP smaller and recompile.
                           However, as the orthogonality of the computed vectors is
                           proportional to 1/MINRGP, the user should be aware that
                           he might be trading in precision when he decreases MINRGP.
                     =-3:  Problem in DLARRB when refining a single eigenvalue
                           after the Rayleigh correction was rejected.
                     = 5:  The Rayleigh Quotient Iteration failed to converge to
                           full accuracy in MAXITR steps.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           June 2016

       Contributors:
           Beresford Parlett, University of California, Berkeley, USA
            Jim Demmel, University of California, Berkeley, USA
            Inderjit Dhillon, University of Texas, Austin, USA
            Osni Marques, LBNL/NERSC, USA
            Christof Voemel, University of California, Berkeley, USA

   subroutine dlartv (integer N, double precision, dimension( * ) X, integer INCX, double
       precision, dimension( * ) Y, integer INCY, double precision, dimension( * ) C, double
       precision, dimension( * ) S, integer INCC)
       DLARTV applies a vector of plane rotations with real cosines and real sines to the
       elements of a pair of vectors.

       Purpose:

            DLARTV applies a vector of real plane rotations to elements of the
            real vectors x and y. For i = 1,2,...,n

               ( x(i) ) := (  c(i)  s(i) ) ( x(i) )
               ( y(i) )    ( -s(i)  c(i) ) ( y(i) )

       Parameters:
           N

                     N is INTEGER
                     The number of plane rotations to be applied.

           X

                     X is DOUBLE PRECISION array,
                                    dimension (1+(N-1)*INCX)
                     The vector x.

           INCX

                     INCX is INTEGER
                     The increment between elements of X. INCX > 0.

           Y

                     Y is DOUBLE PRECISION array,
                                    dimension (1+(N-1)*INCY)
                     The vector y.

           INCY

                     INCY is INTEGER
                     The increment between elements of Y. INCY > 0.

           C

                     C is DOUBLE PRECISION array, dimension (1+(N-1)*INCC)
                     The cosines of the plane rotations.

           S

                     S is DOUBLE PRECISION array, dimension (1+(N-1)*INCC)
                     The sines of the plane rotations.

           INCC

                     INCC is INTEGER
                     The increment between elements of C and S. INCC > 0.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlaswp (integer N, double precision, dimension( lda, * ) A, integer LDA, integer
       K1, integer K2, integer, dimension( * ) IPIV, integer INCX)
       DLASWP performs a series of row interchanges on a general rectangular matrix.

       Purpose:

            DLASWP performs a series of row interchanges on the matrix A.
            One row interchange is initiated for each of rows K1 through K2 of A.

       Parameters:
           N

                     N is INTEGER
                     The number of columns of the matrix A.

           A

                     A is DOUBLE PRECISION array, dimension (LDA,N)
                     On entry, the matrix of column dimension N to which the row
                     interchanges will be applied.
                     On exit, the permuted matrix.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.

           K1

                     K1 is INTEGER
                     The first element of IPIV for which a row interchange will
                     be done.

           K2

                     K2 is INTEGER
                     (K2-K1+1) is the number of elements of IPIV for which a row
                     interchange will be done.

           IPIV

                     IPIV is INTEGER array, dimension (K1+(K2-K1)*abs(INCX))
                     The vector of pivot indices. Only the elements in positions
                     K1 through K1+(K2-K1)*abs(INCX) of IPIV are accessed.
                     IPIV(K1+(K-K1)*abs(INCX)) = L implies rows K and L are to be
                     interchanged.

           INCX

                     INCX is INTEGER
                     The increment between successive values of IPIV. If INCX
                     is negative, the pivots are applied in reverse order.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           June 2017

       Further Details:

             Modified by
              R. C. Whaley, Computer Science Dept., Univ. of Tenn., Knoxville, USA

   subroutine dlat2s (character UPLO, integer N, double precision, dimension( lda, * ) A, integer
       LDA, real, dimension( ldsa, * ) SA, integer LDSA, integer INFO)
       DLAT2S converts a double-precision triangular matrix to a single-precision triangular
       matrix.

       Purpose:

            DLAT2S converts a DOUBLE PRECISION triangular matrix, SA, to a SINGLE
            PRECISION triangular matrix, A.

            RMAX is the overflow for the SINGLE PRECISION arithmetic
            DLAS2S checks that all the entries of A are between -RMAX and
            RMAX. If not the conversion is aborted and a flag is raised.

            This is an auxiliary routine so there is no argument checking.

       Parameters:
           UPLO

                     UPLO is CHARACTER*1
                     = 'U':  A is upper triangular;
                     = 'L':  A is lower triangular.

           N

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

           A

                     A is DOUBLE PRECISION array, dimension (LDA,N)
                     On entry, the N-by-N triangular coefficient matrix A.

           LDA

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

           SA

                     SA is REAL array, dimension (LDSA,N)
                     Only the UPLO part of SA is referenced.  On exit, if INFO=0,
                     the N-by-N coefficient matrix SA; if INFO>0, the content of
                     the UPLO part of SA is unspecified.

           LDSA

                     LDSA is INTEGER
                     The leading dimension of the array SA.  LDSA >= max(1,M).

           INFO

                     INFO is INTEGER
                     = 0:  successful exit.
                     = 1:  an entry of the matrix A is greater than the SINGLE
                           PRECISION overflow threshold, in this case, the content
                           of the UPLO part of SA in exit is unspecified.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlatbs (character UPLO, character TRANS, character DIAG, character NORMIN, integer
       N, integer KD, double precision, dimension( ldab, * ) AB, integer LDAB, double precision,
       dimension( * ) X, double precision SCALE, double precision, dimension( * ) CNORM, integer
       INFO)
       DLATBS solves a triangular banded system of equations.

       Purpose:

            DLATBS solves one of the triangular systems

               A *x = s*b  or  A**T*x = s*b

            with scaling to prevent overflow, where A is an upper or lower
            triangular band matrix.  Here A**T denotes the transpose of A, x and b
            are n-element vectors, and s is a scaling factor, usually less than
            or equal to 1, chosen so that the components of x will be less than
            the overflow threshold.  If the unscaled problem will not cause
            overflow, the Level 2 BLAS routine DTBSV is called.  If the matrix A
            is singular (A(j,j) = 0 for some j), then s is set to 0 and a
            non-trivial solution to A*x = 0 is returned.

       Parameters:
           UPLO

                     UPLO is CHARACTER*1
                     Specifies whether the matrix A is upper or lower triangular.
                     = 'U':  Upper triangular
                     = 'L':  Lower triangular

           TRANS

                     TRANS is CHARACTER*1
                     Specifies the operation applied to A.
                     = 'N':  Solve A * x = s*b  (No transpose)
                     = 'T':  Solve A**T* x = s*b  (Transpose)
                     = 'C':  Solve A**T* x = s*b  (Conjugate transpose = Transpose)

           DIAG

                     DIAG is CHARACTER*1
                     Specifies whether or not the matrix A is unit triangular.
                     = 'N':  Non-unit triangular
                     = 'U':  Unit triangular

           NORMIN

                     NORMIN is CHARACTER*1
                     Specifies whether CNORM has been set or not.
                     = 'Y':  CNORM contains the column norms on entry
                     = 'N':  CNORM is not set on entry.  On exit, the norms will
                             be computed and stored in CNORM.

           N

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

           KD

                     KD is INTEGER
                     The number of subdiagonals or superdiagonals in the
                     triangular matrix A.  KD >= 0.

           AB

                     AB is DOUBLE PRECISION array, dimension (LDAB,N)
                     The upper or lower triangular band matrix A, stored in the
                     first KD+1 rows of the array. The j-th column of A is stored
                     in the j-th column of the array AB as follows:
                     if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
                     if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).

           LDAB

                     LDAB is INTEGER
                     The leading dimension of the array AB.  LDAB >= KD+1.

           X

                     X is DOUBLE PRECISION array, dimension (N)
                     On entry, the right hand side b of the triangular system.
                     On exit, X is overwritten by the solution vector x.

           SCALE

                     SCALE is DOUBLE PRECISION
                     The scaling factor s for the triangular system
                        A * x = s*b  or  A**T* x = s*b.
                     If SCALE = 0, the matrix A is singular or badly scaled, and
                     the vector x is an exact or approximate solution to A*x = 0.

           CNORM

                     CNORM is DOUBLE PRECISION array, dimension (N)

                     If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
                     contains the norm of the off-diagonal part of the j-th column
                     of A.  If TRANS = 'N', CNORM(j) must be greater than or equal
                     to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
                     must be greater than or equal to the 1-norm.

                     If NORMIN = 'N', CNORM is an output argument and CNORM(j)
                     returns the 1-norm of the offdiagonal part of the j-th column
                     of A.

           INFO

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

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

       Further Details:

             A rough bound on x is computed; if that is less than overflow, DTBSV
             is called, otherwise, specific code is used which checks for possible
             overflow or divide-by-zero at every operation.

             A columnwise scheme is used for solving A*x = b.  The basic algorithm
             if A is lower triangular is

                  x[1:n] := b[1:n]
                  for j = 1, ..., n
                       x(j) := x(j) / A(j,j)
                       x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
                  end

             Define bounds on the components of x after j iterations of the loop:
                M(j) = bound on x[1:j]
                G(j) = bound on x[j+1:n]
             Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.

             Then for iteration j+1 we have
                M(j+1) <= G(j) / | A(j+1,j+1) |
                G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
                       <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )

             where CNORM(j+1) is greater than or equal to the infinity-norm of
             column j+1 of A, not counting the diagonal.  Hence

                G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
                             1<=i<=j
             and

                |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
                                              1<=i< j

             Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTBSV if the
             reciprocal of the largest M(j), j=1,..,n, is larger than
             max(underflow, 1/overflow).

             The bound on x(j) is also used to determine when a step in the
             columnwise method can be performed without fear of overflow.  If
             the computed bound is greater than a large constant, x is scaled to
             prevent overflow, but if the bound overflows, x is set to 0, x(j) to
             1, and scale to 0, and a non-trivial solution to A*x = 0 is found.

             Similarly, a row-wise scheme is used to solve A**T*x = b.  The basic
             algorithm for A upper triangular is

                  for j = 1, ..., n
                       x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j)
                  end

             We simultaneously compute two bounds
                  G(j) = bound on ( b(i) - A[1:i-1,i]**T * x[1:i-1] ), 1<=i<=j
                  M(j) = bound on x(i), 1<=i<=j

             The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
             add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
             Then the bound on x(j) is

                  M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |

                       <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
                                 1<=i<=j

             and we can safely call DTBSV if 1/M(n) and 1/G(n) are both greater
             than max(underflow, 1/overflow).

   subroutine dlatdf (integer IJOB, integer N, double precision, dimension( ldz, * ) Z, integer
       LDZ, double precision, dimension( * ) RHS, double precision RDSUM, double precision
       RDSCAL, integer, dimension( * ) IPIV, integer, dimension( * ) JPIV)
       DLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a
       contribution to the reciprocal Dif-estimate.

       Purpose:

            DLATDF uses the LU factorization of the n-by-n matrix Z computed by
            DGETC2 and computes a contribution to the reciprocal Dif-estimate
            by solving Z * x = b for x, and choosing the r.h.s. b such that
            the norm of x is as large as possible. On entry RHS = b holds the
            contribution from earlier solved sub-systems, and on return RHS = x.

            The factorization of Z returned by DGETC2 has the form Z = P*L*U*Q,
            where P and Q are permutation matrices. L is lower triangular with
            unit diagonal elements and U is upper triangular.

       Parameters:
           IJOB

                     IJOB is INTEGER
                     IJOB = 2: First compute an approximative null-vector e
                         of Z using DGECON, e is normalized and solve for
                         Zx = +-e - f with the sign giving the greater value
                         of 2-norm(x). About 5 times as expensive as Default.
                     IJOB .ne. 2: Local look ahead strategy where all entries of
                         the r.h.s. b is chosen as either +1 or -1 (Default).

           N

                     N is INTEGER
                     The number of columns of the matrix Z.

           Z

                     Z is DOUBLE PRECISION array, dimension (LDZ, N)
                     On entry, the LU part of the factorization of the n-by-n
                     matrix Z computed by DGETC2:  Z = P * L * U * Q

           LDZ

                     LDZ is INTEGER
                     The leading dimension of the array Z.  LDA >= max(1, N).

           RHS

                     RHS is DOUBLE PRECISION array, dimension (N)
                     On entry, RHS contains contributions from other subsystems.
                     On exit, RHS contains the solution of the subsystem with
                     entries acoording to the value of IJOB (see above).

           RDSUM

                     RDSUM is DOUBLE PRECISION
                     On entry, the sum of squares of computed contributions to
                     the Dif-estimate under computation by DTGSYL, where the
                     scaling factor RDSCAL (see below) has been factored out.
                     On exit, the corresponding sum of squares updated with the
                     contributions from the current sub-system.
                     If TRANS = 'T' RDSUM is not touched.
                     NOTE: RDSUM only makes sense when DTGSY2 is called by STGSYL.

           RDSCAL

                     RDSCAL is DOUBLE PRECISION
                     On entry, scaling factor used to prevent overflow in RDSUM.
                     On exit, RDSCAL is updated w.r.t. the current contributions
                     in RDSUM.
                     If TRANS = 'T', RDSCAL is not touched.
                     NOTE: RDSCAL only makes sense when DTGSY2 is called by
                           DTGSYL.

           IPIV

                     IPIV is INTEGER array, dimension (N).
                     The pivot indices; for 1 <= i <= N, row i of the
                     matrix has been interchanged with row IPIV(i).

           JPIV

                     JPIV is INTEGER array, dimension (N).
                     The pivot indices; for 1 <= j <= N, column j of the
                     matrix has been interchanged with column JPIV(j).

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           June 2016

       Further Details:
           This routine is a further developed implementation of algorithm BSOLVE in [1] using
           complete pivoting in the LU factorization.

       Contributors:
           Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901
           87 Umea, Sweden.

       References:

             [1] Bo Kagstrom and Lars Westin,
                 Generalized Schur Methods with Condition Estimators for
                 Solving the Generalized Sylvester Equation, IEEE Transactions
                 on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751.

             [2] Peter Poromaa,
                 On Efficient and Robust Estimators for the Separation
                 between two Regular Matrix Pairs with Applications in
                 Condition Estimation. Report IMINF-95.05, Departement of
                 Computing Science, Umea University, S-901 87 Umea, Sweden, 1995.

   subroutine dlatps (character UPLO, character TRANS, character DIAG, character NORMIN, integer
       N, double precision, dimension( * ) AP, double precision, dimension( * ) X, double
       precision SCALE, double precision, dimension( * ) CNORM, integer INFO)
       DLATPS solves a triangular system of equations with the matrix held in packed storage.

       Purpose:

            DLATPS solves one of the triangular systems

               A *x = s*b  or  A**T*x = s*b

            with scaling to prevent overflow, where A is an upper or lower
            triangular matrix stored in packed form.  Here A**T denotes the
            transpose of A, x and b are n-element vectors, and s is a scaling
            factor, usually less than or equal to 1, chosen so that the
            components of x will be less than the overflow threshold.  If the
            unscaled problem will not cause overflow, the Level 2 BLAS routine
            DTPSV is called. If the matrix A is singular (A(j,j) = 0 for some j),
            then s is set to 0 and a non-trivial solution to A*x = 0 is returned.

       Parameters:
           UPLO

                     UPLO is CHARACTER*1
                     Specifies whether the matrix A is upper or lower triangular.
                     = 'U':  Upper triangular
                     = 'L':  Lower triangular

           TRANS

                     TRANS is CHARACTER*1
                     Specifies the operation applied to A.
                     = 'N':  Solve A * x = s*b  (No transpose)
                     = 'T':  Solve A**T* x = s*b  (Transpose)
                     = 'C':  Solve A**T* x = s*b  (Conjugate transpose = Transpose)

           DIAG

                     DIAG is CHARACTER*1
                     Specifies whether or not the matrix A is unit triangular.
                     = 'N':  Non-unit triangular
                     = 'U':  Unit triangular

           NORMIN

                     NORMIN is CHARACTER*1
                     Specifies whether CNORM has been set or not.
                     = 'Y':  CNORM contains the column norms on entry
                     = 'N':  CNORM is not set on entry.  On exit, the norms will
                             be computed and stored in CNORM.

           N

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

           AP

                     AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
                     The upper or lower triangular matrix A, packed columnwise in
                     a linear array.  The j-th column of A is stored in the array
                     AP as follows:
                     if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
                     if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.

           X

                     X is DOUBLE PRECISION array, dimension (N)
                     On entry, the right hand side b of the triangular system.
                     On exit, X is overwritten by the solution vector x.

           SCALE

                     SCALE is DOUBLE PRECISION
                     The scaling factor s for the triangular system
                        A * x = s*b  or  A**T* x = s*b.
                     If SCALE = 0, the matrix A is singular or badly scaled, and
                     the vector x is an exact or approximate solution to A*x = 0.

           CNORM

                     CNORM is DOUBLE PRECISION array, dimension (N)

                     If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
                     contains the norm of the off-diagonal part of the j-th column
                     of A.  If TRANS = 'N', CNORM(j) must be greater than or equal
                     to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
                     must be greater than or equal to the 1-norm.

                     If NORMIN = 'N', CNORM is an output argument and CNORM(j)
                     returns the 1-norm of the offdiagonal part of the j-th column
                     of A.

           INFO

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

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

       Further Details:

             A rough bound on x is computed; if that is less than overflow, DTPSV
             is called, otherwise, specific code is used which checks for possible
             overflow or divide-by-zero at every operation.

             A columnwise scheme is used for solving A*x = b.  The basic algorithm
             if A is lower triangular is

                  x[1:n] := b[1:n]
                  for j = 1, ..., n
                       x(j) := x(j) / A(j,j)
                       x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
                  end

             Define bounds on the components of x after j iterations of the loop:
                M(j) = bound on x[1:j]
                G(j) = bound on x[j+1:n]
             Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.

             Then for iteration j+1 we have
                M(j+1) <= G(j) / | A(j+1,j+1) |
                G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
                       <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )

             where CNORM(j+1) is greater than or equal to the infinity-norm of
             column j+1 of A, not counting the diagonal.  Hence

                G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
                             1<=i<=j
             and

                |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
                                              1<=i< j

             Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTPSV if the
             reciprocal of the largest M(j), j=1,..,n, is larger than
             max(underflow, 1/overflow).

             The bound on x(j) is also used to determine when a step in the
             columnwise method can be performed without fear of overflow.  If
             the computed bound is greater than a large constant, x is scaled to
             prevent overflow, but if the bound overflows, x is set to 0, x(j) to
             1, and scale to 0, and a non-trivial solution to A*x = 0 is found.

             Similarly, a row-wise scheme is used to solve A**T*x = b.  The basic
             algorithm for A upper triangular is

                  for j = 1, ..., n
                       x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j)
                  end

             We simultaneously compute two bounds
                  G(j) = bound on ( b(i) - A[1:i-1,i]**T * x[1:i-1] ), 1<=i<=j
                  M(j) = bound on x(i), 1<=i<=j

             The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
             add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
             Then the bound on x(j) is

                  M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |

                       <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
                                 1<=i<=j

             and we can safely call DTPSV if 1/M(n) and 1/G(n) are both greater
             than max(underflow, 1/overflow).

   subroutine dlatrd (character UPLO, integer N, integer NB, double precision, dimension( lda, *
       ) A, integer LDA, double precision, dimension( * ) E, double precision, dimension( * )
       TAU, double precision, dimension( ldw, * ) W, integer LDW)
       DLATRD reduces the first nb rows and columns of a symmetric/Hermitian matrix A to real
       tridiagonal form by an orthogonal similarity transformation.

       Purpose:

            DLATRD reduces NB rows and columns of a real symmetric matrix A to
            symmetric tridiagonal form by an orthogonal similarity
            transformation Q**T * A * Q, and returns the matrices V and W which are
            needed to apply the transformation to the unreduced part of A.

            If UPLO = 'U', DLATRD reduces the last NB rows and columns of a
            matrix, of which the upper triangle is supplied;
            if UPLO = 'L', DLATRD reduces the first NB rows and columns of a
            matrix, of which the lower triangle is supplied.

            This is an auxiliary routine called by DSYTRD.

       Parameters:
           UPLO

                     UPLO is CHARACTER*1
                     Specifies whether the upper or lower triangular part of the
                     symmetric matrix A is stored:
                     = 'U': Upper triangular
                     = 'L': Lower triangular

           N

                     N is INTEGER
                     The order of the matrix A.

           NB

                     NB is INTEGER
                     The number of rows and columns to be reduced.

           A

                     A is DOUBLE PRECISION array, dimension (LDA,N)
                     On entry, the symmetric matrix A.  If UPLO = 'U', the leading
                     n-by-n upper triangular part of A contains the upper
                     triangular part of the matrix A, and the strictly lower
                     triangular part of A is not referenced.  If UPLO = 'L', the
                     leading n-by-n lower triangular part of A contains the lower
                     triangular part of the matrix A, and the strictly upper
                     triangular part of A is not referenced.
                     On exit:
                     if UPLO = 'U', the last NB columns have been reduced to
                       tridiagonal form, with the diagonal elements overwriting
                       the diagonal elements of A; the elements above the diagonal
                       with the array TAU, represent the orthogonal matrix Q as a
                       product of elementary reflectors;
                     if UPLO = 'L', the first NB columns have been reduced to
                       tridiagonal form, with the diagonal elements overwriting
                       the diagonal elements of A; the elements below the diagonal
                       with the array TAU, represent the  orthogonal matrix Q as a
                       product of elementary reflectors.
                     See Further Details.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= (1,N).

           E

                     E is DOUBLE PRECISION array, dimension (N-1)
                     If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal
                     elements of the last NB columns of the reduced matrix;
                     if UPLO = 'L', E(1:nb) contains the subdiagonal elements of
                     the first NB columns of the reduced matrix.

           TAU

                     TAU is DOUBLE PRECISION array, dimension (N-1)
                     The scalar factors of the elementary reflectors, stored in
                     TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.
                     See Further Details.

           W

                     W is DOUBLE PRECISION array, dimension (LDW,NB)
                     The n-by-nb matrix W required to update the unreduced part
                     of A.

           LDW

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

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

       Further Details:

             If UPLO = 'U', the matrix Q is represented as a product of elementary
             reflectors

                Q = H(n) H(n-1) . . . H(n-nb+1).

             Each H(i) has the form

                H(i) = I - tau * v * v**T

             where tau is a real scalar, and v is a real vector with
             v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i),
             and tau in TAU(i-1).

             If UPLO = 'L', the matrix Q is represented as a product of elementary
             reflectors

                Q = H(1) H(2) . . . H(nb).

             Each H(i) has the form

                H(i) = I - tau * v * v**T

             where tau is a real scalar, and v is a real vector with
             v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
             and tau in TAU(i).

             The elements of the vectors v together form the n-by-nb matrix V
             which is needed, with W, to apply the transformation to the unreduced
             part of the matrix, using a symmetric rank-2k update of the form:
             A := A - V*W**T - W*V**T.

             The contents of A on exit are illustrated by the following examples
             with n = 5 and nb = 2:

             if UPLO = 'U':                       if UPLO = 'L':

               (  a   a   a   v4  v5 )              (  d                  )
               (      a   a   v4  v5 )              (  1   d              )
               (          a   1   v5 )              (  v1  1   a          )
               (              d   1  )              (  v1  v2  a   a      )
               (                  d  )              (  v1  v2  a   a   a  )

             where d denotes a diagonal element of the reduced matrix, a denotes
             an element of the original matrix that is unchanged, and vi denotes
             an element of the vector defining H(i).

   subroutine dlatrs (character UPLO, character TRANS, character DIAG, character NORMIN, integer
       N, double precision, dimension( lda, * ) A, integer LDA, double precision, dimension( * )
       X, double precision SCALE, double precision, dimension( * ) CNORM, integer INFO)
       DLATRS solves a triangular system of equations with the scale factor set to prevent
       overflow.

       Purpose:

            DLATRS solves one of the triangular systems

               A *x = s*b  or  A**T *x = s*b

            with scaling to prevent overflow.  Here A is an upper or lower
            triangular matrix, A**T denotes the transpose of A, x and b are
            n-element vectors, and s is a scaling factor, usually less than
            or equal to 1, chosen so that the components of x will be less than
            the overflow threshold.  If the unscaled problem will not cause
            overflow, the Level 2 BLAS routine DTRSV is called.  If the matrix A
            is singular (A(j,j) = 0 for some j), then s is set to 0 and a
            non-trivial solution to A*x = 0 is returned.

       Parameters:
           UPLO

                     UPLO is CHARACTER*1
                     Specifies whether the matrix A is upper or lower triangular.
                     = 'U':  Upper triangular
                     = 'L':  Lower triangular

           TRANS

                     TRANS is CHARACTER*1
                     Specifies the operation applied to A.
                     = 'N':  Solve A * x = s*b  (No transpose)
                     = 'T':  Solve A**T* x = s*b  (Transpose)
                     = 'C':  Solve A**T* x = s*b  (Conjugate transpose = Transpose)

           DIAG

                     DIAG is CHARACTER*1
                     Specifies whether or not the matrix A is unit triangular.
                     = 'N':  Non-unit triangular
                     = 'U':  Unit triangular

           NORMIN

                     NORMIN is CHARACTER*1
                     Specifies whether CNORM has been set or not.
                     = 'Y':  CNORM contains the column norms on entry
                     = 'N':  CNORM is not set on entry.  On exit, the norms will
                             be computed and stored in CNORM.

           N

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

           A

                     A is DOUBLE PRECISION array, dimension (LDA,N)
                     The triangular matrix A.  If UPLO = 'U', the leading n by n
                     upper triangular part of the array A contains the upper
                     triangular matrix, and the strictly lower triangular part of
                     A is not referenced.  If UPLO = 'L', the leading n by n lower
                     triangular part of the array A contains the lower triangular
                     matrix, and the strictly upper triangular part of A is not
                     referenced.  If DIAG = 'U', the diagonal elements of A are
                     also not referenced and are assumed to be 1.

           LDA

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

           X

                     X is DOUBLE PRECISION array, dimension (N)
                     On entry, the right hand side b of the triangular system.
                     On exit, X is overwritten by the solution vector x.

           SCALE

                     SCALE is DOUBLE PRECISION
                     The scaling factor s for the triangular system
                        A * x = s*b  or  A**T* x = s*b.
                     If SCALE = 0, the matrix A is singular or badly scaled, and
                     the vector x is an exact or approximate solution to A*x = 0.

           CNORM

                     CNORM is DOUBLE PRECISION array, dimension (N)

                     If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
                     contains the norm of the off-diagonal part of the j-th column
                     of A.  If TRANS = 'N', CNORM(j) must be greater than or equal
                     to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
                     must be greater than or equal to the 1-norm.

                     If NORMIN = 'N', CNORM is an output argument and CNORM(j)
                     returns the 1-norm of the offdiagonal part of the j-th column
                     of A.

           INFO

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

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

       Further Details:

             A rough bound on x is computed; if that is less than overflow, DTRSV
             is called, otherwise, specific code is used which checks for possible
             overflow or divide-by-zero at every operation.

             A columnwise scheme is used for solving A*x = b.  The basic algorithm
             if A is lower triangular is

                  x[1:n] := b[1:n]
                  for j = 1, ..., n
                       x(j) := x(j) / A(j,j)
                       x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
                  end

             Define bounds on the components of x after j iterations of the loop:
                M(j) = bound on x[1:j]
                G(j) = bound on x[j+1:n]
             Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.

             Then for iteration j+1 we have
                M(j+1) <= G(j) / | A(j+1,j+1) |
                G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
                       <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )

             where CNORM(j+1) is greater than or equal to the infinity-norm of
             column j+1 of A, not counting the diagonal.  Hence

                G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
                             1<=i<=j
             and

                |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
                                              1<=i< j

             Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTRSV if the
             reciprocal of the largest M(j), j=1,..,n, is larger than
             max(underflow, 1/overflow).

             The bound on x(j) is also used to determine when a step in the
             columnwise method can be performed without fear of overflow.  If
             the computed bound is greater than a large constant, x is scaled to
             prevent overflow, but if the bound overflows, x is set to 0, x(j) to
             1, and scale to 0, and a non-trivial solution to A*x = 0 is found.

             Similarly, a row-wise scheme is used to solve A**T*x = b.  The basic
             algorithm for A upper triangular is

                  for j = 1, ..., n
                       x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j)
                  end

             We simultaneously compute two bounds
                  G(j) = bound on ( b(i) - A[1:i-1,i]**T * x[1:i-1] ), 1<=i<=j
                  M(j) = bound on x(i), 1<=i<=j

             The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
             add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
             Then the bound on x(j) is

                  M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |

                       <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
                                 1<=i<=j

             and we can safely call DTRSV if 1/M(n) and 1/G(n) are both greater
             than max(underflow, 1/overflow).

   subroutine dlauu2 (character UPLO, integer N, double precision, dimension( lda, * ) A, integer
       LDA, integer INFO)
       DLAUU2 computes the product UUH or LHL, where U and L are upper or lower triangular
       matrices (unblocked algorithm).

       Purpose:

            DLAUU2 computes the product U * U**T or L**T * L, where the triangular
            factor U or L is stored in the upper or lower triangular part of
            the array A.

            If UPLO = 'U' or 'u' then the upper triangle of the result is stored,
            overwriting the factor U in A.
            If UPLO = 'L' or 'l' then the lower triangle of the result is stored,
            overwriting the factor L in A.

            This is the unblocked form of the algorithm, calling Level 2 BLAS.

       Parameters:
           UPLO

                     UPLO is CHARACTER*1
                     Specifies whether the triangular factor stored in the array A
                     is upper or lower triangular:
                     = 'U':  Upper triangular
                     = 'L':  Lower triangular

           N

                     N is INTEGER
                     The order of the triangular factor U or L.  N >= 0.

           A

                     A is DOUBLE PRECISION array, dimension (LDA,N)
                     On entry, the triangular factor U or L.
                     On exit, if UPLO = 'U', the upper triangle of A is
                     overwritten with the upper triangle of the product U * U**T;
                     if UPLO = 'L', the lower triangle of A is overwritten with
                     the lower triangle of the product L**T * L.

           LDA

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

           INFO

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

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine dlauum (character UPLO, integer N, double precision, dimension( lda, * ) A, integer
       LDA, integer INFO)
       DLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular
       matrices (blocked algorithm).

       Purpose:

            DLAUUM computes the product U * U**T or L**T * L, where the triangular
            factor U or L is stored in the upper or lower triangular part of
            the array A.

            If UPLO = 'U' or 'u' then the upper triangle of the result is stored,
            overwriting the factor U in A.
            If UPLO = 'L' or 'l' then the lower triangle of the result is stored,
            overwriting the factor L in A.

            This is the blocked form of the algorithm, calling Level 3 BLAS.

       Parameters:
           UPLO

                     UPLO is CHARACTER*1
                     Specifies whether the triangular factor stored in the array A
                     is upper or lower triangular:
                     = 'U':  Upper triangular
                     = 'L':  Lower triangular

           N

                     N is INTEGER
                     The order of the triangular factor U or L.  N >= 0.

           A

                     A is DOUBLE PRECISION array, dimension (LDA,N)
                     On entry, the triangular factor U or L.
                     On exit, if UPLO = 'U', the upper triangle of A is
                     overwritten with the upper triangle of the product U * U**T;
                     if UPLO = 'L', the lower triangle of A is overwritten with
                     the lower triangle of the product L**T * L.

           LDA

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

           INFO

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

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

   subroutine drscl (integer N, double precision SA, double precision, dimension( * ) SX, integer
       INCX)
       DRSCL multiplies a vector by the reciprocal of a real scalar.

       Purpose:

            DRSCL multiplies an n-element real vector x by the real scalar 1/a.
            This is done without overflow or underflow as long as
            the final result x/a does not overflow or underflow.

       Parameters:
           N

                     N is INTEGER
                     The number of components of the vector x.

           SA

                     SA is DOUBLE PRECISION
                     The scalar a which is used to divide each component of x.
                     SA must be >= 0, or the subroutine will divide by zero.

           SX

                     SX is DOUBLE PRECISION array, dimension
                                    (1+(N-1)*abs(INCX))
                     The n-element vector x.

           INCX

                     INCX is INTEGER
                     The increment between successive values of the vector SX.
                     > 0:  SX(1) = X(1) and SX(1+(i-1)*INCX) = x(i),     1< i<= n

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           November 2017

   subroutine dtprfb (character SIDE, character TRANS, character DIRECT, character STOREV,
       integer M, integer N, integer K, integer L, double precision, dimension( ldv, * ) V,
       integer LDV, double precision, dimension( ldt, * ) T, integer LDT, double precision,
       dimension( lda, * ) A, integer LDA, double precision, dimension( ldb, * ) B, integer LDB,
       double precision, dimension( ldwork, * ) WORK, integer LDWORK)
       DTPRFB applies a real or complex 'triangular-pentagonal' blocked reflector to a real or
       complex matrix, which is composed of two blocks.

       Purpose:

            DTPRFB applies a real "triangular-pentagonal" block reflector H or its
            transpose H**T to a real matrix C, which is composed of two
            blocks A and B, either from the left or right.

       Parameters:
           SIDE

                     SIDE is CHARACTER*1
                     = 'L': apply H or H**T from the Left
                     = 'R': apply H or H**T from the Right

           TRANS

                     TRANS is CHARACTER*1
                     = 'N': apply H (No transpose)
                     = 'T': apply H**T (Transpose)

           DIRECT

                     DIRECT is CHARACTER*1
                     Indicates how H is formed from a product of elementary
                     reflectors
                     = 'F': H = H(1) H(2) . . . H(k) (Forward)
                     = 'B': H = H(k) . . . H(2) H(1) (Backward)

           STOREV

                     STOREV is CHARACTER*1
                     Indicates how the vectors which define the elementary
                     reflectors are stored:
                     = 'C': Columns
                     = 'R': Rows

           M

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

           N

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

           K

                     K is INTEGER
                     The order of the matrix T, i.e. the number of elementary
                     reflectors whose product defines the block reflector.
                     K >= 0.

           L

                     L is INTEGER
                     The order of the trapezoidal part of V.
                     K >= L >= 0.  See Further Details.

           V

                     V is DOUBLE PRECISION array, dimension
                                           (LDV,K) if STOREV = 'C'
                                           (LDV,M) if STOREV = 'R' and SIDE = 'L'
                                           (LDV,N) if STOREV = 'R' and SIDE = 'R'
                     The pentagonal matrix V, which contains the elementary reflectors
                     H(1), H(2), ..., H(K).  See Further Details.

           LDV

                     LDV is INTEGER
                     The leading dimension of the array V.
                     If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
                     if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
                     if STOREV = 'R', LDV >= K.

           T

                     T is DOUBLE PRECISION array, dimension (LDT,K)
                     The triangular K-by-K matrix T in the representation of the
                     block reflector.

           LDT

                     LDT is INTEGER
                     The leading dimension of the array T.
                     LDT >= K.

           A

                     A is DOUBLE PRECISION array, dimension
                     (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R'
                     On entry, the K-by-N or M-by-K matrix A.
                     On exit, A is overwritten by the corresponding block of
                     H*C or H**T*C or C*H or C*H**T.  See Further Details.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.
                     If SIDE = 'L', LDC >= max(1,K);
                     If SIDE = 'R', LDC >= max(1,M).

           B

                     B is DOUBLE PRECISION array, dimension (LDB,N)
                     On entry, the M-by-N matrix B.
                     On exit, B is overwritten by the corresponding block of
                     H*C or H**T*C or C*H or C*H**T.  See Further Details.

           LDB

                     LDB is INTEGER
                     The leading dimension of the array B.
                     LDB >= max(1,M).

           WORK

                     WORK is DOUBLE PRECISION array, dimension
                     (LDWORK,N) if SIDE = 'L',
                     (LDWORK,K) if SIDE = 'R'.

           LDWORK

                     LDWORK is INTEGER
                     The leading dimension of the array WORK.
                     If SIDE = 'L', LDWORK >= K;
                     if SIDE = 'R', LDWORK >= M.

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

       Further Details:

             The matrix C is a composite matrix formed from blocks A and B.
             The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K,
             and if SIDE = 'L', A is of size K-by-N.

             If SIDE = 'R' and DIRECT = 'F', C = [A B].

             If SIDE = 'L' and DIRECT = 'F', C = [A]
                                                 [B].

             If SIDE = 'R' and DIRECT = 'B', C = [B A].

             If SIDE = 'L' and DIRECT = 'B', C = [B]
                                                 [A].

             The pentagonal matrix V is composed of a rectangular block V1 and a
             trapezoidal block V2.  The size of the trapezoidal block is determined by
             the parameter L, where 0<=L<=K.  If L=K, the V2 block of V is triangular;
             if L=0, there is no trapezoidal block, thus V = V1 is rectangular.

             If DIRECT = 'F' and STOREV = 'C':  V = [V1]
                                                    [V2]
                - V2 is upper trapezoidal (first L rows of K-by-K upper triangular)

             If DIRECT = 'F' and STOREV = 'R':  V = [V1 V2]

                - V2 is lower trapezoidal (first L columns of K-by-K lower triangular)

             If DIRECT = 'B' and STOREV = 'C':  V = [V2]
                                                    [V1]
                - V2 is lower trapezoidal (last L rows of K-by-K lower triangular)

             If DIRECT = 'B' and STOREV = 'R':  V = [V2 V1]

                - V2 is upper trapezoidal (last L columns of K-by-K upper triangular)

             If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K.

             If STOREV = 'C' and SIDE = 'R', V is N-by-K with V2 L-by-K.

             If STOREV = 'R' and SIDE = 'L', V is K-by-M with V2 K-by-L.

             If STOREV = 'R' and SIDE = 'R', V is K-by-N with V2 K-by-L.

   subroutine slatrd (character UPLO, integer N, integer NB, real, dimension( lda, * ) A, integer
       LDA, real, dimension( * ) E, real, dimension( * ) TAU, real, dimension( ldw, * ) W,
       integer LDW)
       SLATRD reduces the first nb rows and columns of a symmetric/Hermitian matrix A to real
       tridiagonal form by an orthogonal similarity transformation.

       Purpose:

            SLATRD reduces NB rows and columns of a real symmetric matrix A to
            symmetric tridiagonal form by an orthogonal similarity
            transformation Q**T * A * Q, and returns the matrices V and W which are
            needed to apply the transformation to the unreduced part of A.

            If UPLO = 'U', SLATRD reduces the last NB rows and columns of a
            matrix, of which the upper triangle is supplied;
            if UPLO = 'L', SLATRD reduces the first NB rows and columns of a
            matrix, of which the lower triangle is supplied.

            This is an auxiliary routine called by SSYTRD.

       Parameters:
           UPLO

                     UPLO is CHARACTER*1
                     Specifies whether the upper or lower triangular part of the
                     symmetric matrix A is stored:
                     = 'U': Upper triangular
                     = 'L': Lower triangular

           N

                     N is INTEGER
                     The order of the matrix A.

           NB

                     NB is INTEGER
                     The number of rows and columns to be reduced.

           A

                     A is REAL array, dimension (LDA,N)
                     On entry, the symmetric matrix A.  If UPLO = 'U', the leading
                     n-by-n upper triangular part of A contains the upper
                     triangular part of the matrix A, and the strictly lower
                     triangular part of A is not referenced.  If UPLO = 'L', the
                     leading n-by-n lower triangular part of A contains the lower
                     triangular part of the matrix A, and the strictly upper
                     triangular part of A is not referenced.
                     On exit:
                     if UPLO = 'U', the last NB columns have been reduced to
                       tridiagonal form, with the diagonal elements overwriting
                       the diagonal elements of A; the elements above the diagonal
                       with the array TAU, represent the orthogonal matrix Q as a
                       product of elementary reflectors;
                     if UPLO = 'L', the first NB columns have been reduced to
                       tridiagonal form, with the diagonal elements overwriting
                       the diagonal elements of A; the elements below the diagonal
                       with the array TAU, represent the  orthogonal matrix Q as a
                       product of elementary reflectors.
                     See Further Details.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= (1,N).

           E

                     E is REAL array, dimension (N-1)
                     If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal
                     elements of the last NB columns of the reduced matrix;
                     if UPLO = 'L', E(1:nb) contains the subdiagonal elements of
                     the first NB columns of the reduced matrix.

           TAU

                     TAU is REAL array, dimension (N-1)
                     The scalar factors of the elementary reflectors, stored in
                     TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.
                     See Further Details.

           W

                     W is REAL array, dimension (LDW,NB)
                     The n-by-nb matrix W required to update the unreduced part
                     of A.

           LDW

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

       Author:
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date:
           December 2016

       Further Details:

             If UPLO = 'U', the matrix Q is represented as a product of elementary
             reflectors

                Q = H(n) H(n-1) . . . H(n-nb+1).

             Each H(i) has the form

                H(i) = I - tau * v * v**T

             where tau is a real scalar, and v is a real vector with
             v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i),
             and tau in TAU(i-1).

             If UPLO = 'L', the matrix Q is represented as a product of elementary
             reflectors

                Q = H(1) H(2) . . . H(nb).

             Each H(i) has the form

                H(i) = I - tau * v * v**T

             where tau is a real scalar, and v is a real vector with
             v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
             and tau in TAU(i).

             The elements of the vectors v together form the n-by-nb matrix V
             which is needed, with W, to apply the transformation to the unreduced
             part of the matrix, using a symmetric rank-2k update of the form:
             A := A - V*W**T - W*V**T.

             The contents of A on exit are illustrated by the following examples
             with n = 5 and nb = 2:

             if UPLO = 'U':                       if UPLO = 'L':

               (  a   a   a   v4  v5 )              (  d                  )
               (      a   a   v4  v5 )              (  1   d              )
               (          a   1   v5 )              (  v1  1   a          )
               (              d   1  )              (  v1  v2  a   a      )
               (                  d  )              (  v1  v2  a   a   a  )

             where d denotes a diagonal element of the reduced matrix, a denotes
             an element of the original matrix that is unchanged, and vi denotes
             an element of the vector defining H(i).

Author

       Generated automatically by Doxygen for LAPACK from the source code.