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

NAME

       complex16OTHERauxiliary - complex16

SYNOPSIS

   Functions
       subroutine clag2z (M, N, SA, LDSA, A, LDA, INFO)
           CLAG2Z converts a complex single precision matrix to a complex double precision
           matrix.
       double precision function dzsum1 (N, CX, INCX)
           DZSUM1 forms the 1-norm of the complex vector using the true absolute value.
       integer function ilazlc (M, N, A, LDA)
           ILAZLC scans a matrix for its last non-zero column.
       integer function ilazlr (M, N, A, LDA)
           ILAZLR scans a matrix for its last non-zero row.
       subroutine zdrscl (N, SA, SX, INCX)
           ZDRSCL multiplies a vector by the reciprocal of a real scalar.
       subroutine zlabrd (M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y, LDY)
           ZLABRD reduces the first nb rows and columns of a general matrix to a bidiagonal form.
       subroutine zlacgv (N, X, INCX)
           ZLACGV conjugates a complex vector.
       subroutine zlacn2 (N, V, X, EST, KASE, ISAVE)
           ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for
           evaluating matrix-vector products.
       subroutine zlacon (N, V, X, EST, KASE)
           ZLACON estimates the 1-norm of a square matrix, using reverse communication for
           evaluating matrix-vector products.
       subroutine zlacp2 (UPLO, M, N, A, LDA, B, LDB)
           ZLACP2 copies all or part of a real two-dimensional array to a complex array.
       subroutine zlacpy (UPLO, M, N, A, LDA, B, LDB)
           ZLACPY copies all or part of one two-dimensional array to another.
       subroutine zlacrm (M, N, A, LDA, B, LDB, C, LDC, RWORK)
           ZLACRM multiplies a complex matrix by a square real matrix.
       subroutine zlacrt (N, CX, INCX, CY, INCY, C, S)
           ZLACRT performs a linear transformation of a pair of complex vectors.
       complex *16 function zladiv (X, Y)
           ZLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
       subroutine zlaein (RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK, EPS3, SMLNUM, INFO)
           ZLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by
           inverse iteration.
       subroutine zlaev2 (A, B, C, RT1, RT2, CS1, SN1)
           ZLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian
           matrix.
       subroutine zlag2c (M, N, A, LDA, SA, LDSA, INFO)
           ZLAG2C converts a complex double precision matrix to a complex single precision
           matrix.
       subroutine zlags2 (UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV, SNV, CSQ, SNQ)
           ZLAGS2
       subroutine zlagtm (TRANS, N, NRHS, ALPHA, DL, D, DU, X, LDX, BETA, B, LDB)
           ZLAGTM 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 zlahqr (WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)
           ZLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,
           using the double-shift/single-shift QR algorithm.
       subroutine zlahr2 (N, K, NB, A, LDA, TAU, T, LDT, Y, LDY)
           ZLAHR2 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 zlaic1 (JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
           ZLAIC1 applies one step of incremental condition estimation.
       double precision function zlangt (NORM, N, DL, D, DU)
           ZLANGT 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 zlanhb (NORM, UPLO, N, K, AB, LDAB, WORK)
           ZLANHB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,
           or the element of largest absolute value of a Hermitian band matrix.
       double precision function zlanhp (NORM, UPLO, N, AP, WORK)
           ZLANHP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,
           or the element of largest absolute value of a complex Hermitian matrix supplied in
           packed form.
       double precision function zlanhs (NORM, N, A, LDA, WORK)
           ZLANHS 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 zlanht (NORM, N, D, E)
           ZLANHT returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,
           or the element of largest absolute value of a complex Hermitian tridiagonal matrix.
       double precision function zlansb (NORM, UPLO, N, K, AB, LDAB, WORK)
           ZLANSB 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 zlansp (NORM, UPLO, N, AP, WORK)
           ZLANSP 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 zlantb (NORM, UPLO, DIAG, N, K, AB, LDAB, WORK)
           ZLANTB 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 zlantp (NORM, UPLO, DIAG, N, AP, WORK)
           ZLANTP 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 zlantr (NORM, UPLO, DIAG, M, N, A, LDA, WORK)
           ZLANTR 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 zlapll (N, X, INCX, Y, INCY, SSMIN)
           ZLAPLL measures the linear dependence of two vectors.
       subroutine zlapmr (FORWRD, M, N, X, LDX, K)
           ZLAPMR rearranges rows of a matrix as specified by a permutation vector.
       subroutine zlapmt (FORWRD, M, N, X, LDX, K)
           ZLAPMT performs a forward or backward permutation of the columns of a matrix.
       subroutine zlaqhb (UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED)
           ZLAQHB scales a Hermitian band matrix, using scaling factors computed by cpbequ.
       subroutine zlaqhp (UPLO, N, AP, S, SCOND, AMAX, EQUED)
           ZLAQHP scales a Hermitian matrix stored in packed form.
       subroutine zlaqp2 (M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, WORK)
           ZLAQP2 computes a QR factorization with column pivoting of the matrix block.
       subroutine zlaqps (M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, VN2, AUXV, F, LDF)
           ZLAQPS computes a step of QR factorization with column pivoting of a real m-by-n
           matrix A by using BLAS level 3.
       subroutine zlaqr0 (WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK,
           INFO)
           ZLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices
           from the Schur decomposition.
       subroutine zlaqr1 (N, H, LDH, S1, S2, V)
           ZLAQR1 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 zlaqr2 (WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND,
           SH, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)
           ZLAQR2 performs the unitary similarity transformation of a Hessenberg matrix to detect
           and deflate fully converged eigenvalues from a trailing principal submatrix
           (aggressive early deflation).
       subroutine zlaqr3 (WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND,
           SH, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)
           ZLAQR3 performs the unitary similarity transformation of a Hessenberg matrix to detect
           and deflate fully converged eigenvalues from a trailing principal submatrix
           (aggressive early deflation).
       subroutine zlaqr4 (WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK,
           INFO)
           ZLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices
           from the Schur decomposition.
       subroutine zlaqr5 (WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S, H, LDH, ILOZ, IHIZ, Z,
           LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH)
           ZLAQR5 performs a single small-bulge multi-shift QR sweep.
       subroutine zlaqsb (UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED)
           ZLAQSB scales a symmetric/Hermitian band matrix, using scaling factors computed by
           spbequ.
       subroutine zlaqsp (UPLO, N, AP, S, SCOND, AMAX, EQUED)
           ZLAQSP scales a symmetric/Hermitian matrix in packed storage, using scaling factors
           computed by sppequ.
       subroutine zlar1v (N, B1, BN, LAMBDA, D, L, LD, LLD, PIVMIN, GAPTOL, Z, WANTNC, NEGCNT,
           ZTZ, MINGMA, R, ISUPPZ, NRMINV, RESID, RQCORR, WORK)
           ZLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1
           through bn of the tridiagonal matrix LDLT - λI.
       subroutine zlar2v (N, X, Y, Z, INCX, C, S, INCC)
           ZLAR2V applies a vector of plane rotations with real cosines and complex sines from
           both sides to a sequence of 2-by-2 symmetric/Hermitian matrices.
       subroutine zlarcm (M, N, A, LDA, B, LDB, C, LDC, RWORK)
           ZLARCM copies all or part of a real two-dimensional array to a complex array.
       subroutine zlarf (SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
           ZLARF applies an elementary reflector to a general rectangular matrix.
       subroutine zlarfb (SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK,
           LDWORK)
           ZLARFB applies a block reflector or its conjugate-transpose to a general rectangular
           matrix.
       subroutine zlarfb_gett (IDENT, M, N, K, T, LDT, A, LDA, B, LDB, WORK, LDWORK)
           ZLARFB_GETT
       subroutine zlarfg (N, ALPHA, X, INCX, TAU)
           ZLARFG generates an elementary reflector (Householder matrix).
       subroutine zlarfgp (N, ALPHA, X, INCX, TAU)
           ZLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
       subroutine zlarft (DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
           ZLARFT forms the triangular factor T of a block reflector H = I - vtvH
       subroutine zlarfx (SIDE, M, N, V, TAU, C, LDC, WORK)
           ZLARFX applies an elementary reflector to a general rectangular matrix, with loop
           unrolling when the reflector has order ≤ 10.
       subroutine zlarfy (UPLO, N, V, INCV, TAU, C, LDC, WORK)
           ZLARFY
       subroutine zlargv (N, X, INCX, Y, INCY, C, INCC)
           ZLARGV generates a vector of plane rotations with real cosines and complex sines.
       subroutine zlarnv (IDIST, ISEED, N, X)
           ZLARNV returns a vector of random numbers from a uniform or normal distribution.
       subroutine zlarrv (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)
           ZLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and
           the eigenvalues of L D LT.
       subroutine zlartv (N, X, INCX, Y, INCY, C, S, INCC)
           ZLARTV applies a vector of plane rotations with real cosines and complex sines to the
           elements of a pair of vectors.
       subroutine zlascl (TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
           ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
       subroutine zlaset (UPLO, M, N, ALPHA, BETA, A, LDA)
           ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to
           given values.
       subroutine zlasr (SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
           ZLASR applies a sequence of plane rotations to a general rectangular matrix.
       subroutine zlaswp (N, A, LDA, K1, K2, IPIV, INCX)
           ZLASWP performs a series of row interchanges on a general rectangular matrix.
       subroutine zlat2c (UPLO, N, A, LDA, SA, LDSA, INFO)
           ZLAT2C converts a double complex triangular matrix to a complex triangular matrix.
       subroutine zlatbs (UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE, CNORM, INFO)
           ZLATBS solves a triangular banded system of equations.
       subroutine zlatdf (IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, JPIV)
           ZLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes
           a contribution to the reciprocal Dif-estimate.
       subroutine zlatps (UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE, CNORM, INFO)
           ZLATPS solves a triangular system of equations with the matrix held in packed storage.
       subroutine zlatrd (UPLO, N, NB, A, LDA, E, TAU, W, LDW)
           ZLATRD reduces the first nb rows and columns of a symmetric/Hermitian matrix A to real
           tridiagonal form by an unitary similarity transformation.
       subroutine zlatrs (UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
           ZLATRS solves a triangular system of equations with the scale factor set to prevent
           overflow.
       subroutine zlauu2 (UPLO, N, A, LDA, INFO)
           ZLAUU2 computes the product UUH or LHL, where U and L are upper or lower triangular
           matrices (unblocked algorithm).
       subroutine zlauum (UPLO, N, A, LDA, INFO)
           ZLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular
           matrices (blocked algorithm).
       subroutine zrot (N, CX, INCX, CY, INCY, C, S)
           ZROT applies a plane rotation with real cosine and complex sine to a pair of complex
           vectors.
       subroutine zspmv (UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
           ZSPMV computes a matrix-vector product for complex vectors using a complex symmetric
           packed matrix
       subroutine zspr (UPLO, N, ALPHA, X, INCX, AP)
           ZSPR performs the symmetrical rank-1 update of a complex symmetric packed matrix.
       subroutine ztprfb (SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V, LDV, T, LDT, A, LDA, B,
           LDB, WORK, LDWORK)
           ZTPRFB applies a complex 'triangular-pentagonal' block reflector to a complex matrix,
           which is composed of two blocks.

Detailed Description

       This is the group of complex16 other auxiliary routines

Function Documentation

   subroutine clag2z (integer M, integer N, complex, dimension( ldsa, * ) SA, integer LDSA,
       complex*16, dimension( lda, * ) A, integer LDA, integer INFO)
       CLAG2Z converts a complex single precision matrix to a complex double precision matrix.

       Purpose:

            CLAG2Z converts a COMPLEX matrix, SA, to a COMPLEX*16 matrix, A.

            Note that while it is possible to overflow while converting
            from double to single, it is not possible to overflow when
            converting from single to double.

            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.

           SA

                     SA is COMPLEX array, dimension (LDSA,N)
                     On entry, the M-by-N coefficient matrix SA.

           LDSA

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

           A

                     A is COMPLEX*16 array, dimension (LDA,N)
                     On exit, the M-by-N coefficient matrix A.

           LDA

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

           INFO

                     INFO is INTEGER
                     = 0:  successful exit

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   double precision function dzsum1 (integer N, complex*16, dimension( * ) CX, integer INCX)
       DZSUM1 forms the 1-norm of the complex vector using the true absolute value.

       Purpose:

            DZSUM1 takes the sum of the absolute values of a complex
            vector and returns a double precision result.

            Based on DZASUM from the Level 1 BLAS.
            The change is to use the 'genuine' absolute value.

       Parameters
           N

                     N is INTEGER
                     The number of elements in the vector CX.

           CX

                     CX is COMPLEX*16 array, dimension (N)
                     The vector whose elements will be summed.

           INCX

                     INCX is INTEGER
                     The spacing between successive values of CX.  INCX > 0.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Contributors:
           Nick Higham for use with ZLACON.

   integer function ilazlc (integer M, integer N, complex*16, dimension( lda, * ) A, integer LDA)
       ILAZLC scans a matrix for its last non-zero column.

       Purpose:

            ILAZLC scans A for its last non-zero column.

       Parameters
           M

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

           N

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

           A

                     A is COMPLEX*16 array, dimension (LDA,N)
                     The m by n matrix A.

           LDA

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   integer function ilazlr (integer M, integer N, complex*16, dimension( lda, * ) A, integer LDA)
       ILAZLR scans a matrix for its last non-zero row.

       Purpose:

            ILAZLR scans A for its last non-zero row.

       Parameters
           M

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

           N

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

           A

                     A is COMPLEX*16 array, dimension (LDA,N)
                     The m by n matrix A.

           LDA

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine zdrscl (integer N, double precision SA, complex*16, dimension( * ) SX, integer
       INCX)
       ZDRSCL multiplies a vector by the reciprocal of a real scalar.

       Purpose:

            ZDRSCL multiplies an n-element complex 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 COMPLEX*16 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.

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

       Purpose:

            ZLABRD reduces the first NB rows and columns of a complex general
            m by n matrix A to upper or lower real bidiagonal form by a unitary
            transformation Q**H * 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 ZGEBRD

       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 COMPLEX*16 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 unitary
                       matrix Q as a product of elementary reflectors; and
                       elements above the diagonal in the first NB rows, with the
                       array TAUP, represent the unitary 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 unitary
                       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 unitary 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 COMPLEX*16 array, dimension (NB)
                     The scalar factors of the elementary reflectors which
                     represent the unitary matrix Q. See Further Details.

           TAUP

                     TAUP is COMPLEX*16 array, dimension (NB)
                     The scalar factors of the elementary reflectors which
                     represent the unitary matrix P. See Further Details.

           X

                     X is COMPLEX*16 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 COMPLEX*16 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.

       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**H  and G(i) = I - taup * u * u**H

             where tauq and taup are complex scalars, and v and u are complex
             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**H 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**H - X*U**H.

             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 zlacgv (integer N, complex*16, dimension( * ) X, integer INCX)
       ZLACGV conjugates a complex vector.

       Purpose:

            ZLACGV conjugates a complex vector of length N.

       Parameters
           N

                     N is INTEGER
                     The length of the vector X.  N >= 0.

           X

                     X is COMPLEX*16 array, dimension
                                    (1+(N-1)*abs(INCX))
                     On entry, the vector of length N to be conjugated.
                     On exit, X is overwritten with conjg(X).

           INCX

                     INCX is INTEGER
                     The spacing between successive elements of X.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

       Purpose:

            ZLACN2 estimates the 1-norm of a square, complex 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 COMPLEX*16 array, dimension (N)
                    On the final return, V = A*W,  where  EST = norm(V)/norm(W)
                    (W is not returned).

           X

                     X is COMPLEX*16 array, dimension (N)
                    On an intermediate return, X should be overwritten by
                          A * X,   if KASE=1,
                          A**H * X,  if KASE=2,
                    where A**H is the conjugate transpose of A, and ZLACN2 must be
                    re-called with all the other parameters unchanged.

           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 ZLACN2.
                    On exit, EST is an estimate (a lower bound) for norm(A).

           KASE

                     KASE is INTEGER
                    On the initial call to ZLACN2, 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**H * X.
                    On the final return from ZLACN2, KASE will again be 0.

           ISAVE

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Further Details:

             Originally named CONEST, dated March 16, 1988.

             Last modified:  April, 1999

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

                ZLACON     ZLACN2
                 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 zlacon (integer N, complex*16, dimension( n ) V, complex*16, dimension( n ) X,
       double precision EST, integer KASE)
       ZLACON estimates the 1-norm of a square matrix, using reverse communication for evaluating
       matrix-vector products.

       Purpose:

            ZLACON estimates the 1-norm of a square, complex 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 COMPLEX*16 array, dimension (N)
                    On the final return, V = A*W,  where  EST = norm(V)/norm(W)
                    (W is not returned).

           X

                     X is COMPLEX*16 array, dimension (N)
                    On an intermediate return, X should be overwritten by
                          A * X,   if KASE=1,
                          A**H * X,  if KASE=2,
                    where A**H is the conjugate transpose of A, and ZLACON must be
                    re-called with all the other parameters unchanged.

           EST

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

           KASE

                     KASE is INTEGER
                    On the initial call to ZLACON, 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**H * X.
                    On the final return from ZLACON, KASE will again be 0.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Further Details:
           Originally named CONEST, dated March 16, 1988.
            Last modified: April, 1999

       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 zlacp2 (character UPLO, integer M, integer N, double precision, dimension( lda, * )
       A, integer LDA, complex*16, dimension( ldb, * ) B, integer LDB)
       ZLACP2 copies all or part of a real two-dimensional array to a complex array.

       Purpose:

            ZLACP2 copies all or part of a real two-dimensional matrix A to a
            complex matrix B.

       Parameters
           UPLO

                     UPLO is CHARACTER*1
                     Specifies the part of the matrix A to be copied to B.
                     = 'U':      Upper triangular part
                     = 'L':      Lower triangular part
                     Otherwise:  All of the matrix A

           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.

           A

                     A is DOUBLE PRECISION array, dimension (LDA,N)
                     The m by n matrix A.  If UPLO = 'U', only the upper trapezium
                     is accessed; if UPLO = 'L', only the lower trapezium is
                     accessed.

           LDA

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

           B

                     B is COMPLEX*16 array, dimension (LDB,N)
                     On exit, B = A in the locations specified by UPLO.

           LDB

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine zlacpy (character UPLO, integer M, integer N, complex*16, dimension( lda, * ) A,
       integer LDA, complex*16, dimension( ldb, * ) B, integer LDB)
       ZLACPY copies all or part of one two-dimensional array to another.

       Purpose:

            ZLACPY copies all or part of a two-dimensional matrix A to another
            matrix B.

       Parameters
           UPLO

                     UPLO is CHARACTER*1
                     Specifies the part of the matrix A to be copied to B.
                     = 'U':      Upper triangular part
                     = 'L':      Lower triangular part
                     Otherwise:  All of the matrix A

           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.

           A

                     A is COMPLEX*16 array, dimension (LDA,N)
                     The m by n matrix A.  If UPLO = 'U', only the upper trapezium
                     is accessed; if UPLO = 'L', only the lower trapezium is
                     accessed.

           LDA

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

           B

                     B is COMPLEX*16 array, dimension (LDB,N)
                     On exit, B = A in the locations specified by UPLO.

           LDB

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine zlacrm (integer M, integer N, complex*16, dimension( lda, * ) A, integer LDA,
       double precision, dimension( ldb, * ) B, integer LDB, complex*16, dimension( ldc, * ) C,
       integer LDC, double precision, dimension( * ) RWORK)
       ZLACRM multiplies a complex matrix by a square real matrix.

       Purpose:

            ZLACRM performs a very simple matrix-matrix multiplication:
                     C := A * B,
            where A is M by N and complex; B is N by N and real;
            C is M by N and complex.

       Parameters
           M

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

           N

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

           A

                     A is COMPLEX*16 array, dimension (LDA, N)
                     On entry, A contains the M by N matrix A.

           LDA

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

           B

                     B is DOUBLE PRECISION array, dimension (LDB, N)
                     On entry, B contains the N by N matrix B.

           LDB

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

           C

                     C is COMPLEX*16 array, dimension (LDC, N)
                     On exit, C contains the M by N matrix C.

           LDC

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

           RWORK

                     RWORK is DOUBLE PRECISION array, dimension (2*M*N)

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine zlacrt (integer N, complex*16, dimension( * ) CX, integer INCX, complex*16,
       dimension( * ) CY, integer INCY, complex*16 C, complex*16 S)
       ZLACRT performs a linear transformation of a pair of complex vectors.

       Purpose:

            ZLACRT performs the operation

               (  c  s )( x )  ==> ( x )
               ( -s  c )( y )      ( y )

            where c and s are complex and the vectors x and y are complex.

       Parameters
           N

                     N is INTEGER
                     The number of elements in the vectors CX and CY.

           CX

                     CX is COMPLEX*16 array, dimension (N)
                     On input, the vector x.
                     On output, CX is overwritten with c*x + s*y.

           INCX

                     INCX is INTEGER
                     The increment between successive values of CX.  INCX <> 0.

           CY

                     CY is COMPLEX*16 array, dimension (N)
                     On input, the vector y.
                     On output, CY is overwritten with -s*x + c*y.

           INCY

                     INCY is INTEGER
                     The increment between successive values of CY.  INCY <> 0.

           C

                     C is COMPLEX*16

           S

                     S is COMPLEX*16
                     C and S define the matrix
                        [  C   S  ].
                        [ -S   C  ]

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   complex*16 function zladiv (complex*16 X, complex*16 Y)
       ZLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.

       Purpose:

            ZLADIV := X / Y, where X and Y are complex.  The computation of X / Y
            will not overflow on an intermediary step unless the results
            overflows.

       Parameters
           X

                     X is COMPLEX*16

           Y

                     Y is COMPLEX*16
                     The complex scalars X and Y.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine zlaein (logical RIGHTV, logical NOINIT, integer N, complex*16, dimension( ldh, * )
       H, integer LDH, complex*16 W, complex*16, dimension( * ) V, complex*16, dimension( ldb, *
       ) B, integer LDB, double precision, dimension( * ) RWORK, double precision EPS3, double
       precision SMLNUM, integer INFO)
       ZLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by
       inverse iteration.

       Purpose:

            ZLAEIN uses inverse iteration to find a right or left eigenvector
            corresponding to the eigenvalue W of a complex 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 V
                     = .FALSE.: initial vector supplied in V.

           N

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

           H

                     H is COMPLEX*16 array, dimension (LDH,N)
                     The upper Hessenberg matrix H.

           LDH

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

           W

                     W is COMPLEX*16
                     The eigenvalue of H whose corresponding right or left
                     eigenvector is to be computed.

           V

                     V is COMPLEX*16 array, dimension (N)
                     On entry, if NOINIT = .FALSE., V must contain a starting
                     vector for inverse iteration; otherwise V need not be set.
                     On exit, V contains the computed eigenvector, 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|.

           B

                     B is COMPLEX*16 array, dimension (LDB,N)

           LDB

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

           RWORK

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

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     = 1:  inverse iteration did not converge; V is set to the
                           last iterate.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine zlaev2 (complex*16 A, complex*16 B, complex*16 C, double precision RT1, double
       precision RT2, double precision CS1, complex*16 SN1)
       ZLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.

       Purpose:

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

            [ CS1  CONJG(SN1) ] [    A     B ] [ CS1 -CONJG(SN1) ] = [ RT1  0  ]
            [-SN1     CS1     ] [ CONJG(B) C ] [ SN1     CS1     ]   [  0  RT2 ].

       Parameters
           A

                     A is COMPLEX*16
                    The (1,1) element of the 2-by-2 matrix.

           B

                     B is COMPLEX*16
                    The (1,2) element and the conjugate of the (2,1) element of
                    the 2-by-2 matrix.

           C

                     C is COMPLEX*16
                    The (2,2) element of the 2-by-2 matrix.

           RT1

                     RT1 is DOUBLE PRECISION
                    The eigenvalue of larger absolute value.

           RT2

                     RT2 is DOUBLE PRECISION
                    The eigenvalue of smaller absolute value.

           CS1

                     CS1 is DOUBLE PRECISION

           SN1

                     SN1 is COMPLEX*16
                    The vector (CS1, SN1) is a unit right eigenvector for RT1.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Further Details:

             RT1 is accurate to a few ulps barring over/underflow.

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

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

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

   subroutine zlag2c (integer M, integer N, complex*16, dimension( lda, * ) A, integer LDA,
       complex, dimension( ldsa, * ) SA, integer LDSA, integer INFO)
       ZLAG2C converts a complex double precision matrix to a complex single precision matrix.

       Purpose:

            ZLAG2C converts a COMPLEX*16 matrix, SA, to a COMPLEX matrix, A.

            RMAX is the overflow for the SINGLE PRECISION arithmetic
            ZLAG2C 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 COMPLEX*16 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 COMPLEX 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.

   subroutine zlags2 (logical UPPER, double precision A1, complex*16 A2, double precision A3,
       double precision B1, complex*16 B2, double precision B3, double precision CSU, complex*16
       SNU, double precision CSV, complex*16 SNV, double precision CSQ, complex*16 SNQ)
       ZLAGS2

       Purpose:

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

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

            or if ( .NOT.UPPER ) then

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

              U = (   CSU    SNU ), V = (  CSV    SNV ),
                  ( -SNU**H  CSU )      ( -SNV**H CSV )

              Q = (   CSQ    SNQ )
                  ( -SNQ**H  CSQ )

            The rows of the transformed A and B are parallel. Moreover, if the
            input 2-by-2 matrix A is not zero, then the transformed (1,1) entry
            of A is not zero. If the input matrices A and B are both not zero,
            then the transformed (2,2) element of B is not zero, except when the
            first rows of input A and B are parallel and the second rows are
            zero.

       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 COMPLEX*16

           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 COMPLEX*16

           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 COMPLEX*16
                     The desired unitary matrix U.

           CSV

                     CSV is DOUBLE PRECISION

           SNV

                     SNV is COMPLEX*16
                     The desired unitary matrix V.

           CSQ

                     CSQ is DOUBLE PRECISION

           SNQ

                     SNQ is COMPLEX*16
                     The desired unitary matrix Q.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine zlagtm (character TRANS, integer N, integer NRHS, double precision ALPHA,
       complex*16, dimension( * ) DL, complex*16, dimension( * ) D, complex*16, dimension( * )
       DU, complex*16, dimension( ldx, * ) X, integer LDX, double precision BETA, complex*16,
       dimension( ldb, * ) B, integer LDB)
       ZLAGTM 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:

            ZLAGTM 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**T * X + beta * B
                     = 'C':  Conjugate transpose, B := alpha * A**H * X + beta * B

           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 COMPLEX*16 array, dimension (N-1)
                     The (n-1) sub-diagonal elements of T.

           D

                     D is COMPLEX*16 array, dimension (N)
                     The diagonal elements of T.

           DU

                     DU is COMPLEX*16 array, dimension (N-1)
                     The (n-1) super-diagonal elements of T.

           X

                     X is COMPLEX*16 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 COMPLEX*16 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.

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

       Purpose:

               ZLAHQR is an auxiliary routine called by CHSEQR to update the
               eigenvalues and Schur decomposition already computed by CHSEQR, 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 triangular in rows and
                     columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1).
                     ZLAHQR 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 COMPLEX*16 array, dimension (LDH,N)
                     On entry, the upper Hessenberg matrix H.
                     On exit, if INFO is zero and if WANTT is .TRUE., then H
                     is upper triangular in rows and columns ILO:IHI.  If INFO
                     is zero and if WANTT is .FALSE., then the contents of H
                     are unspecified on exit.  The output state of H in case
                     INF is positive is below under the description of INFO.

           LDH

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

           W

                     W is COMPLEX*16 array, dimension (N)
                     The computed eigenvalues ILO to IHI are stored in the
                     corresponding elements of W. If WANTT is .TRUE., the
                     eigenvalues are stored in the same order as on the diagonal
                     of the Schur form returned in H, with W(i) = H(i,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 COMPLEX*16 array, dimension (LDZ,N)
                     If WANTZ is .TRUE., on entry Z must contain the current
                     matrix Z of transformations accumulated by CHSEQR, 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
                      > 0:   if INFO = i, ZLAHQR failed to compute all the
                             eigenvalues ILO to IHI in a total of 30 iterations
                             per eigenvalue; elements i+1:ihi of W contain
                             those eigenvalues which have been successfully
                             computed.

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

                             If INFO > 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 triangular in
                             rows and columns INFO+1 through IHI.

                             If INFO > 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.

       Contributors:

                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 ZLAHQR 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 zlahr2 (integer N, integer K, integer NB, complex*16, dimension( lda, * ) A,
       integer LDA, complex*16, dimension( nb ) TAU, complex*16, dimension( ldt, nb ) T, integer
       LDT, complex*16, dimension( ldy, nb ) Y, integer LDY)
       ZLAHR2 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:

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

            This is an auxiliary routine called by ZGEHRD.

       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 COMPLEX*16 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 COMPLEX*16 array, dimension (NB)
                     The scalar factors of the elementary reflectors. See Further
                     Details.

           T

                     T is COMPLEX*16 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 COMPLEX*16 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.

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

             where tau is a complex scalar, and v is a complex 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**H) * (A - Y*V**H).

             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 ZLAHRD
             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 ZLAHRD routine. (This
             subroutine is not backward compatible with LAPACK-3.0's ZLAHRD.)

       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 zlaic1 (integer JOB, integer J, complex*16, dimension( j ) X, double precision
       SEST, complex*16, dimension( j ) W, complex*16 GAMMA, double precision SESTPR, complex*16
       S, complex*16 C)
       ZLAIC1 applies one step of incremental condition estimation.

       Purpose:

            ZLAIC1 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 ZLAIC1 computes sestpr, s, c such that
            the vector
                            [ s*x ]
                     xhat = [  c  ]
            is an approximate singular vector of
                            [ L       0  ]
                     Lhat = [ w**H 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]**H and sestpr**2 is an eigenpair of the system

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

            where  alpha =  x**H * 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 COMPLEX*16 array, dimension (J)
                     The j-vector x.

           SEST

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

           W

                     W is COMPLEX*16 array, dimension (J)
                     The j-vector w.

           GAMMA

                     GAMMA is COMPLEX*16
                     The diagonal element gamma.

           SESTPR

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

           S

                     S is COMPLEX*16
                     Sine needed in forming xhat.

           C

                     C is COMPLEX*16
                     Cosine needed in forming xhat.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   double precision function zlangt (character NORM, integer N, complex*16, dimension( * ) DL,
       complex*16, dimension( * ) D, complex*16, dimension( * ) DU)
       ZLANGT 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:

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

       Returns
           ZLANGT

               ZLANGT = ( 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 ZLANGT as described
                     above.

           N

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

           DL

                     DL is COMPLEX*16 array, dimension (N-1)
                     The (n-1) sub-diagonal elements of A.

           D

                     D is COMPLEX*16 array, dimension (N)
                     The diagonal elements of A.

           DU

                     DU is COMPLEX*16 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.

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

       Purpose:

            ZLANHB  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 hermitian band matrix A,  with k super-diagonals.

       Returns
           ZLANHB

               ZLANHB = ( 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 ZLANHB 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
                     = 'L':  Lower triangular

           N

                     N is INTEGER
                     The order of the matrix A.  N >= 0.  When N = 0, ZLANHB 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 COMPLEX*16 array, dimension (LDAB,N)
                     The upper or lower triangle of the hermitian 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 the imaginary parts of the diagonal elements need
                     not be set and are assumed to be zero.

           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.

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

       Purpose:

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

       Returns
           ZLANHP

               ZLANHP = ( 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 ZLANHP as described
                     above.

           UPLO

                     UPLO is CHARACTER*1
                     Specifies whether the upper or lower triangular part of the
                     hermitian 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, ZLANHP is
                     set to zero.

           AP

                     AP is COMPLEX*16 array, dimension (N*(N+1)/2)
                     The upper or lower triangle of the hermitian 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 the  imaginary parts of the diagonal elements need
                     not be set and are assumed to be zero.

           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.

   double precision function zlanhs (character NORM, integer N, complex*16, dimension( lda, * )
       A, integer LDA, double precision, dimension( * ) WORK)
       ZLANHS 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:

            ZLANHS  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
           ZLANHS

               ZLANHS = ( 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 ZLANHS as described
                     above.

           N

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

           A

                     A is COMPLEX*16 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.

   double precision function zlanht (character NORM, integer N, double precision, dimension( * )
       D, complex*16, dimension( * ) E)
       ZLANHT returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or
       the element of largest absolute value of a complex Hermitian tridiagonal matrix.

       Purpose:

            ZLANHT  returns the value of the one norm,  or the Frobenius norm, or
            the  infinity norm,  or the  element of  largest absolute value  of a
            complex Hermitian tridiagonal matrix A.

       Returns
           ZLANHT

               ZLANHT = ( 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 ZLANHT as described
                     above.

           N

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

           D

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

           E

                     E is COMPLEX*16 array, dimension (N-1)
                     The (n-1) sub-diagonal or super-diagonal elements of A.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   double precision function zlansb (character NORM, character UPLO, integer N, integer K,
       complex*16, dimension( ldab, * ) AB, integer LDAB, double precision, dimension( * ) WORK)
       ZLANSB 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:

            ZLANSB  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
           ZLANSB

               ZLANSB = ( 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 ZLANSB 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, ZLANSB 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 COMPLEX*16 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.

   double precision function zlansp (character NORM, character UPLO, integer N, complex*16,
       dimension( * ) AP, double precision, dimension( * ) WORK)
       ZLANSP 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:

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

       Returns
           ZLANSP

               ZLANSP = ( 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 ZLANSP 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, ZLANSP is
                     set to zero.

           AP

                     AP is COMPLEX*16 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.

   double precision function zlantb (character NORM, character UPLO, character DIAG, integer N,
       integer K, complex*16, dimension( ldab, * ) AB, integer LDAB, double precision, dimension(
       * ) WORK)
       ZLANTB 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:

            ZLANTB  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
           ZLANTB

               ZLANTB = ( 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 ZLANTB 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, ZLANTB 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 COMPLEX*16 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.

   double precision function zlantp (character NORM, character UPLO, character DIAG, integer N,
       complex*16, dimension( * ) AP, double precision, dimension( * ) WORK)
       ZLANTP 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:

            ZLANTP  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
           ZLANTP

               ZLANTP = ( 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 ZLANTP 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, ZLANTP is
                     set to zero.

           AP

                     AP is COMPLEX*16 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.

   double precision function zlantr (character NORM, character UPLO, character DIAG, integer M,
       integer N, complex*16, dimension( lda, * ) A, integer LDA, double precision, dimension( *
       ) WORK)
       ZLANTR 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:

            ZLANTR  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
           ZLANTR

               ZLANTR = ( 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 ZLANTR 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, ZLANTR 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, ZLANTR is set to zero.

           A

                     A is COMPLEX*16 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.

   subroutine zlapll (integer N, complex*16, dimension( * ) X, integer INCX, complex*16,
       dimension( * ) Y, integer INCY, double precision SSMIN)
       ZLAPLL 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 COMPLEX*16 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 COMPLEX*16 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.

   subroutine zlapmr (logical FORWRD, integer M, integer N, complex*16, dimension( ldx, * ) X,
       integer LDX, integer, dimension( * ) K)
       ZLAPMR rearranges rows of a matrix as specified by a permutation vector.

       Purpose:

            ZLAPMR 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 COMPLEX*16 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.

   subroutine zlapmt (logical FORWRD, integer M, integer N, complex*16, dimension( ldx, * ) X,
       integer LDX, integer, dimension( * ) K)
       ZLAPMT performs a forward or backward permutation of the columns of a matrix.

       Purpose:

            ZLAPMT 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 COMPLEX*16 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.

   subroutine zlaqhb (character UPLO, integer N, integer KD, complex*16, dimension( ldab, * ) AB,
       integer LDAB, double precision, dimension( * ) S, double precision SCOND, double precision
       AMAX, character EQUED)
       ZLAQHB scales a Hermitian band matrix, using scaling factors computed by cpbequ.

       Purpose:

            ZLAQHB equilibrates a Hermitian 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 COMPLEX*16 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**H *U or A = L*L**H 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.

   subroutine zlaqhp (character UPLO, integer N, complex*16, dimension( * ) AP, double precision,
       dimension( * ) S, double precision SCOND, double precision AMAX, character EQUED)
       ZLAQHP scales a Hermitian matrix stored in packed form.

       Purpose:

            ZLAQHP equilibrates a Hermitian 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
                     Hermitian 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 COMPLEX*16 array, dimension (N*(N+1)/2)
                     On entry, the upper or lower triangle of the Hermitian 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.

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

       Purpose:

            ZLAQP2 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 array, dimension (N)

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       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 zlaqps (integer M, integer N, integer OFFSET, integer NB, integer KB, complex*16,
       dimension( lda, * ) A, integer LDA, integer, dimension( * ) JPVT, complex*16, dimension( *
       ) TAU, double precision, dimension( * ) VN1, double precision, dimension( * ) VN2,
       complex*16, dimension( * ) AUXV, complex*16, dimension( ldf, * ) F, integer LDF)
       ZLAQPS computes a step of QR factorization with column pivoting of a real m-by-n matrix A
       by using BLAS level 3.

       Purpose:

            ZLAQPS computes a step of QR factorization with column pivoting
            of a complex 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 array, dimension (NB)
                     Auxiliary vector.

           F

                     F is COMPLEX*16 array, dimension (LDF,NB)
                     Matrix F**H = L * Y**H * 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.

       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 zlaqr0 (logical WANTT, logical WANTZ, integer N, integer ILO, integer IHI,
       complex*16, dimension( ldh, * ) H, integer LDH, complex*16, dimension( * ) W, integer
       ILOZ, integer IHIZ, complex*16, dimension( ldz, * ) Z, integer LDZ, complex*16, dimension(
       * ) WORK, integer LWORK, integer INFO)
       ZLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from
       the Schur decomposition.

       Purpose:

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

               Optionally Z may be postmultiplied into an input unitary
               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 unitary matrix Q:  A = Q*H*Q**H = (QZ)*H*(QZ)**H.

       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 triangular in rows
                      and columns 1:ILO-1 and IHI+1:N and, if ILO > 1,
                      H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
                      previous call to ZGEBAL, and then passed to ZGEHRD when the
                      matrix output by ZGEBAL is reduced to Hessenberg form.
                      Otherwise, ILO and IHI should be set to 1 and N,
                      respectively.  If N > 0, then 1 <= ILO <= IHI <= N.
                      If N = 0, then ILO = 1 and IHI = 0.

           H

                     H is COMPLEX*16 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 triangular matrix T from the Schur
                      decomposition (the Schur form). If INFO = 0 and WANT is
                      .FALSE., then the contents of H are unspecified on exit.
                      (The output value of H when INFO > 0 is given under the
                      description of INFO below.)

                      This subroutine may explicitly set H(i,j) = 0 for i > 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 >= max(1,N).

           W

                     W is COMPLEX*16 array, dimension (N)
                      The computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored
                      in W(ILO:IHI). 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 W(i) = H(i,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 COMPLEX*16 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 > 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 >= MAX(1,IHIZ).  Otherwise, LDZ >= 1.

           WORK

                     WORK is COMPLEX*16 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 >= 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 ZLAQR0 does a workspace query.
                      In this case, ZLAQR0 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
                        > 0:  if INFO = i, ZLAQR0 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 > 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 > 0 and WANTT is .TRUE., then on exit

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

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

                           If INFO > 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 unitary matrix in (*) (regard-
                           less of the value of WANTT.)

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       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 zlaqr1 (integer N, complex*16, dimension( ldh, * ) H, integer LDH, complex*16 S1,
       complex*16 S2, complex*16, dimension( * ) V)
       ZLAQR1 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, ZLAQR1 sets v to a
                 scalar multiple of the first column of the product

                 (*)  K = (H - s1*I)*(H - s2*I)

                 scaling to avoid overflows and most underflows.

                 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 COMPLEX*16 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 >= N

           S1

                     S1 is COMPLEX*16

           S2

                     S2 is COMPLEX*16

                     S1 and S2 are the shifts defining K in (*) above.

           V

                     V is COMPLEX*16 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.

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

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

       Purpose:

               ZLAQR2 is identical to ZLAQR3 except that it avoids
               recursion by calling ZLAHQR instead of ZLAQR4.

               Aggressive early deflation:

               ZLAQR2 accepts as input an upper Hessenberg matrix
               H and performs an unitary 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 unitary 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 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 unitary matrix Z is updated so
                     so that the unitary 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 unitary 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 <= NW <= (KBOT-KTOP+1).

           H

                     H is COMPLEX*16 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 a unitary
                     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 <= 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 <= ILOZ <= IHIZ <= N.

           Z

                     Z is COMPLEX*16 array, dimension (LDZ,N)
                     IF WANTZ is .TRUE., then on output, the unitary
                     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 <= 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.

           SH

                     SH is COMPLEX*16 array, dimension (KBOT)
                     On output, approximate eigenvalues that may
                     be used for shifts are stored in SH(KBOT-ND-NS+1)
                     through SR(KBOT-ND).  Converged eigenvalues are
                     stored in SH(KBOT-ND+1) through SH(KBOT).

           V

                     V is COMPLEX*16 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 <= LDV

           NH

                     NH is INTEGER
                     The number of columns of T.  NH >= NW.

           T

                     T is COMPLEX*16 array, dimension (LDT,NW)

           LDT

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

           NV

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

           WV

                     WV is COMPLEX*16 array, dimension (LDWV,NW)

           LDWV

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

           WORK

                     WORK is COMPLEX*16 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; ZLAQR2
                     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.

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

   subroutine zlaqr3 (logical WANTT, logical WANTZ, integer N, integer KTOP, integer KBOT,
       integer NW, complex*16, dimension( ldh, * ) H, integer LDH, integer ILOZ, integer IHIZ,
       complex*16, dimension( ldz, * ) Z, integer LDZ, integer NS, integer ND, complex*16,
       dimension( * ) SH, complex*16, dimension( ldv, * ) V, integer LDV, integer NH, complex*16,
       dimension( ldt, * ) T, integer LDT, integer NV, complex*16, dimension( ldwv, * ) WV,
       integer LDWV, complex*16, dimension( * ) WORK, integer LWORK)
       ZLAQR3 performs the unitary 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:

               ZLAQR3 accepts as input an upper Hessenberg matrix
               H and performs an unitary 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 unitary 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 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 unitary matrix Z is updated so
                     so that the unitary 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 unitary 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 <= NW <= (KBOT-KTOP+1).

           H

                     H is COMPLEX*16 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 a unitary
                     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 <= 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 <= ILOZ <= IHIZ <= N.

           Z

                     Z is COMPLEX*16 array, dimension (LDZ,N)
                     IF WANTZ is .TRUE., then on output, the unitary
                     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 <= 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.

           SH

                     SH is COMPLEX*16 array, dimension (KBOT)
                     On output, approximate eigenvalues that may
                     be used for shifts are stored in SH(KBOT-ND-NS+1)
                     through SR(KBOT-ND).  Converged eigenvalues are
                     stored in SH(KBOT-ND+1) through SH(KBOT).

           V

                     V is COMPLEX*16 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 <= LDV

           NH

                     NH is INTEGER
                     The number of columns of T.  NH >= NW.

           T

                     T is COMPLEX*16 array, dimension (LDT,NW)

           LDT

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

           NV

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

           WV

                     WV is COMPLEX*16 array, dimension (LDWV,NW)

           LDWV

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

           WORK

                     WORK is COMPLEX*16 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; ZLAQR3
                     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.

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

   subroutine zlaqr4 (logical WANTT, logical WANTZ, integer N, integer ILO, integer IHI,
       complex*16, dimension( ldh, * ) H, integer LDH, complex*16, dimension( * ) W, integer
       ILOZ, integer IHIZ, complex*16, dimension( ldz, * ) Z, integer LDZ, complex*16, dimension(
       * ) WORK, integer LWORK, integer INFO)
       ZLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from
       the Schur decomposition.

       Purpose:

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

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

               Optionally Z may be postmultiplied into an input unitary
               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 unitary matrix Q:  A = Q*H*Q**H = (QZ)*H*(QZ)**H.

       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 triangular in rows
                      and columns 1:ILO-1 and IHI+1:N and, if ILO > 1,
                      H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
                      previous call to ZGEBAL, and then passed to ZGEHRD when the
                      matrix output by ZGEBAL is reduced to Hessenberg form.
                      Otherwise, ILO and IHI should be set to 1 and N,
                      respectively.  If N > 0, then 1 <= ILO <= IHI <= N.
                      If N = 0, then ILO = 1 and IHI = 0.

           H

                     H is COMPLEX*16 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 triangular matrix T from the Schur
                      decomposition (the Schur form). If INFO = 0 and WANT is
                      .FALSE., then the contents of H are unspecified on exit.
                      (The output value of H when INFO > 0 is given under the
                      description of INFO below.)

                      This subroutine may explicitly set H(i,j) = 0 for i > 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 >= max(1,N).

           W

                     W is COMPLEX*16 array, dimension (N)
                      The computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored
                      in W(ILO:IHI). 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 W(i) = H(i,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 COMPLEX*16 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 > 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 >= MAX(1,IHIZ).  Otherwise, LDZ >= 1.

           WORK

                     WORK is COMPLEX*16 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 >= 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 ZLAQR4 does a workspace query.
                      In this case, ZLAQR4 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
                        > 0:  if INFO = i, ZLAQR4 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 > 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 > 0 and WANTT is .TRUE., then on exit

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

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

                           If INFO > 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 unitary matrix in (*) (regard-
                           less of the value of WANTT.)

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       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 zlaqr5 (logical WANTT, logical WANTZ, integer KACC22, integer N, integer KTOP,
       integer KBOT, integer NSHFTS, complex*16, dimension( * ) S, complex*16, dimension( ldh, *
       ) H, integer LDH, integer ILOZ, integer IHIZ, complex*16, dimension( ldz, * ) Z, integer
       LDZ, complex*16, dimension( ldv, * ) V, integer LDV, complex*16, dimension( ldu, * ) U,
       integer LDU, integer NV, complex*16, dimension( ldwv, * ) WV, integer LDWV, integer NH,
       complex*16, dimension( ldwh, * ) WH, integer LDWH)
       ZLAQR5 performs a single small-bulge multi-shift QR sweep.

       Purpose:

               ZLAQR5, called by ZLAQR0, performs a
               single small-bulge multi-shift QR sweep.

       Parameters
           WANTT

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

           WANTZ

                     WANTZ is LOGICAL
                        WANTZ = .true. if the unitary 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: ZLAQR5 does not accumulate reflections and does not
                        use matrix-matrix multiply to update far-from-diagonal
                        matrix entries.
                   = 1: ZLAQR5 accumulates reflections and uses matrix-matrix
                        multiply to update the far-from-diagonal matrix entries.
                   = 2: Same as KACC22 = 1. This option used to enable exploiting
                        the 2-by-2 structure during matrix multiplications, but
                        this is no longer supported.

           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.

           S

                     S is COMPLEX*16 array, dimension (NSHFTS)
                        S contains the shifts of origin that define the multi-
                        shift QR sweep.  On output S may be reordered.

           H

                     H is COMPLEX*16 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 >= 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 <= ILOZ <= IHIZ <= N

           Z

                     Z is COMPLEX*16 array, dimension (LDZ,IHIZ)
                        If WANTZ = .TRUE., then the QR Sweep unitary
                        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 >= N.

           V

                     V is COMPLEX*16 array, dimension (LDV,NSHFTS/2)

           LDV

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

           U

                     U is COMPLEX*16 array, dimension (LDU,2*NSHFTS)

           LDU

                     LDU is INTEGER
                        LDU is the leading dimension of U just as declared in the
                        in the calling subroutine.  LDU >= 2*NSHFTS.

           NV

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

           WV

                     WV is COMPLEX*16 array, dimension (LDWV,2*NSHFTS)

           LDWV

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

           NH

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

           WH

                     WH is COMPLEX*16 array, dimension (LDWH,NH)

           LDWH

                     LDWH is INTEGER
                        Leading dimension of WH just as declared in the
                        calling procedure.  LDWH >= 2*NSHFTS.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

       Lars Karlsson, Daniel Kressner, and Bruno Lang

       Thijs Steel, Department of Computer science, KU Leuven, Belgium

       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.

       Lars Karlsson, Daniel Kressner, and Bruno Lang, Optimally packed chains of bulges in
       multishift QR algorithms. ACM Trans. Math. Softw. 40, 2, Article 12 (February 2014).

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

       Purpose:

            ZLAQSB 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 COMPLEX*16 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**H *U or A = L*L**H 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.

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

       Purpose:

            ZLAQSP 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 COMPLEX*16 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.

   subroutine zlar1v (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, complex*16, 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)
       ZLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1
       through bn of the tridiagonal matrix LDLT - λI.

       Purpose:

            ZLAR1V 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 COMPLEX*16 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.

       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 zlar2v (integer N, complex*16, dimension( * ) X, complex*16, dimension( * ) Y,
       complex*16, dimension( * ) Z, integer INCX, double precision, dimension( * ) C,
       complex*16, dimension( * ) S, integer INCC)
       ZLAR2V applies a vector of plane rotations with real cosines and complex sines from both
       sides to a sequence of 2-by-2 symmetric/Hermitian matrices.

       Purpose:

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

               (       x(i)  z(i) ) :=
               ( conjg(z(i)) y(i) )

                 (  c(i) conjg(s(i)) ) (       x(i)  z(i) ) ( c(i) -conjg(s(i)) )
                 ( -s(i)       c(i)  ) ( conjg(z(i)) y(i) ) ( s(i)        c(i)  )

       Parameters
           N

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

           X

                     X is COMPLEX*16 array, dimension (1+(N-1)*INCX)
                     The vector x; the elements of x are assumed to be real.

           Y

                     Y is COMPLEX*16 array, dimension (1+(N-1)*INCX)
                     The vector y; the elements of y are assumed to be real.

           Z

                     Z is COMPLEX*16 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 COMPLEX*16 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.

   subroutine zlarcm (integer M, integer N, double precision, dimension( lda, * ) A, integer LDA,
       complex*16, dimension( ldb, * ) B, integer LDB, complex*16, dimension( ldc, * ) C, integer
       LDC, double precision, dimension( * ) RWORK)
       ZLARCM copies all or part of a real two-dimensional array to a complex array.

       Purpose:

            ZLARCM performs a very simple matrix-matrix multiplication:
                     C := A * B,
            where A is M by M and real; B is M by N and complex;
            C is M by N and complex.

       Parameters
           M

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

           N

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

           A

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

           LDA

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

           B

                     B is COMPLEX*16 array, dimension (LDB, N)
                     On entry, B contains the M by N matrix B.

           LDB

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

           C

                     C is COMPLEX*16 array, dimension (LDC, N)
                     On exit, C contains the M by N matrix C.

           LDC

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

           RWORK

                     RWORK is DOUBLE PRECISION array, dimension (2*M*N)

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine zlarf (character SIDE, integer M, integer N, complex*16, dimension( * ) V, integer
       INCV, complex*16 TAU, complex*16, dimension( ldc, * ) C, integer LDC, complex*16,
       dimension( * ) WORK)
       ZLARF applies an elementary reflector to a general rectangular matrix.

       Purpose:

            ZLARF applies a complex elementary reflector H to a complex M-by-N
            matrix C, from either the left or the right. H is represented in the
            form

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

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

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

            To apply H**H, supply conjg(tau) instead
            tau.

       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 COMPLEX*16 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 COMPLEX*16
                     The value tau in the representation of H.

           C

                     C is COMPLEX*16 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 COMPLEX*16 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.

   subroutine zlarfb (character SIDE, character TRANS, character DIRECT, character STOREV,
       integer M, integer N, integer K, complex*16, dimension( ldv, * ) V, integer LDV,
       complex*16, dimension( ldt, * ) T, integer LDT, complex*16, dimension( ldc, * ) C, integer
       LDC, complex*16, dimension( ldwork, * ) WORK, integer LDWORK)
       ZLARFB applies a block reflector or its conjugate-transpose to a general rectangular
       matrix.

       Purpose:

            ZLARFB applies a complex block reflector H or its transpose H**H to a
            complex M-by-N matrix C, from either the left or the right.

       Parameters
           SIDE

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

           TRANS

                     TRANS is CHARACTER*1
                     = 'N': apply H (No transpose)
                     = 'C': apply H**H (Conjugate 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).
                     If SIDE = 'L', M >= K >= 0;
                     if SIDE = 'R', N >= K >= 0.

           V

                     V is COMPLEX*16 array, dimension
                                           (LDV,K) if STOREV = 'C'
                                           (LDV,M) if STOREV = 'R' and SIDE = 'L'
                                           (LDV,N) if STOREV = 'R' and SIDE = 'R'
                     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 COMPLEX*16 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 COMPLEX*16 array, dimension (LDC,N)
                     On entry, the M-by-N matrix C.
                     On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H.

           LDC

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

           WORK

                     WORK is COMPLEX*16 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.

       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 zlarfb_gett (character IDENT, integer M, integer N, integer K, complex*16,
       dimension( ldt, * ) T, integer LDT, complex*16, dimension( lda, * ) A, integer LDA,
       complex*16, dimension( ldb, * ) B, integer LDB, complex*16, dimension( ldwork, * ) WORK,
       integer LDWORK)
       ZLARFB_GETT

       Purpose:

            ZLARFB_GETT applies a complex Householder block reflector H from the
            left to a complex (K+M)-by-N  'triangular-pentagonal' matrix
            composed of two block matrices: an upper trapezoidal K-by-N matrix A
            stored in the array A, and a rectangular M-by-(N-K) matrix B, stored
            in the array B. The block reflector H is stored in a compact
            WY-representation, where the elementary reflectors are in the
            arrays A, B and T. See Further Details section.

       Parameters
           IDENT

                     IDENT is CHARACTER*1
                     If IDENT = not 'I', or not 'i', then V1 is unit
                        lower-triangular and stored in the left K-by-K block of
                        the input matrix A,
                     If IDENT = 'I' or 'i', then  V1 is an identity matrix and
                        not stored.
                     See Further Details section.

           M

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

           N

                     N is INTEGER
                     The number of columns of the matrices A and B.
                     N >= 0.

           K

                     K is INTEGER
                     The number or rows of the matrix A.
                     K is also order of the matrix T, i.e. the number of
                     elementary reflectors whose product defines the block
                     reflector. 0 <= K <= N.

           T

                     T is COMPLEX*16 array, dimension (LDT,K)
                     The upper-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 COMPLEX*16 array, dimension (LDA,N)

                     On entry:
                      a) In the K-by-N upper-trapezoidal part A: input matrix A.
                      b) In the columns below the diagonal: columns of V1
                         (ones are not stored on the diagonal).

                     On exit:
                       A is overwritten by rectangular K-by-N product H*A.

                     See Further Details section.

           LDA

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

           B

                     B is COMPLEX*16 array, dimension (LDB,N)

                     On entry:
                       a) In the M-by-(N-K) right block: input matrix B.
                       b) In the M-by-N left block: columns of V2.

                     On exit:
                       B is overwritten by rectangular M-by-N product H*B.

                     See Further Details section.

           LDB

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

           WORK

                     WORK is COMPLEX*16 array,
                     dimension (LDWORK,max(K,N-K))

           LDWORK

                     LDWORK is INTEGER
                     The leading dimension of the array WORK. LDWORK>=max(1,K).

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Contributors:

            November 2020, Igor Kozachenko,
                           Computer Science Division,
                           University of California, Berkeley

       Further Details:

               (1) Description of the Algebraic Operation.

               The matrix A is a K-by-N matrix composed of two column block
               matrices, A1, which is K-by-K, and A2, which is K-by-(N-K):
               A = ( A1, A2 ).
               The matrix B is an M-by-N matrix composed of two column block
               matrices, B1, which is M-by-K, and B2, which is M-by-(N-K):
               B = ( B1, B2 ).

               Perform the operation:

                  ( A_out ) := H * ( A_in ) = ( I - V * T * V**H ) * ( A_in ) =
                  ( B_out )        ( B_in )                          ( B_in )
                             = ( I - ( V1 ) * T * ( V1**H, V2**H ) ) * ( A_in )
                                     ( V2 )                            ( B_in )
                On input:

               a) ( A_in )  consists of two block columns:
                  ( B_in )

                  ( A_in ) = (( A1_in ) ( A2_in )) = (( A1_in ) ( A2_in ))
                  ( B_in )   (( B1_in ) ( B2_in ))   ((     0 ) ( B2_in )),

                  where the column blocks are:

                  (  A1_in )  is a K-by-K upper-triangular matrix stored in the
                              upper triangular part of the array A(1:K,1:K).
                  (  B1_in )  is an M-by-K rectangular ZERO matrix and not stored.

                  ( A2_in )  is a K-by-(N-K) rectangular matrix stored
                             in the array A(1:K,K+1:N).
                  ( B2_in )  is an M-by-(N-K) rectangular matrix stored
                             in the array B(1:M,K+1:N).

               b) V = ( V1 )
                      ( V2 )

                  where:
                  1) if IDENT == 'I',V1 is a K-by-K identity matrix, not stored;
                  2) if IDENT != 'I',V1 is a K-by-K unit lower-triangular matrix,
                     stored in the lower-triangular part of the array
                     A(1:K,1:K) (ones are not stored),
                  and V2 is an M-by-K rectangular stored the array B(1:M,1:K),
                            (because on input B1_in is a rectangular zero
                             matrix that is not stored and the space is
                             used to store V2).

               c) T is a K-by-K upper-triangular matrix stored
                  in the array T(1:K,1:K).

               On output:

               a) ( A_out ) consists of two  block columns:
                  ( B_out )

                  ( A_out ) = (( A1_out ) ( A2_out ))
                  ( B_out )   (( B1_out ) ( B2_out )),

                  where the column blocks are:

                  ( A1_out )  is a K-by-K square matrix, or a K-by-K
                              upper-triangular matrix, if V1 is an
                              identity matrix. AiOut is stored in
                              the array A(1:K,1:K).
                  ( B1_out )  is an M-by-K rectangular matrix stored
                              in the array B(1:M,K:N).

                  ( A2_out )  is a K-by-(N-K) rectangular matrix stored
                              in the array A(1:K,K+1:N).
                  ( B2_out )  is an M-by-(N-K) rectangular matrix stored
                              in the array B(1:M,K+1:N).

               The operation above can be represented as the same operation
               on each block column:

                  ( A1_out ) := H * ( A1_in ) = ( I - V * T * V**H ) * ( A1_in )
                  ( B1_out )        (     0 )                          (     0 )

                  ( A2_out ) := H * ( A2_in ) = ( I - V * T * V**H ) * ( A2_in )
                  ( B2_out )        ( B2_in )                          ( B2_in )

               If IDENT != 'I':

                  The computation for column block 1:

                  A1_out: = A1_in - V1*T*(V1**H)*A1_in

                  B1_out: = - V2*T*(V1**H)*A1_in

                  The computation for column block 2, which exists if N > K:

                  A2_out: = A2_in - V1*T*( (V1**H)*A2_in + (V2**H)*B2_in )

                  B2_out: = B2_in - V2*T*( (V1**H)*A2_in + (V2**H)*B2_in )

               If IDENT == 'I':

                  The operation for column block 1:

                  A1_out: = A1_in - V1*T*A1_in

                  B1_out: = - V2*T*A1_in

                  The computation for column block 2, which exists if N > K:

                  A2_out: = A2_in - T*( A2_in + (V2**H)*B2_in )

                  B2_out: = B2_in - V2*T*( A2_in + (V2**H)*B2_in )

               (2) Description of the Algorithmic Computation.

               In the first step, we compute column block 2, i.e. A2 and B2.
               Here, we need to use the K-by-(N-K) rectangular workspace
               matrix W2 that is of the same size as the matrix A2.
               W2 is stored in the array WORK(1:K,1:(N-K)).

               In the second step, we compute column block 1, i.e. A1 and B1.
               Here, we need to use the K-by-K square workspace matrix W1
               that is of the same size as the as the matrix A1.
               W1 is stored in the array WORK(1:K,1:K).

               NOTE: Hence, in this routine, we need the workspace array WORK
               only of size WORK(1:K,1:max(K,N-K)) so it can hold both W2 from
               the first step and W1 from the second step.

               Case (A), when V1 is unit lower-triangular, i.e. IDENT != 'I',
               more computations than in the Case (B).

               if( IDENT != 'I' ) then
                if ( N > K ) then
                  (First Step - column block 2)
                  col2_(1) W2: = A2
                  col2_(2) W2: = (V1**H) * W2 = (unit_lower_tr_of_(A1)**H) * W2
                  col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2
                  col2_(4) W2: = T * W2
                  col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2
                  col2_(6) W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2
                  col2_(7) A2: = A2 - W2
                else
                  (Second Step - column block 1)
                  col1_(1) W1: = A1
                  col1_(2) W1: = (V1**H) * W1 = (unit_lower_tr_of_(A1)**H) * W1
                  col1_(3) W1: = T * W1
                  col1_(4) B1: = - V2 * W1 = - B1 * W1
                  col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1
                  col1_(6) square A1: = A1 - W1
                end if
               end if

               Case (B), when V1 is an identity matrix, i.e. IDENT == 'I',
               less computations than in the Case (A)

               if( IDENT == 'I' ) then
                if ( N > K ) then
                  (First Step - column block 2)
                  col2_(1) W2: = A2
                  col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2
                  col2_(4) W2: = T * W2
                  col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2
                  col2_(7) A2: = A2 - W2
                else
                  (Second Step - column block 1)
                  col1_(1) W1: = A1
                  col1_(3) W1: = T * W1
                  col1_(4) B1: = - V2 * W1 = - B1 * W1
                  col1_(6) upper-triangular_of_(A1): = A1 - W1
                end if
               end if

               Combine these cases (A) and (B) together, this is the resulting
               algorithm:

               if ( N > K ) then

                 (First Step - column block 2)

                 col2_(1)  W2: = A2
                 if( IDENT != 'I' ) then
                   col2_(2)  W2: = (V1**H) * W2
                                 = (unit_lower_tr_of_(A1)**H) * W2
                 end if
                 col2_(3)  W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2]
                 col2_(4)  W2: = T * W2
                 col2_(5)  B2: = B2 - V2 * W2 = B2 - B1 * W2
                 if( IDENT != 'I' ) then
                   col2_(6)    W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2
                 end if
                 col2_(7) A2: = A2 - W2

               else

               (Second Step - column block 1)

                 col1_(1) W1: = A1
                 if( IDENT != 'I' ) then
                   col1_(2) W1: = (V1**H) * W1
                               = (unit_lower_tr_of_(A1)**H) * W1
                 end if
                 col1_(3) W1: = T * W1
                 col1_(4) B1: = - V2 * W1 = - B1 * W1
                 if( IDENT != 'I' ) then
                   col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1
                   col1_(6_a) below_diag_of_(A1): =  - below_diag_of_(W1)
                 end if
                 col1_(6_b) up_tr_of_(A1): = up_tr_of_(A1) - up_tr_of_(W1)

               end if

   subroutine zlarfg (integer N, complex*16 ALPHA, complex*16, dimension( * ) X, integer INCX,
       complex*16 TAU)
       ZLARFG generates an elementary reflector (Householder matrix).

       Purpose:

            ZLARFG generates a complex elementary reflector H of order n, such
            that

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

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

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

            where tau is a complex scalar and v is a complex (n-1)-element
            vector. Note that H is not hermitian.

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

            Otherwise  1 <= real(tau) <= 2  and  abs(tau-1) <= 1 .

       Parameters
           N

                     N is INTEGER
                     The order of the elementary reflector.

           ALPHA

                     ALPHA is COMPLEX*16
                     On entry, the value alpha.
                     On exit, it is overwritten with the value beta.

           X

                     X is COMPLEX*16 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 COMPLEX*16
                     The value tau.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine zlarfgp (integer N, complex*16 ALPHA, complex*16, dimension( * ) X, integer INCX,
       complex*16 TAU)
       ZLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.

       Purpose:

            ZLARFGP generates a complex elementary reflector H of order n, such
            that

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

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

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

            where tau is a complex scalar and v is a complex (n-1)-element
            vector. Note that H is not hermitian.

            If the elements of x are all zero and alpha is real, 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 COMPLEX*16
                     On entry, the value alpha.
                     On exit, it is overwritten with the value beta.

           X

                     X is COMPLEX*16 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 COMPLEX*16
                     The value tau.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

       Purpose:

            ZLARFT forms the triangular factor T of a complex 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**H

            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**H * 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 COMPLEX*16 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 COMPLEX*16 array, dimension (K)
                     TAU(i) must contain the scalar factor of the elementary
                     reflector H(i).

           T

                     T is COMPLEX*16 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.

       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 zlarfx (character SIDE, integer M, integer N, complex*16, dimension( * ) V,
       complex*16 TAU, complex*16, dimension( ldc, * ) C, integer LDC, complex*16, dimension( * )
       WORK)
       ZLARFX applies an elementary reflector to a general rectangular matrix, with loop
       unrolling when the reflector has order ≤ 10.

       Purpose:

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

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

            where tau is a complex scalar and v is a complex 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 COMPLEX*16 array, dimension (M) if SIDE = 'L'
                                                   or (N) if SIDE = 'R'
                     The vector v in the representation of H.

           TAU

                     TAU is COMPLEX*16
                     The value tau in the representation of H.

           C

                     C is COMPLEX*16 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 COMPLEX*16 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.

   subroutine zlarfy (character UPLO, integer N, complex*16, dimension( * ) V, integer INCV,
       complex*16 TAU, complex*16, dimension( ldc, * ) C, integer LDC, complex*16, dimension( * )
       WORK)
       ZLARFY

       Purpose:

            ZLARFY applies an elementary reflector, or Householder matrix, H,
            to an n x n Hermitian matrix C, from both the left and the right.

            H is represented in the form

               H = I - tau * v * v'

            where  tau  is a scalar and  v  is a vector.

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

       Parameters
           UPLO

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

           N

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

           V

                     V is COMPLEX*16 array, dimension
                             (1 + (N-1)*abs(INCV))
                     The vector v as described above.

           INCV

                     INCV is INTEGER
                     The increment between successive elements of v.  INCV must
                     not be zero.

           TAU

                     TAU is COMPLEX*16
                     The value tau as described above.

           C

                     C is COMPLEX*16 array, dimension (LDC, N)
                     On entry, the matrix C.
                     On exit, C is overwritten by H * C * H'.

           LDC

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

           WORK

                     WORK is COMPLEX*16 array, dimension (N)

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

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

       Purpose:

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

               (        c(i)   s(i) ) ( x(i) ) = ( r(i) )
               ( -conjg(s(i))  c(i) ) ( y(i) ) = (   0  )

               where c(i)**2 + ABS(s(i))**2 = 1

            The following conventions are used (these are the same as in ZLARTG,
            but differ from the BLAS1 routine ZROTG):
               If y(i)=0, then c(i)=1 and s(i)=0.
               If x(i)=0, then c(i)=0 and s(i) is chosen so that r(i) is real.

       Parameters
           N

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

           X

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

           INCX

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

           Y

                     Y is COMPLEX*16 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.

       Further Details:

             6-6-96 - Modified with a new algorithm by W. Kahan and J. Demmel

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

   subroutine zlarnv (integer IDIST, integer, dimension( 4 ) ISEED, integer N, complex*16,
       dimension( * ) X)
       ZLARNV returns a vector of random numbers from a uniform or normal distribution.

       Purpose:

            ZLARNV returns a vector of n random complex numbers from a uniform or
            normal distribution.

       Parameters
           IDIST

                     IDIST is INTEGER
                     Specifies the distribution of the random numbers:
                     = 1:  real and imaginary parts each uniform (0,1)
                     = 2:  real and imaginary parts each uniform (-1,1)
                     = 3:  real and imaginary parts each normal (0,1)
                     = 4:  uniformly distributed on the disc abs(z) < 1
                     = 5:  uniformly distributed on the circle abs(z) = 1

           ISEED

                     ISEED is INTEGER array, dimension (4)
                     On entry, the seed of the random number generator; the array
                     elements must be between 0 and 4095, and ISEED(4) must be
                     odd.
                     On exit, the seed is updated.

           N

                     N is INTEGER
                     The number of random numbers to be generated.

           X

                     X is COMPLEX*16 array, dimension (N)
                     The generated random numbers.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Further Details:

             This routine calls the auxiliary routine DLARUV to generate random
             real numbers from a uniform (0,1) distribution, in batches of up to
             128 using vectorisable code. The Box-Muller method is used to
             transform numbers from a uniform to a normal distribution.

   subroutine zlarrv (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,
       complex*16, dimension( ldz, * ) Z, integer LDZ, integer, dimension( * ) ISUPPZ, double
       precision, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)
       ZLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the
       eigenvalues of L D LT.

       Purpose:

            ZLARRV 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. Needed to compute gaps on the left or right
                     end of the extremal eigenvalues in the desired RANGE.

           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 < 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 COMPLEX*16 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 ZLARRV.
                     < 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.

       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 zlartv (integer N, complex*16, dimension( * ) X, integer INCX, complex*16,
       dimension( * ) Y, integer INCY, double precision, dimension( * ) C, complex*16, dimension(
       * ) S, integer INCC)
       ZLARTV applies a vector of plane rotations with real cosines and complex sines to the
       elements of a pair of vectors.

       Purpose:

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

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

       Parameters
           N

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

           X

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

           INCX

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

           Y

                     Y is COMPLEX*16 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 COMPLEX*16 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.

   subroutine zlascl (character TYPE, integer KL, integer KU, double precision CFROM, double
       precision CTO, integer M, integer N, complex*16, dimension( lda, * ) A, integer LDA,
       integer INFO)
       ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.

       Purpose:

            ZLASCL multiplies the M by N complex matrix A by the real scalar
            CTO/CFROM.  This is done without over/underflow as long as the final
            result CTO*A(I,J)/CFROM does not over/underflow. TYPE specifies that
            A may be full, upper triangular, lower triangular, upper Hessenberg,
            or banded.

       Parameters
           TYPE

                     TYPE is CHARACTER*1
                     TYPE indices the storage type of the input matrix.
                     = 'G':  A is a full matrix.
                     = 'L':  A is a lower triangular matrix.
                     = 'U':  A is an upper triangular matrix.
                     = 'H':  A is an upper Hessenberg matrix.
                     = 'B':  A is a symmetric band matrix with lower bandwidth KL
                             and upper bandwidth KU and with the only the lower
                             half stored.
                     = 'Q':  A is a symmetric band matrix with lower bandwidth KL
                             and upper bandwidth KU and with the only the upper
                             half stored.
                     = 'Z':  A is a band matrix with lower bandwidth KL and upper
                             bandwidth KU. See ZGBTRF for storage details.

           KL

                     KL is INTEGER
                     The lower bandwidth of A.  Referenced only if TYPE = 'B',
                     'Q' or 'Z'.

           KU

                     KU is INTEGER
                     The upper bandwidth of A.  Referenced only if TYPE = 'B',
                     'Q' or 'Z'.

           CFROM

                     CFROM is DOUBLE PRECISION

           CTO

                     CTO is DOUBLE PRECISION

                     The matrix A is multiplied by CTO/CFROM. A(I,J) is computed
                     without over/underflow if the final result CTO*A(I,J)/CFROM
                     can be represented without over/underflow.  CFROM must be
                     nonzero.

           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.

           A

                     A is COMPLEX*16 array, dimension (LDA,N)
                     The matrix to be multiplied by CTO/CFROM.  See TYPE for the
                     storage type.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.
                     If TYPE = 'G', 'L', 'U', 'H', LDA >= max(1,M);
                        TYPE = 'B', LDA >= KL+1;
                        TYPE = 'Q', LDA >= KU+1;
                        TYPE = 'Z', LDA >= 2*KL+KU+1.

           INFO

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine zlaset (character UPLO, integer M, integer N, complex*16 ALPHA, complex*16 BETA,
       complex*16, dimension( lda, * ) A, integer LDA)
       ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to
       given values.

       Purpose:

            ZLASET initializes a 2-D array A to BETA on the diagonal and
            ALPHA on the offdiagonals.

       Parameters
           UPLO

                     UPLO is CHARACTER*1
                     Specifies the part of the matrix A to be set.
                     = 'U':      Upper triangular part is set. The lower triangle
                                 is unchanged.
                     = 'L':      Lower triangular part is set. The upper triangle
                                 is unchanged.
                     Otherwise:  All of the matrix A is set.

           M

                     M is INTEGER
                     On entry, M specifies the number of rows of A.

           N

                     N is INTEGER
                     On entry, N specifies the number of columns of A.

           ALPHA

                     ALPHA is COMPLEX*16
                     All the offdiagonal array elements are set to ALPHA.

           BETA

                     BETA is COMPLEX*16
                     All the diagonal array elements are set to BETA.

           A

                     A is COMPLEX*16 array, dimension (LDA,N)
                     On entry, the m by n matrix A.
                     On exit, A(i,j) = ALPHA, 1 <= i <= m, 1 <= j <= n, i.ne.j;
                              A(i,i) = BETA , 1 <= i <= min(m,n)

           LDA

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine zlasr (character SIDE, character PIVOT, character DIRECT, integer M, integer N,
       double precision, dimension( * ) C, double precision, dimension( * ) S, complex*16,
       dimension( lda, * ) A, integer LDA)
       ZLASR applies a sequence of plane rotations to a general rectangular matrix.

       Purpose:

            ZLASR applies a sequence of real plane rotations to a complex matrix
            A, from either the left or the right.

            When SIDE = 'L', the transformation takes the form

               A := P*A

            and when SIDE = 'R', the transformation takes the form

               A := A*P**T

            where P is an orthogonal matrix consisting of a sequence of z plane
            rotations, with z = M when SIDE = 'L' and z = N when SIDE = 'R',
            and P**T is the transpose of P.

            When DIRECT = 'F' (Forward sequence), then

               P = P(z-1) * ... * P(2) * P(1)

            and when DIRECT = 'B' (Backward sequence), then

               P = P(1) * P(2) * ... * P(z-1)

            where P(k) is a plane rotation matrix defined by the 2-by-2 rotation

               R(k) = (  c(k)  s(k) )
                    = ( -s(k)  c(k) ).

            When PIVOT = 'V' (Variable pivot), the rotation is performed
            for the plane (k,k+1), i.e., P(k) has the form

               P(k) = (  1                                            )
                      (       ...                                     )
                      (              1                                )
                      (                   c(k)  s(k)                  )
                      (                  -s(k)  c(k)                  )
                      (                                1              )
                      (                                     ...       )
                      (                                            1  )

            where R(k) appears as a rank-2 modification to the identity matrix in
            rows and columns k and k+1.

            When PIVOT = 'T' (Top pivot), the rotation is performed for the
            plane (1,k+1), so P(k) has the form

               P(k) = (  c(k)                    s(k)                 )
                      (         1                                     )
                      (              ...                              )
                      (                     1                         )
                      ( -s(k)                    c(k)                 )
                      (                                 1             )
                      (                                      ...      )
                      (                                             1 )

            where R(k) appears in rows and columns 1 and k+1.

            Similarly, when PIVOT = 'B' (Bottom pivot), the rotation is
            performed for the plane (k,z), giving P(k) the form

               P(k) = ( 1                                             )
                      (      ...                                      )
                      (             1                                 )
                      (                  c(k)                    s(k) )
                      (                         1                     )
                      (                              ...              )
                      (                                     1         )
                      (                 -s(k)                    c(k) )

            where R(k) appears in rows and columns k and z.  The rotations are
            performed without ever forming P(k) explicitly.

       Parameters
           SIDE

                     SIDE is CHARACTER*1
                     Specifies whether the plane rotation matrix P is applied to
                     A on the left or the right.
                     = 'L':  Left, compute A := P*A
                     = 'R':  Right, compute A:= A*P**T

           PIVOT

                     PIVOT is CHARACTER*1
                     Specifies the plane for which P(k) is a plane rotation
                     matrix.
                     = 'V':  Variable pivot, the plane (k,k+1)
                     = 'T':  Top pivot, the plane (1,k+1)
                     = 'B':  Bottom pivot, the plane (k,z)

           DIRECT

                     DIRECT is CHARACTER*1
                     Specifies whether P is a forward or backward sequence of
                     plane rotations.
                     = 'F':  Forward, P = P(z-1)*...*P(2)*P(1)
                     = 'B':  Backward, P = P(1)*P(2)*...*P(z-1)

           M

                     M is INTEGER
                     The number of rows of the matrix A.  If m <= 1, an immediate
                     return is effected.

           N

                     N is INTEGER
                     The number of columns of the matrix A.  If n <= 1, an
                     immediate return is effected.

           C

                     C is DOUBLE PRECISION array, dimension
                             (M-1) if SIDE = 'L'
                             (N-1) if SIDE = 'R'
                     The cosines c(k) of the plane rotations.

           S

                     S is DOUBLE PRECISION array, dimension
                             (M-1) if SIDE = 'L'
                             (N-1) if SIDE = 'R'
                     The sines s(k) of the plane rotations.  The 2-by-2 plane
                     rotation part of the matrix P(k), R(k), has the form
                     R(k) = (  c(k)  s(k) )
                            ( -s(k)  c(k) ).

           A

                     A is COMPLEX*16 array, dimension (LDA,N)
                     The M-by-N matrix A.  On exit, A is overwritten by P*A if
                     SIDE = 'R' or by A*P**T if SIDE = 'L'.

           LDA

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

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine zlaswp (integer N, complex*16, dimension( lda, * ) A, integer LDA, integer K1,
       integer K2, integer, dimension( * ) IPIV, integer INCX)
       ZLASWP performs a series of row interchanges on a general rectangular matrix.

       Purpose:

            ZLASWP 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 COMPLEX*16 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.

       Further Details:

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

   subroutine zlat2c (character UPLO, integer N, complex*16, dimension( lda, * ) A, integer LDA,
       complex, dimension( ldsa, * ) SA, integer LDSA, integer INFO)
       ZLAT2C converts a double complex triangular matrix to a complex triangular matrix.

       Purpose:

            ZLAT2C converts a COMPLEX*16 triangular matrix, SA, to a COMPLEX
            triangular matrix, A.

            RMAX is the overflow for the SINGLE PRECISION arithmetic
            ZLAT2C 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 COMPLEX*16 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 COMPLEX 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.

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

       Purpose:

            ZLATBS solves one of the triangular systems

               A * x = s*b,  A**T * x = s*b,  or  A**H * 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 ZTBSV 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**H * x = s*b  (Conjugate 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 COMPLEX*16 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 COMPLEX*16 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,  A**T * x = s*b,  or  A**H * 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.

       Further Details:

             A rough bound on x is computed; if that is less than overflow, ZTBSV
             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 ZTBSV 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  or
             A**H *x = b.  The basic algorithm for A upper triangular is

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

             We simultaneously compute two bounds
                  G(j) = bound on ( b(i) - A[1:i-1,i]' * 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 ZTBSV if 1/M(n) and 1/G(n) are both greater
             than max(underflow, 1/overflow).

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

       Purpose:

            ZLATDF computes the contribution to the reciprocal Dif-estimate
            by solving for x in Z * x = b, where b is chosen such that the norm
            of x is as large as possible. It is assumed that LU decomposition
            of Z has been computed by ZGETC2. On entry RHS = f holds the
            contribution from earlier solved sub-systems, and on return RHS = x.

            The factorization of Z returned by ZGETC2 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 ZGECON, 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 COMPLEX*16 array, dimension (LDZ, N)
                     On entry, the LU part of the factorization of the n-by-n
                     matrix Z computed by ZGETC2:  Z = P * L * U * Q

           LDZ

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

           RHS

                     RHS is COMPLEX*16 array, dimension (N).
                     On entry, RHS contains contributions from other subsystems.
                     On exit, RHS contains the solution of the subsystem with
                     entries according 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 ZTGSYL, 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 ZTGSY2 is called by CTGSYL.

           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 ZTGSY2 is called by
                     ZTGSYL.

           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.

       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 UMINF-95.05,
           Department of Computing Science, Umea University, S-901 87 Umea, Sweden, 1995.

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

       Purpose:

            ZLATPS solves one of the triangular systems

               A * x = s*b,  A**T * x = s*b,  or  A**H * 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, A**H denotes the conjugate 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 ZTPSV 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**H * x = s*b  (Conjugate 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 COMPLEX*16 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 COMPLEX*16 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,  A**T * x = s*b,  or  A**H * 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.

       Further Details:

             A rough bound on x is computed; if that is less than overflow, ZTPSV
             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 ZTPSV 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  or
             A**H *x = b.  The basic algorithm for A upper triangular is

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

             We simultaneously compute two bounds
                  G(j) = bound on ( b(i) - A[1:i-1,i]' * 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 ZTPSV if 1/M(n) and 1/G(n) are both greater
             than max(underflow, 1/overflow).

   subroutine zlatrd (character UPLO, integer N, integer NB, complex*16, dimension( lda, * ) A,
       integer LDA, double precision, dimension( * ) E, complex*16, dimension( * ) TAU,
       complex*16, dimension( ldw, * ) W, integer LDW)
       ZLATRD reduces the first nb rows and columns of a symmetric/Hermitian matrix A to real
       tridiagonal form by an unitary similarity transformation.

       Purpose:

            ZLATRD reduces NB rows and columns of a complex Hermitian matrix A to
            Hermitian tridiagonal form by a unitary similarity
            transformation Q**H * 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', ZLATRD reduces the last NB rows and columns of a
            matrix, of which the upper triangle is supplied;
            if UPLO = 'L', ZLATRD reduces the first NB rows and columns of a
            matrix, of which the lower triangle is supplied.

            This is an auxiliary routine called by ZHETRD.

       Parameters
           UPLO

                     UPLO is CHARACTER*1
                     Specifies whether the upper or lower triangular part of the
                     Hermitian 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 COMPLEX*16 array, dimension (LDA,N)
                     On entry, the Hermitian 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 unitary 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  unitary matrix Q as a
                       product of elementary reflectors.
                     See Further Details.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(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 COMPLEX*16 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 COMPLEX*16 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.

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

             where tau is a complex scalar, and v is a complex 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**H

             where tau is a complex scalar, and v is a complex 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 Hermitian rank-2k update of the form:
             A := A - V*W**H - W*V**H.

             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 zlatrs (character UPLO, character TRANS, character DIAG, character NORMIN, integer
       N, complex*16, dimension( lda, * ) A, integer LDA, complex*16, dimension( * ) X, double
       precision SCALE, double precision, dimension( * ) CNORM, integer INFO)
       ZLATRS solves a triangular system of equations with the scale factor set to prevent
       overflow.

       Purpose:

            ZLATRS solves one of the triangular systems

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

            with scaling to prevent overflow.  Here A is an upper or lower
            triangular matrix, A**T denotes the transpose of A, A**H denotes the
            conjugate 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
            ZTRSV 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**H * x = s*b  (Conjugate 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 COMPLEX*16 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 COMPLEX*16 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,  A**T * x = s*b,  or  A**H * 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.

       Further Details:

             A rough bound on x is computed; if that is less than overflow, ZTRSV
             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 ZTRSV 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  or
             A**H *x = b.  The basic algorithm for A upper triangular is

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

             We simultaneously compute two bounds
                  G(j) = bound on ( b(i) - A[1:i-1,i]' * 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 ZTRSV if 1/M(n) and 1/G(n) are both greater
             than max(underflow, 1/overflow).

   subroutine zlauu2 (character UPLO, integer N, complex*16, dimension( lda, * ) A, integer LDA,
       integer INFO)
       ZLAUU2 computes the product UUH or LHL, where U and L are upper or lower triangular
       matrices (unblocked algorithm).

       Purpose:

            ZLAUU2 computes the product U * U**H or L**H * 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 COMPLEX*16 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**H;
                     if UPLO = 'L', the lower triangle of A is overwritten with
                     the lower triangle of the product L**H * 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.

   subroutine zlauum (character UPLO, integer N, complex*16, dimension( lda, * ) A, integer LDA,
       integer INFO)
       ZLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular
       matrices (blocked algorithm).

       Purpose:

            ZLAUUM computes the product U * U**H or L**H * 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 COMPLEX*16 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**H;
                     if UPLO = 'L', the lower triangle of A is overwritten with
                     the lower triangle of the product L**H * 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.

   subroutine zrot (integer N, complex*16, dimension( * ) CX, integer INCX, complex*16,
       dimension( * ) CY, integer INCY, double precision C, complex*16 S)
       ZROT applies a plane rotation with real cosine and complex sine to a pair of complex
       vectors.

       Purpose:

            ZROT   applies a plane rotation, where the cos (C) is real and the
            sin (S) is complex, and the vectors CX and CY are complex.

       Parameters
           N

                     N is INTEGER
                     The number of elements in the vectors CX and CY.

           CX

                     CX is COMPLEX*16 array, dimension (N)
                     On input, the vector X.
                     On output, CX is overwritten with C*X + S*Y.

           INCX

                     INCX is INTEGER
                     The increment between successive values of CX.  INCX <> 0.

           CY

                     CY is COMPLEX*16 array, dimension (N)
                     On input, the vector Y.
                     On output, CY is overwritten with -CONJG(S)*X + C*Y.

           INCY

                     INCY is INTEGER
                     The increment between successive values of CY.  INCX <> 0.

           C

                     C is DOUBLE PRECISION

           S

                     S is COMPLEX*16
                     C and S define a rotation
                        [  C          S  ]
                        [ -conjg(S)   C  ]
                     where C*C + S*CONJG(S) = 1.0.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine zspmv (character UPLO, integer N, complex*16 ALPHA, complex*16, dimension( * ) AP,
       complex*16, dimension( * ) X, integer INCX, complex*16 BETA, complex*16, dimension( * ) Y,
       integer INCY)
       ZSPMV computes a matrix-vector product for complex vectors using a complex symmetric
       packed matrix

       Purpose:

            ZSPMV  performs the matrix-vector operation

               y := alpha*A*x + beta*y,

            where alpha and beta are scalars, x and y are n element vectors and
            A is an n by n symmetric matrix, supplied in packed form.

       Parameters
           UPLO

                     UPLO is CHARACTER*1
                      On entry, UPLO specifies whether the upper or lower
                      triangular part of the matrix A is supplied in the packed
                      array AP as follows:

                         UPLO = 'U' or 'u'   The upper triangular part of A is
                                             supplied in AP.

                         UPLO = 'L' or 'l'   The lower triangular part of A is
                                             supplied in AP.

                      Unchanged on exit.

           N

                     N is INTEGER
                      On entry, N specifies the order of the matrix A.
                      N must be at least zero.
                      Unchanged on exit.

           ALPHA

                     ALPHA is COMPLEX*16
                      On entry, ALPHA specifies the scalar alpha.
                      Unchanged on exit.

           AP

                     AP is COMPLEX*16 array, dimension at least
                      ( ( N*( N + 1 ) )/2 ).
                      Before entry, with UPLO = 'U' or 'u', the array AP must
                      contain the upper triangular part of the symmetric matrix
                      packed sequentially, column by column, so that AP( 1 )
                      contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
                      and a( 2, 2 ) respectively, and so on.
                      Before entry, with UPLO = 'L' or 'l', the array AP must
                      contain the lower triangular part of the symmetric matrix
                      packed sequentially, column by column, so that AP( 1 )
                      contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
                      and a( 3, 1 ) respectively, and so on.
                      Unchanged on exit.

           X

                     X is COMPLEX*16 array, dimension at least
                      ( 1 + ( N - 1 )*abs( INCX ) ).
                      Before entry, the incremented array X must contain the N-
                      element vector x.
                      Unchanged on exit.

           INCX

                     INCX is INTEGER
                      On entry, INCX specifies the increment for the elements of
                      X. INCX must not be zero.
                      Unchanged on exit.

           BETA

                     BETA is COMPLEX*16
                      On entry, BETA specifies the scalar beta. When BETA is
                      supplied as zero then Y need not be set on input.
                      Unchanged on exit.

           Y

                     Y is COMPLEX*16 array, dimension at least
                      ( 1 + ( N - 1 )*abs( INCY ) ).
                      Before entry, the incremented array Y must contain the n
                      element vector y. On exit, Y is overwritten by the updated
                      vector y.

           INCY

                     INCY is INTEGER
                      On entry, INCY specifies the increment for the elements of
                      Y. INCY must not be zero.
                      Unchanged on exit.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine zspr (character UPLO, integer N, complex*16 ALPHA, complex*16, dimension( * ) X,
       integer INCX, complex*16, dimension( * ) AP)
       ZSPR performs the symmetrical rank-1 update of a complex symmetric packed matrix.

       Purpose:

            ZSPR    performs the symmetric rank 1 operation

               A := alpha*x*x**H + A,

            where alpha is a complex scalar, x is an n element vector and A is an
            n by n symmetric matrix, supplied in packed form.

       Parameters
           UPLO

                     UPLO is CHARACTER*1
                      On entry, UPLO specifies whether the upper or lower
                      triangular part of the matrix A is supplied in the packed
                      array AP as follows:

                         UPLO = 'U' or 'u'   The upper triangular part of A is
                                             supplied in AP.

                         UPLO = 'L' or 'l'   The lower triangular part of A is
                                             supplied in AP.

                      Unchanged on exit.

           N

                     N is INTEGER
                      On entry, N specifies the order of the matrix A.
                      N must be at least zero.
                      Unchanged on exit.

           ALPHA

                     ALPHA is COMPLEX*16
                      On entry, ALPHA specifies the scalar alpha.
                      Unchanged on exit.

           X

                     X is COMPLEX*16 array, dimension at least
                      ( 1 + ( N - 1 )*abs( INCX ) ).
                      Before entry, the incremented array X must contain the N-
                      element vector x.
                      Unchanged on exit.

           INCX

                     INCX is INTEGER
                      On entry, INCX specifies the increment for the elements of
                      X. INCX must not be zero.
                      Unchanged on exit.

           AP

                     AP is COMPLEX*16 array, dimension at least
                      ( ( N*( N + 1 ) )/2 ).
                      Before entry, with  UPLO = 'U' or 'u', the array AP must
                      contain the upper triangular part of the symmetric matrix
                      packed sequentially, column by column, so that AP( 1 )
                      contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
                      and a( 2, 2 ) respectively, and so on. On exit, the array
                      AP is overwritten by the upper triangular part of the
                      updated matrix.
                      Before entry, with UPLO = 'L' or 'l', the array AP must
                      contain the lower triangular part of the symmetric matrix
                      packed sequentially, column by column, so that AP( 1 )
                      contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
                      and a( 3, 1 ) respectively, and so on. On exit, the array
                      AP is overwritten by the lower triangular part of the
                      updated matrix.
                      Note that the imaginary parts of the diagonal elements need
                      not be set, they are assumed to be zero, and on exit they
                      are set to zero.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

   subroutine ztprfb (character SIDE, character TRANS, character DIRECT, character STOREV,
       integer M, integer N, integer K, integer L, complex*16, dimension( ldv, * ) V, integer
       LDV, complex*16, dimension( ldt, * ) T, integer LDT, complex*16, dimension( lda, * ) A,
       integer LDA, complex*16, dimension( ldb, * ) B, integer LDB, complex*16, dimension(
       ldwork, * ) WORK, integer LDWORK)
       ZTPRFB applies a complex 'triangular-pentagonal' block reflector to a complex matrix,
       which is composed of two blocks.

       Purpose:

            ZTPRFB applies a complex 'triangular-pentagonal' block reflector H or its
            conjugate transpose H**H to a complex 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**H from the Left
                     = 'R': apply H or H**H from the Right

           TRANS

                     TRANS is CHARACTER*1
                     = 'N': apply H (No transpose)
                     = 'C': apply H**H (Conjugate 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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**H*C or C*H or C*H**H.  See Further Details.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.
                     If SIDE = 'L', LDA >= max(1,K);
                     If SIDE = 'R', LDA >= max(1,M).

           B

                     B is COMPLEX*16 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**H*C or C*H or C*H**H.  See Further Details.

           LDB

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

           WORK

                     WORK is COMPLEX*16 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.

       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.

Author

       Generated automatically by Doxygen for LAPACK from the source code.